An experimental study on nonlinear wave dynamics for freak waves over an uneven bottom

The effect of the non-uniform bathymetry on the nonlinear wave dynamics for the freak wave is investigated experimentally, with emphasis on the interrelations between different nonlinear behaviors resulting from various geometric parameters and spectral analysis. Both the frequency modulation and the nonlinear phase coupling can be provoked by the decreasing water depth and weakened after the top peak of the bar, the nonlinear exhibition for transferring energy to high-frequency contents over shoal supports that frequency modulation can reflect nonlinear phase coupling well. The consistent change of the Hilbert energy spectrum and the bicoherence shows that the main nonlinear interaction in the process of wave propagation in shallowing water is quadratic nonlinearity. In addition, the geometric study is conducted to investigate the effect of the water depth on the parametric variations, the research results show that the mean asymmetry and kurtosis change abruptly when the wave approaches the top peak of the bar. As the wave propagates along the water flume, freak waves can be generated at various locations, however, they appear more frequently as waves propagate close to the shallowest water depth, and the maximum probability of occurrence for a freak wave can be up to about 1%.


Introduction
Freak waves pose a significant threat to human lives and structures in open seas as well as coastline (Kharif and Pelinovsky, 2003;Teng and Ning, 2009;Zhao et al., 2010;Ma et al., 2013;Li et al., 2015;Fu et al., 2021;Zhang and Benoit, 2021). There are different definitions for the freak wave, of which, wave heights are twice the significant wave height, or crest heights exceed 1.25 times the significant wave height often the common ones (Dysthe et al., 2008;Trulsen et al., 2012). Generally, when waves propagate over a shoal unidirectionally, the modulation instability becomes weaker (Shukla et al., 2006) as the water depth decreases, and when the relative water depth kh is smaller than 1.36, k is the wave number and h is the water depth, waves become stable (Baldock and Swan, 1996;Kharif and Pelinovsky, 2003;Whittaker et al., 2016). Therefore, when waves propagate toward the coast with varying topography, the mechanism for generating the freak wave could be different. As waves evolve from deeper to shallower water, the shallow water effect may be dominant, and the wave shape, wave spectrum, wave height, etc. will change a lot; in the meanwhile, decreasing wavelength and increasing wave steepness appear, and the freak wave may be provoked by shoaling effect (Doering and Bowen, 1995;Trulsen et al., 2012;Gramstad et al., 2013;Trulsen et al., 2020), studying the nonlinear characteristics and shallow effect which are related to the large waves deeply are of great significance to coastal and harbor engineering (Christou et al., 2008;Sergeeva et al., 2011;Zou and Peng, 2011;Mahmoudof et al., 2016;Ma et al., 2017;Gao et al., 2020;Mahmoudof and Hajivalie, 2021;Zhang and Benoit, 2021).
Some researchers have investigated the dependence of the geometric parameters on the nonlinear parameter Ursell number as waves propagated over shallow water (Pelinovsky and Sergeeva, 2006;Peng et al., 2009). Under different slopes, Chen et al. (2018a) analyzed the same relationship as Peng et al. (2009)'s; and in the meanwhile, Chen et al. (2018a) also showed that the skewness is larger in the shallower water depth than that in deep water. When waves propagate over variable water depths, the nonlinear interactions are attributed to not only the geometric parameters but the energy transformation or spectral shape. Beji and Battjes (1993) studied the spectral changes during the shoaling process experimentally and elucidated the generation of high-frequency energy in the process. The variation of the spectral shape depends on the comprehensive nonlinear interactions. When waves evolve over the uneven bathymetry from deeper to shallower water depth, the quartet wave interactions weaken and the triad wave interactions are dominant gradually, in the meanwhile, the energy of harmonics is also enhanced.
Bispectrum analysis is suited to study the nonlinear behaviors of shoaling waves (Elgar and Guza, 1985), and some researchers have studied the local triad wave interactions by the bispectrum method (Elgar and Guza, 1985;Doering and Bowen, 1995;Crawford and Hay, 2003;Christou et al., 2008). Dong et al. (2008) investigated the nonlinear phase coupling as a wave propagating over the curvilinear bottom, including the effect of wave breaking on the nonlinear phase coupling relationship, they elucidated that the degree of the quadratic phase coupling enhances as the water depth decreases. Cherneva and Soares (2007) analyzed the storms in the field by bispectrum and showed that components in a broad band take part in the second-order nonlinear interactions of the freak wave. Ma et al. (2014) carried out the bicoherence analysis for the segment of the wave series containing the freak wave and found that the higher harmonics are induced by the local wave-wave interactions which are exhibited by the obvious variations of the triad wave interactions and the energy transformations. Chen et al. (2018b) employed the bispectrum to investigate the effect of the bottom slope on the nonlinear triad wave interactions, and they found that a steeper slope can increase the interactions.
When the water depth varies from deeper to shallower (but not the extremely shallow water depth), the nonlinear effects including the resonance nonlinear interactions are needed to reveal the comprehensive mechanisms for the nonlinear wave dynamics. Bispectrum provides a good result for the analysis of the triad wave interactions when the water depth decreases, but the whole nonlinear interactions are hard to be exhibited by this method. Therefore, for the non-stationary and non-linear water waves, some researchers have analyzed the nonlinear interactions in the process of waves propagating by the Hilbert-Huang transform (HHT) (Huang et al., 1998). Based on HHT, Veltcheva (2001) analyzed the characteristics of wave groupiness and wave energy for the nearshore field data, the results showed that the small difference between intrinsic mode functions could be considered a necessary condition for the existence of wave groups, in the meanwhile, they verified that the shoaling could change the wave group structure. Wu and Yao (2004) provided that a smaller interwave frequency modulation can indicate the freak wave occurrence with a larger limiting wave steepness. Veltcheva and Soares (2016) investigated the intrawave frequency modulation in the vicinity of the extreme wave for the field data, they found that the intrawave frequency modulation increases as the asymmetry of the wave grows. The above studies verified that frequency modulation can indicate nonlinearity appropriately. Based on the dependent relationship between the frequency modulation and the nonlinearity, He et al. (2022) investigated and proposed a breaking criterion based on the intrawave frequency modulation, in the meanwhile, they also demonstrated the strong dependence of the intrawave frequency modulation on the local wave steepness.
Up to now, HHT is employed to investigate the nonlinear characteristics of waves evolving in the constant water depth mainly, for waves propagating over the uneven bathymetry, the wavelet-based bicoherence is frequently used, but HHT is used less. Considering the good applicability of HHT, in the present study, the HHT based on the ensemble empirical mode function (EEMD) (Wu and Huang, 2008) and the Hilbert transform (HT) and the wavelet-based bicoherence are combined and employed to investigate the nonlinear characteristics of the freak wave propagating over the shoaling bathymetry.
The paper is organized as follows: In Section 2, the introduction to the data processing methods is described, and the experimental set-up and wave conditions are introduced. The study results are presented in Section 3. Finally, the discussion and conclusion are delivered in Section 4.

Hilbert-Huang transform based on EEMD
According to Huang et al. (1998), the HHT consists of empirical mode decomposition (EMD) and the HT. To improve the intermittency problem and avoid mode fixing, the EMD is replaced by the EEMD in the front part of the HHT in the present study. The EEMD procedure is introduced briefly as follows based on the EMD. He et al. 10.3389/fmars.2023.1150896 Frontiers in Marine Science frontiersin.org For recorded data x(t), it can be decomposed into some intrinsic mode functions (IMFs) and a residual R by the EMD (Huang et al., 1998), the decomposed results can be denoted as follows where c j (t) is the j th IMF. For each IMF, like the j th IMF, the Hilbert transform, y j (t), can be obtained, as where P denotes the Cauchy principal value. With this definition, an analytic signal Z(t) can be obtained based on the complex conjugate pair c j (t) and y j (t), as Of which, a j (t) and q j (t) are the instantaneous amplitude and the instantaneous phase, respectively. Based on the instantaneous phase, the instantaneous frequency w j (t) can be defined as Finally, the recorded data x(t) can be described as The equation (6) represents the amplitude and the frequency as functions of time in the form of three-dimensional drawings, which is the Hilbert amplitude spectrum. When the amplitude is replaced by the squared amplitude, the Hilbert energy spectrum is produced.
When the EEMD is used, many so-called "artificial" signals can be comprised by adding random white noise into x(t) many times, and the k th "artificial" signal can be obtained, as in which, n k (t) is the white noise added. After applying the EMD procedure to x k (t), x k (t) can be denoted as Repeating steps (7) and (8) NE (Wu and Huang (2008) suggested a value of 300) times but for different white noise each time, based on EEMD, the j th IMF component c j (t) can be presented as Following the equations (2) through (6), the Hilbert energy spectrum based on the EEMD and HT is obtained. The amplitude of the white noise added is about 0.2 standard deviation of that of the recorded data (Wu and Huang, 2008).

Wavelet-based bicoherence
For recorded data x(t), the continuous wavelet transform W(a, t) is defined as follows (Dong et al., 2008;Dong et al., 2014;Ma et al., 2017;Chen et al., 2018b): in which, the (*) denotes the complex conjugate, and y a , t (t) represents a family of functions called wavelets. These functions are constructed by translating in timet and dilation with scale a (the scale a is the reciprocal of frequency, namely, f = 1/a.) the "mother" wavelet function y (t) (Dong et al., 2008) which is assumed by Morlet wavelet in the present study. Therefore, y a , t (t) is defined as To study the triad phase coupling interactions (Ma et al., 2017), the wavelet-based bispectrum is defined as where the (*) denotes the complex conjugate and T is the time of duration, and f 1 , f 2, and f must satisfy the frequency relationship: Usually, based on the bispectrum, the normalized squared wavelet bicoherence b 2 can be obtained to measure the degree of the nonlinear phase coupling interactions (Dong et al., 2008) b The value of b 2 is between 0 and 1, when b 2 is equal to 0, the coupling relationship between triads of waves is random phase; when b 2 is equal to 1, the amount of the coupling relationship is the maximum. To ensure the two-dimensional condition wave field, the flume is divided into two parts longitudinally with widths of 0.8 m and 2.2 m by a steel wall, and the part with a width of 0.8 m is selected as the working section. The wavemaker and the absorbing device are arranged at both ends of the flume, respectively. In the experiment, the average position of the wave paddle is set to 0.0 m, and the propagation direction of waves denoted by the capital letter C is marked in Figure 1A. Wave gauges are installed at 17 different locations, the data collected by the first wave gauge which is located 3.0 m away from the zero-point location is used to obtain the initial parameters, the installation locations of other wave gauges are sketched in Figure 1. Prior experiments have shown that the reflection factor is less than 5%, therefore, the reflection effect in the present study is ignored.

Experimental setup and wave conditions
In the present study, a set of irregular waves are generated based on the JONSWAP spectrum (Goda, 1999) with variable significant wave heights H s and peak periods T p , the peak enhancement factor is 3.3. In the experiment, the sample frequency f s is 100 Hz, and the time duration T is 200 s. To conduct statistical analysis better, twenty realizations with different random phases for each initial spectrum are generated. Therefore, 140 data realizations are obtained. In addition, a time of about 15 minutes is needed to guarantee a still water surface between two realizations. Four peak periods are used in the tests, then the value of k p h ranges from 1 to 3, which satisfies the finite water depth condition. Table 1 presents the details of the conditions for experimental cases. In which, f p is the peak frequency; k p corresponds to the wavenumber for f p and H max is the maximal wave height.

Results
Evolutions of wave parameters over the smooth breakwater Water depth has a significant effect on the wave profile, Figure 2 presents the evolutions of the wave surface elevations of case RJ 7 . When the wave propagates from the 0.45-m water depth to the 0.15-m water depth, sharper crests and flatter troughs make the wave profiles asymmetric. The main reasons for the transformation of the wave surface elevation are the nonlinear interactions as waves evolve over the shoaling bottom. It can be seen that the deeply asymmetric wave shapes still sustain after waves break, including the region of the whole top crest of the breakwater. When waves propagate to the de-shoaling region, the water depth increases gradually, and more symmetric wave profiles can be observed ( Figure 2H). Comparing the wave profiles of the time series in the same water depths but for shoaling and de-shoaling regions, respectively, it can be seen that when waves evolve close to the top crest of the breakwater, more asymmetric wave shapes and larger crests in Figure 2D than that in Figure 2H are found. The main reason for the difference in the time series results from the quadratic nonlinear interaction and wave breaking induced by the water depth, which also affects significantly the other parameters for water waves.
The nonlinear interactions are affected significantly by the topography for waves with different wavelengths and periods, the Ursell number, U r , is always used to denote the nonlinearity and where L is the wavelength based on the mean period. Figure 3 presents the relationship between the maximum Ursell number U rmax , and the maximum crest height h m , the dashed lines in the figure are the fitting curves, and the fitted equations are marked in the figure. As shown in Figure 3A, for the cases with different significant wave heights, the maximum crest heights are dependent on the maximum Ursell number linearly, in the meanwhile, the slopes for the growth curves are very close. This illustrates that the maximum Ursell numbers can well reflect the variations of the maximum crest heights with water depth. Figure 3B describes the relationship between the amplification factors AU rmax and Ah m . The amplification factor AU rmax can be obtained by the maximum Ursell number along the path dividing the one at the initial position; and the amplification factor Ah m can be obtained by the maximum crest height along the path dividing the one at the initial position. It can be seen from Figure 3B that there is an exponential dependence between the two parameters, as the maximum U r amplifies gradually, the enlargement of the maximum crest height weakens. The declined dependence for the two expansion factors indicates that the amplification proportion for the maximum crest heights is smaller than that of the maximum Ursell number. This phenomenon reflects that the maximum Ursell number still increases rapidly when the crest height grows moderately, which indicates that when the Ursell number has a Scatters for the relationship between the maximum Ursell number, U rmax , and the maximum crest height, h m (A), and the dependence between the enlargement factors AU rmax and Ah m for the U rmax and the h m (B). He et al. 10.3389/fmars.2023.1150896 Frontiers in Marine Science frontiersin.org large increase with the variable topography, the crest height may grow a little. It is shown that when the wave with a longer wavelength propagates over a shallower water depth, the Ursell number may become bigger according to the equation (15), but the crest height may not grow significantly in the case of weak nonlinearity. For the given imposed spectrum, the evolutions of the mean significant wave heights are shown in Figure 4. It can be seen that the significant wave heights increase up to the maximum as waves approach close to or on the top crest of the breakwater, and then start to decline sharply. After waves propagate across the deshoaling region, the significant wave heights are nearly the same for cases with different wave periods in Figure 4B except for the case RJ 1 whose wave period is 0.8 s. And for this case, a continuing downward trend is found as waves evolve along the wave flume. In the present experiment, the k p h for the case RJ 1 is about 3, and that for the others ranges from 1 to 2 ( Table 1). The spatial variations of the significant wave heights illustrate that the increasing nonlinear interactions attribute to more growth of the significant wave height for cases with longer wavelengths than that for cases with shorter wavelengths. When the wavelength is longer and the water depth is smaller, the triad wave interactions and bound wave interactions increase, and the growing nonlinearity makes the peaks of the significant wave heights. Figure 5 presents the spatial variations of the mean periods, T m , which is the mean value of the periods for all waves in a time series. It is shown that waves with longer periods are affected more obviously by uneven bathymetry. For waves with longer periods, there are two growing trends, both in the shoaling and the deshoaling regions, respectively. The local changes in wave parameters are dependent on the nonlinear dynamic interactions induced by the variable water depths (Zhang and Benoit, 2021). This phenomenon exhibits that the transition variations of the water depth contribute to the local growth in wavelength.
Nonlinearity leads to a departure of the wave field from the Gaussian distribution (Annenkov and Shrira, 2014), which is demonstrated by the higher moments of the wave field, such as skewness and kurtosis, the analyzed results for the two parameters are shown in Figures 6, 7. When waves evolve close to and on the top crest of the breakwater, the mean kurtosis and the mean skewness can achieve their local maxima, the result is the same as that of Ma et al. (2014) and Trulsen et al. (2020). In Chen et al. (2018a)'s study, three slopes (1/15, 1/30, and 1/45) are considered, and the value of the skewness ranges from 0 to about 1.5. The maximum skewness at a slope of 1/20 in the present study is close to that of Chen et al. (2018a). In addition, there are the minima around the end edge of the up-slope and the maxima around the start edge of the down-slope for the mean asymmetry, respectively. The results indicate that the wave profiles exhibit more asymmetric characteristics when the water depth is smaller and/or waves evolve across an abrupt transition of water depth.
For the mean kurtosis, it looks like there is a local minimum around the shallowest part of the lee side of the breakwater at x = 22 m (11 # ) as shown in Figure 6B, which is different from that of Trulsen et al. (2020). The main reason is that the kurtosis and skewness exhibit a squared dependence (Mori and Kobayashi, 1998;Ma et al., 2014), when the skewness is positive, the kurtosis grow absolutely by the skewness; however, when the skewness is negative, the kurtosis still keeps positive. In the present study, the mean skewness is negative at x = 24 m (13 # ), in the meanwhile, a maximum for the mean kurtosis appears, therefore, a local minimum appears between the two local maxima shown at x = 18.8 m (9 # ) and x = 24 m (13 # ), respectively. In addition, in Trulsen et al. (2020), the slope rate for the up-slope is 0.263 and the parameter is 0.05 for this region in the present study. The transition variation in this experiment is milder than that of Trulsen et al. (2020); correspondingly, the nonlinear increases slower. For the 140 realizations, all the freak waves are detected. In Figure 7B, it can be seen that the probability of freak wave occurrence is the highest when the mean kurtosis and the mean skewness achieve the maxima, as waves propagate close to the top crest of the breakwater, the maximum probability of occurrence for the freak wave can be up to about 1%. Figure 8 exhibits the maximum crest heights and wave heights for the freak waves. It is shown that the maximum crest heights and wave heights for cases with H s = 0.025 m change mildly and they fluctuate around the values of 0.04 m and 0.06 m, respectively as shown in Figure 8A. However, they vary significantly for the cases with the H s = 0.05 m, as shown in Figure 8B. The crest heights fluctuate between 0.04 m -0.08 m, and the wave heights fluctuate between 0.055 m -0.12 m. The comparison above illustrates that the B A

FIGURE 4
The evolutions of the mean significant wave heights for cases (A) with H s = 0.025 m and (B) with H s = 0.05 m. He et al. 10.3389/fmars.2023.1150896 Frontiers in Marine Science frontiersin.org significant wave height has an important effect on the local maximum crest heights and the maximum wave heights of the freak waves.
Following Figures 8, 9 shows the ratio, R, between the crest heights of the freak waves and the maximum wave heights. The results show that the values of R for cases with different significant wave heights fluctuate around a similar level, 0.7. It can indicate that the mean ratio between the crest height and the maximum wave height is close to a constant for the freak waves. It is noted that the definition of the ratio R is the same as the horizontal asymmetry factor shown in Bonmarin (1989) if the maximum wave height is the one containing the maximum crest height and is defined by zero-downcross analysis. In their study, the mean horizontal asymmetry factor ranges from 0.69 -0.77 for variable breaking types that occurred in deep water. It can be seen that the mean ratio R for the freak wave in finite water depth is within the scale of the mean horizontal asymmetry factor for breaking waves in deep water.

FIGURE 6
The evolutions of the mean asymmetry (A) and kurtosis (B).

FIGURE 7
The evolutions for the mean skewness (A) and the probability of occurrence for the freak waves (B).

FIGURE 5
The evolutions of the average value for the mean period, T m , for cases (A) with H s = 0.025 m and (B) with H s = 0.05 m. He et al. 10.3389/fmars.2023.1150896 Frontiers in Marine Science frontiersin.org Transformation of spectra over the smooth breakwater From the analysis above, it can be seen that the increasing nonlinearity makes the wave parameters change a lot (Figures 2, 6,  7), especially for the freak waves. In this section, the nonlinear interactions coupling their evolution of the wave over uneven bathymetry are studied in terms of spectra. The information on the frequency modulation can reflect nonlinearity well (Wu and Yao, 2004;Veltcheva and Soares, 2016). Taking the case RJ 7 as an example, the freak wave occurs at the location of x = 17 m (7 # ), and the evolutions of the Hilbert spectra and the wavelet-based bicoherences are presented in Figures 10, 11. In Figure 10A, the spectrum presents that the energy of waves mainly resides in the dominant frequency when the water depth is 0.45 m. Referring to Figure 11A, it also can be seen that the interaction between waves is weak. As the wave evolves along the up-slope, water depth decreases, and the frequency modulation actions appear apparently in Figures 10C, D, especially the latter. It can be found that the frequency modulation can increase up to about 4 Hz at 80 s around when the freak wave occurs.
Corresponding to the increasing frequency modulation, the interactions between waves grow significantly in Figures 11C, D. Comparison of Figures 10C, 11C, as the wave depth decreases, the interactions between primary waves, and that between the primary waves and the high-order harmonics start to grow; it can be reflected by b 2 (0.66, 0.66) = 0.27 and b 2 (1.32, 0.66) = 0.42. The variations in bicoherence illustrate that the nonlinear phase coupling between modes becomes obvious and the energy is transferred to 1.32 Hz and 1.98 Hz with the decreasing water depth. Accordingly, Figure 10C shows the same tendency at t = 10 s, 30 s, and 80 s around. When the water depth is reduced to 0.15 m, the Hilbert spectrum and waveletbased bicoherence are shown in Figures 10D, 11D, respectively. The active nonlinear interactions can be found, and the number of wave modes involved in quadratic phase coupling increases obviously. Figure 11D shows that the energy is transferred to the high-frequency contents significantly as the water depth decreases, for example,  Figure 10D, it can be seen that the frequency modulation acts strongly, especially, at t = 80 s around, 4 Hz component gains obvious energy, which can reflect the occurrence time for the triad wave interactions, corresponding to the components f 1 = 2.16 Hz, f 2 = 1.84 Hz, and f = 4 Hz.

FIGURE 8
The crest heights of the freak waves and the maximum wave heights for cases (A) with H s = 0.025 m and (B) with H s = 0.05 m, t freak is the occurrence time of the freak wave.

FIGURE 9
The ratio between the crest heights of the freak waves and the maximum wave heights for cases ( When the wave moves to the top peak of the bar, the water depth is 0.1 m and the wave breaks. Figure 10E shows that more energy is transferred to the high-frequency components; especially, there's even significant energy transferred to 5 Hz components. The substantial energy that appeared in high-frequency contents is mainly derived from the frequency modulation between wave components. Referring to Figure 11E, when the wave propagates in a depth of 0.1 m, the three-wave interactions appear more prominently, and many wave modes are involved in the nonlinear phase coupling, such as b 2 (0.66, 0.66) = 0.45, b 2 (1.32, 0.66) = 0.61, b 2 (2.64, 0.66) = 0.71, b 2 (3.3, 0.66) = 0.74, b 2 (1.98, 0.66) = 0.65, b 2 (1.98, 1.32) = 0.61, b 2 (2.64, 1.32) = 0.67, b 2 (1.14, 0.74) = 0.67, and b 2 (2.16, 1.84) = 0.75. By comparing the analysis results at 0.15 m water depth, it can be seen that the energy residing in components with 1.32, 1.98, 3.3, 3.96, and 4 Hz keeps increasing, which derives mainly from the interactions between the primary waves, the primary waves and their high-order harmonics, and the extent of the phase coupling is strong. However, there are also energy reductions in some components, such as 2.64, 3.3, and 1.88 Hz by tracing the decrease of b 2 (1.98, 0.66), b 2 (1.98, 1.32) and b 2 (1.14, 0.74), which illustrates that the quadratic phase coupling relationship between components with 1.98 and 0.66 Hz, 1.98 and 1.32 Hz and 1.14 and 0.74 Hz diminishes. It needs to be noted that the energy residing in the component with 3.3 Hz not only increases but decreases at this location, the quadratic phase coupling relationship between components with 2.64 and 0.66 Hz grows and transfers energy to 3.3 Hz; while the nonlinear quadratic phase coupling relationship between components with 1.98 and 1.32 Hz weakens. Referring to Figure 10E, the quadratic phase coupling mainly occurs before 100 s. At the whole top peak of the bar, Evolutions of the three-dimensional Hilbert spectra for the case RJ 7 at stations: (A) 2 # ; (B) 4 # ; (C) 6 # ; (D) 7 # ; (E) 8 # ; (F) 9 # ; (G) 10 # ; (H) 11 # ; (I) 14 # and (J) 17 # ; the freak wave occurs at station 7 # . He et al. 10.3389/fmars.2023.1150896 Figures 10E-G shows that the frequency modulation acts significantly, and the energy of the high-frequency components increases substantially. The significant variations also are exhibited in Figures 11E-G. In the de-shoaling region ( Figures 10H-J), an obvious reduction in frequency modulation is shown, the energy of all wave components is concentrated below 2 Hz, compared with Figures 10A-C, the energy distribution in this region is more uniform with time. The obvious changes are also exhibited in Figures 11H-J, the evolutions of the squared bicoherence spectra indicate that the range of the wave modes which are involved in the nonlinear phase coupling is reduced.
According to the above analysis, the frequency modulation is closely related to the nonlinear quadratic phase coupling as the water depth varies. When the water depth decreases, the quadratic nonlinear interactions strengthen, and the frequency modulation acts strongly. However, when the water depth increases, the quadratic nonlinear phase coupling and the frequency modulation decrease simultaneously.
For the same case, the difference between the energy spectra by FFT at locations x = 17 m (7 # ) and x = 10 m (2 # ) is described in Figure 12 to investigate the energy transformation in the process of the wave growing into a freak wave. There is an obvious increase in the low-and high-frequency contents, this also can be found in Figures 10, 11, this growth is mainly induced by the triad wave interactions as waves evolve over the shallower water depth.
After the extreme event, wave breaking occurs, and the difference between the spectra at the locations x = 18 m (8 # ) and x = 17 m (7 # ) is presented in Figure 13. It can be seen that the energy in the region of the dominant frequency is dissipated a lot after the wave breaking, and the high-frequency energy still increases. The energy variation illustrates that as the wave propagates over the depth of water getting shallower, the triad wave interactions increase, and the energy of the bound waves derived from the sum and the difference of the intrinsic frequencies grows. And when waves evolve coupling with the breaker, the energy of the dominant frequency band is dissipated, in the meanwhile, some bound waves are released to be free waves, these free waves may take part in the triad wave interactions as the water depth continues to decline. The increased energy in high-frequency contents can be referred to in Figures 10E, 11E.
When waves evolve on the top peak of the breakwater, the difference between the spectra at the locations x = 20.8 m (10 # ) and x = 18 m (8 # ) is presented in Figure 14. Energy in the regions of the dominant and low frequencies dissipates sharply, and that in the region of the high frequency still increases. The obvious decrease in energy of the dominant-frequency region may be induced by many small wave breaking. The increase in energy of the high-frequency contents can be found in Figures 10E, G, 11E-G). The triad wave interactions and the frequency modulation are all significant when waves evolve on the top of the breakwater, which indicates that the triad wave interactions depend on the frequency modulation. In the meanwhile, the high-frequency free waves maybe also a factor contributing to the increase of high-frequency energy.

Discussion and conclusion
In the present study, the nonlinear characteristics of the irregular random waves propagating over the uneven bathymetry are discussed. The effects of the varying water depths on the nonlinear interactions are investigated in terms of the wave surface evolutions, geometric parameters, and energy spectra. The shoal effect contributes to the sharp crests and the flat troughs. A value of about 0.7 of the ratio between the maximum crest heights and the maximum wave heights of the wave profiles for the freak waves is found, which illustrates that there is a local intense energy focus on the crest for the freak wave, therefore the extreme event will appear.
As waves evolve with the decreasing water depth, the nonlinear phase coupling increases by the bicoherence spectrum, and significant energy is transferred to high-frequency contents in terms of obvious frequency modulation. After wave breaking occurs, the energy of the dominant-frequency region is dissipated (ignoring the viscous and friction loss) and the energy in the highfrequency region is increased. By analyzing the variations of the nonlinear interactions as the wave evolves over the variable-depth bathymetry, the results show that the frequency modulation can reflect the triad wave interactions very well, and the abrupt increase in the frequency modulation can indicate the appearance of the extreme event.
In the present investigation, even though the nonlinear characteristics of geometric parameters, three-wave interactions, and frequency modulation are studied and some meaningful results are obtained here, some contents need to be learned further due to the limited conditions. More peak factors of the JONSWAP spectrum need to be considered for expanding the scope of the study. The essential relationship between the frequency modulation and the nonlinear phase coupling interaction in shallow water is still not revealed, therefore, more wave conditions are essential to study further to find more characteristics associated with them.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

FIGURE 12
The difference of spectra between the stations 7 # and 2 # for case RJ 7 .

FIGURE 13
The difference of spectra between the stations 8 # and 7 # for case RJ 7 .

FIGURE 14
The difference of spectra between the stations 10 # and 8 # for case RJ 7 . He et al. 10.3389/fmars.2023.1150896 Author contributions YH contributed to the data analysis and the paper writing. GW improved the grammar of the manuscript. HM improved the words. HC designed and conducted the experiment. JL reviewed the manuscript. GD provided the equipment and conditions for the experiment. All authors contributed to the article and approved the submitted version.

Funding
This work was financially supported by the Ocean Youth Talent