Dispersive characteristics of non-linear waves propagating and breaking over a mildly sloping laboratory beach

Abstract The dispersive characteristics of unidirectional irregular waves propagating and breaking over a mildly sloping beach are examined using a highly-resolved laboratory dataset. Cross-spectral analyses are used to determine the cross-shore evolution of (single-valued) dominant wavenumber κ and phase velocity c spectra, and lead to the identification of four different regimes of propagation: I - a linear regime where short waves mostly propagate as free components; II - a shoaling regime where non-linear effects at high harmonics are significant but primary components follow the linear wave dispersion relation; III - a shoaling regime near the mean breaking point location, where amplitude dispersion effects at primary components are important; IV - a surf zone regime, where all components propagate slightly faster than non-dispersive shallow water waves. Bispectral analyses performed onshore of the shoaling region show that the presence of forced energy at high harmonics, which originate from non-linear interactions between triads of frequencies, are responsible for the deviations of wavenumber and phase velocity spectral estimates from predictions by the linear dispersion relation, confirming the findings from previous field-based studies. A Boussinesq approximation of the non-linear energy exchanges between triads is then used to quantify the relative amount of forced energy at high harmonics and explain the differences in dispersion properties observed in the shoaling region between broad and narrow-band spectra. Larger relative amounts of forced energy at high frequencies, which suggest more efficient non-linear energy transfers, are found to be associated with larger deviations of dominant κ and c from predictions by the linear dispersion relation.

relation (Phillips, 1960;Longuet-Higgins and Stewart, 1962;Freilich et al., 1984). As opposed to non-linear 23 resonant interactions between quadruplets in deep water (Hasselmann, 1962), which require very large 24 distances to be effective, non-linear coupling between triads in nearshore areas are non-or near-resonant 25 and can be very efficient in transforming incident wave spectra over just few typical wavelengths (e.g., 26 see Freilich et al., 1984, and the references therein). As both forced (or 'bound') and free components 27 of directionally spread seas can co-exist in a wave field, there is no longer a unique relation between a 28 frequency and wavenumber (e.g. Herbers and Guza, 1994). When forced components dominate over a 29 region of the spectrum, large deviations from predictions by the linear wave dispersion relation can be 30 observed in (single-valued) dominant wavenumber and phase velocity spectra (Thornton and Guza, 1982;31 Freilich et al., 1984; Elgar and Guza, 1985b). In particular near the breaking point or in the surf zone, 32 most wave components of a typical sea-surface spectrum travel at the speed of non-dispersive shallow-water 33 waves (e.g., see Thornton and Guza, 1982; Elgar and Guza, 1985b; Catalán and Haller, 2008; Tissier et al., 34 2011), which is due to the dominance of amplitude dispersion effects over frequency ones (Herbers et al.,35 2002). As noted by Laing (1986), the deviations of measured wave phase speed from predictions by the 36 linear dispersion relation, discussed here for nearshore waves, are quite analogous to those observed in 37 growing seas (e.g., see Ramamonjiarisoa and Coantic, 1976;Mitsuyasu et al., 1979;Crawford et al., 1981; 38 Donelan et al., 1985). In such conditions, growing short-wave fields are dominated by modulated trains 39 of finite amplitude waves to which high-frequency components are bound (Lake and Yuen, 1978;Coantic 40 et al., 1981). 41 In practice, knowledge on the spatial structure of the wave field is generally lacking and the presence of 42 forced energy is therefore difficult to quantify. As forced components at high harmonics are characterised 43 by lower wavenumbers than free components of the same frequency, large errors from depth-inversion algo-  Martins et al., 2020b). This is related to the fact 49 that forced high harmonics are much less attenuated across the vertical than free components of the same and allowed the collection of current velocities at numerous locations across the wave flume (the data from 88 these current meters is not used here). This spatial resolution makes the GLOBEX dataset unique as it 89 allows to characterize and quantify non-linearities at various stages of the waves propagation (from a linear 90 situation up to the surf zone). The wave paddle steering signals included second-order wave generation and 91 the paddle was equipped with an Active Reflection Compensation system to absorb long waves radiated 92 from the beach. The water depth at the wave paddle was 0.85 m for all tests. 93 The present study uses free surface elevation measurements from the 70-min-long irregular wave tests 94 of the A series (A1, A2, A3, see Ruessink et al., 2013). During this series of tests, JONSWAP spectra 95 were imposed, covering moderate to energetic and broad to narrow-banded sea wave conditions, see table   96 1. The wave conditions along the flume for these runs are displayed in Fig. 2, through a range of second 97 and third-order wave parameters. Fig. 2a first shows the root-mean square wave height H rms computed 98 as (8 ζ 2 ) 1/2 where ζ is the free surface elevation and the overbar denotes the time-averaged operator. Note 99 that ζ was high-pass filtered using a cutoff frequency at 0.6 f p , with f p the peak wave frequency, so that 100 the bulk parameters shown in Fig. 2 are computed on the short-wave frequency band only. As conditions 101 were more energetic during A2, the mean breaking point was located farther seaward than during A1 and 102 A3. In the inner surf zone, waves were found to be depth-limited during all runs, which is also evidenced by (1) where H{·} is the Hilbert transform. S k and A s are measures of wave asymmetry along the horizontal and 107 vertical axes respectively and also inform on the energy content at high harmonics (Elgar, 1987 (Thornton and Guza, 1982;Freilich et al., 1984;Elgar and Guza, 1985b).

125
Let C x1, x2 denote the cross-spectrum computed with surface elevation timeseries from two gauges 126 located at positions x 1 and x 2 . The coherence coh(f) and phase φ(f ) spectra computed between x 1 and x 2 127 are then respectively given by where Re{·} and Im{·} are the real and imaginary parts of the cross-spectra respectively and * denotes the 129 complex conjugate. As shown in Fig. 3a with an example of cross-spectral analysis performed during A2 in 130 the deepest region of the flume (x ∼ 10 m), the phase spectra φ(f ) provides a phase lag per frequency that 131 is bounded between -π/2 and π/2. The time delay (in sec) per frequency is obtained from the unwrapped 132 phase φ unw which, in the case of progressive waves propagating in one dimension, is easily retrieved from 133 the phase jumps (e.g., see Fig. 3a). The wavenumber κ(f ) and (cross-shore) phase velocity c(f ) spectra 134 are then readily computed as where ∆x is the spacing between the two wave gauges. As in Herbers et al. (2002), κ refers to the 136 single-valued wavenumber modulus, representative of a mixed sea-state composed of both free and forced 137 components (i.e., for a given angular frequency ω = 2πf , energy in the (ω, k) space is spread across several 138 wavenumbers k). In practice, κ and c provide estimates at x = (x 1 + x 2 )/2 of the dominant wavenumber (in 139 an energy-averaged sense) and the corresponding propagation speed respectively, as shown by the analysis 140 on synthetic data in Appendix A. This aspect is also discussed in Section 5. Cross-spectra were computed 141 using Welch's method and 63 Hann-windowed records of 128 seconds, which were overlapping by 50%. This 142 resulted in each spectral estimate having approximately 70 equivalent degrees of freedom and a spectral 143 resolution of 0.008 Hz.

144
As observed by many authors in the past, the coherence spectra at high harmonics were found to be 145 quite sensitive to the spacing between the two adjacent wave gauges. This is especially true in the deepest 146 parts of the wave flume, where energy levels at these frequencies are quite low and the spacing between 147 the two gauges can represent several wavelengths. An example is provided in Fig. 3b, which shows the 148 coherence spectra computed over the flat section of the flume during A2 with ∆x = 0.93L p , L p being the 149 peak wavelength given by linear wave theory. The coherence remains very high (coh > 0.95) at frequencies 150 between 0.6 and 1.5 f p , which explain over 86% of the short-wave variance at this location. However, valleys 151 in the coherence can be observed around 2.5 f p for this particular spacing configuration and the coherence 152 weakens quickly after 2.8 f p (Fig. 3b). Since using a different spacing slightly shifts the coherence 'valleys', 153 several combinations of wave gauges and spacing were used to obtain spectral estimates at a single location 154 with coherence higher than typically 0.5 at all frequencies (similar idea as that used by Herbers et al.,155 2002, with their field observations). After the removal of data associated with a coherence inferior to 0.5 156 (coh 2 0.25 in Fig. 3b), spectral estimates were ensemble-averaged. When non-linearities are strong, 157 such as near the breaking point or in the surf zone, the coherence remains high up to 4-5 f p , as long as 158 ∆x is taken sufficiently small (typically 0.2-0.4 L p in those regions). In such cases the ensemble-averaging 159 procedure is not really necessary but it was used all along the flume for consistency in the analysis. The power bispectrum of discretely sampled data corresponds to a representation in the frequency 162 domain of its third-order cumulant or moment. As it provides information on the strength of the phase 163 coupling between triads of frequencies f 1 , f 2 and f 1 + f 2 , the bispectrum is a useful and powerful tool for 164 studying non-linear phenomena in ocean waves (Hasselmann et al., 1963;Elgar and Guza, 1985a). Here, 165 we use the definition given by Kim and Powers (1979): where A(f ) are the complex Fourier coefficients and E is an expected, or ensemble-average, value. Similar 167 to the cross-spectrum phase, the biphase is obtained from the bispectrum as It is particularly insightful to recast the biphase as a function of the Fourier coefficients phases θ(f ): Kim et al., 1980;Elgar and Guza, 1985a). The amount of 170 energy transfer between near-resonant components depends on their relative phases, whose information is 171 contained in the biphase (Hasselmann et al., 1963;Kim et al., 1980).

172
Bispectra were computed on the free surface elevation signals down-sampled to 16 Hz by averaging 173 estimates from 126 Hann-windowed records of 128 seconds, which were overlapping by 75%. Statistical  the present modelling approach to A2 and A3 since conditions during A1 are too dispersive.

189
For unidirectional waves propagating shoreward on an alongshore-uniform beach and assuming a weak 190 reflection at the shoreline, a balance between the cross-shore gradient of the energy flux spectrum F (f ), 191 a source term S nl (f ) quantifying the non-linear energy exchanges between triads, and a dissipation term 192 S dis (f ) reads (e.g., Eq. 1 of Herbers et al., 2000) 193 Breaking processes are ignored so that the dissipation term S dis reduces to the energy losses by bottom 194 friction S f r , here simply modelled after Thornton and Guza (1983): where h is the mean water depth and |k| is the wavenumber modulus obtained from the linear wave  Assuming that the wave field is weakly non-linear, weakly dispersive, and that these effects are of similar : Herbers and Burton (1997, their Eq. 22a), the energy in Eq. 9 is assumed to propagate at the shallow 213 water wave speed so that F (f ) = ρgE(f ) √ gh, with ρ the water density and g the acceleration of gravity.

214
After these considerations, Eq. 9 simplifies to Integrating this equation (in space) between the location of the first gauge x 0 and any location x 1 prior 216 to the mean breaking point location yields the following expression for the energy density spectrum at x 1 217 E(x 1 , f ): Finally, the ratio x1 x0 S nl dx/F (x 1 ), which represents the energy flux received (or lost) via non-linear 219 coupling between x 0 and x 1 over the total energy flux at x 1 , is used in the following as an approximation 220 of the relative amount of forced energy at x 1 . Outside the surf zone, where dissipative processes dominate, 221 this estimation was found to be more reliable than the bicoherence, which is often used as a proxy for such 222 an estimation but lacks general consensus upon its definition, see Appendix C for more details. In this section, we present and describe the main results from the cross-spectral and bispectral analyses. links the spatial and temporal information of a linear wave field: In the following, the subscript 'L' is used throughout the manuscript to refer to κ or c values that are solution 234 to Eq. 14 (i.e. κ L and c L = ω/κ L ). For conciseness, the propagation as linear waves is discussed using 235 the results from the most non-linear case (A2), already shown in Fig. 3. The other stages of propagation 236 focus on the differences between broad and narrow spectra, i.e. between A2 and A3. Note that κ and c 237 characteristics during A1 are very similar to those obtained during A2 for similar U r numbers. Finally,

238
given the relatively high spatial resolution, frequency-wavenumber power spectra P (ω, k) (e.g., see Redor   for A3, which explains the significant growth of the third harmonic 3f p observed for that test (Fig. 5b). show a strong coupling with f p , which is the consequence of a broader spectrum for A2 compared to A3. see Elgar and Guza, 1985a).

294
In the infragravity band, wave components are still bound to short-wave groups as indicated by the 295 good match between spectral estimates of phase velocity and the propagation speed of short-wave groups.

296
Energy transfers towards the infragravity band concentrate at frequencies around 0.1-0.2f p (Fig. 5a and 297 5b) and principally originate from strong difference interactions which transfer energy from f p towards that 298 infragravity frequency and components at frequency slightly lower than f p (see the negative imaginary part 299 of B * along the f p anti-diagonal, Fig. 5g and 5h).  In the surf zone (U r taken at 3.2), wavenumber and phase velocity spectra are frequency-independent 329 (Fig. 7e-f). The large differences observed between the measured phase velocity at peak frequency and that , we analyse non-linear energy transfers to 2f p and 3f p in more details here.

362
Wave amplitudes at f p , 2f p and 3f p (a fp , a 2fp and a 3fp respectively) computed from energy spectra 363 modelled with Eq. 13 are compared against observations in Fig. 8. Overall, the Boussinesq approach 364 of Herbers and Burton (1997) for the non-linear energy transfers between triads accurately predicts the 365 growth of second and third harmonics across the shoaling zone for both broad and narrow-band wave tests.

366
The cross-shore evolution of a fp , a 2fp and a 3fp are very well described up to the mean breaking point, with 367 mean absolute percentage errors (MAPE) lower than 5% and 13% for a 2fp and a 3fp respectively. With 368 respect to the amount of (free) energy imposed at the paddle, the narrow-banded conditions during A3 369 promote more efficient energy transfers, with growths of a 2fp and a 3fp by a factor 3 and 10 respectively, 370 while these factors are only between 1.5 and 2 during A2.

371
A closer look at the cross-shore evolution of the source terms (Fig. 9) shows that friction effects are 372 negligible at high frequencies but not around the peak frequency where the energy dissipated through 373 bottom friction is of similar order than the energy lost via non-linear coupling. During A3 (Fig. 9b), the 374 steady increase of S nl (2f p ) indicates a gradual growth of a 2fp in the shoaling region (see also Fig. 8b). This

389
The effect of varying relative amounts of forced energy not only explain the differences observed in κ and 390 c spectra at high frequencies between A2 and A3 but also their variation across frequencies. During A3 for 391 instance, non-linear energy transfers in the short-wave frequency band were predominantly towards 2f p and 392 3f p (Fig. 5h and 6h). At these harmonics, energy is predominantly forced (> 70% around 2f p and > 90% The situation where both free and forced components exist in a wave field is not entirely clear as far 480 as the cross-spectral analysis is concerned. Some studies consider the cross-spectral estimates to be biased 481 towards bound harmonics (e.g., see Lake and Yuen, 1978; Thornton and Guza, 1982), arguing that the 482 travelling distance between two gauges is reduced for bound components compared to free ones. This is 483 evidently exacerbated when using sub-surface pressure sensors or current meters since bound harmonics 484 are much less attenuated and will dominate the spectrum at some depth.

485
To analyse the sensitivity of the cross-spectral analysis to different relative levels of forced energy, some 486 tests are performed here on synthetic data. Starting at a position x 1 , a surface elevation timeseries following 487 a JONSWAP spectrum with similar parameters to the A3 test is generated (Fig. A1). A Gaussian-shaped 488 perturbation around 2 f p is added to represent forced harmonics: for any f around 2 f p , the forced harmonic The results from the cross-spectral analysis performed on the synthetic timeseries at x 1 and x 2 are pro-495 vided in Fig. A2 as the wavenumber and phase speed at 2 f p shown as a function of a 2fp, f orced /a 2fp, f ree .

496
In the absence of forced components (a 2fp, f orced /a 2fp, f ree = 0), components around 2 f p follow the dis- The wave skewness and asymmetry are third-order moments characterizing the wave shape. These 506 parameters are directly related to the energy content at high frequencies of the spectra and bispectra 507 (Elgar, 1987). The fact that the statistical (see Eq. 1 and 2) and the bispectrum-based definitions are 508 theoretically equivalent provide a means to validate the bispectrum's calculations.

509
Due to its symmetry properties, the bispectrum can be uniquely defined in a single octant in the with n > l and n + l < N .

516
The cross-comparison between statistical and bispectrum-based definitions of the wave skewness and 517 asymmetry computed for A3 is provided in Fig. B1 and shows a perfect match between both definitions, 518 thus validating the present computations of the bispectrum. Note that for these comparisons, blocks were 519 not tapered and rectangular windows were used for computing the Fast Fourier Transforms. Using any 520 other types of windows resulted in small differences between the two definitions of both third-order wave 521 parameters.

522
Appendix C. Normalised bispectra (bicoherence) 523 We generally seek a normalisation of the bispectrum so that it takes 1 as value when there is a full 524 coupling between the three components involved and 0 when there is none. In the case of surface gravity 525 waves, this would correspond to a situation with the only presence of forced or free energy respectively. studying non-linearities in the nearshore area but was shown to lead to values greater than 2 (e.g. Elgar

531
and Guza, 1985a). In the present study, this was also the case, but only when merging across frequencies 532 was performed beforehand. Nonetheless, a slightly different normalisation that was proposed in Hagihira    Table 1: Significant wave height (Hs) and peak wave frequency (fp) for the three tests considered in this study. The peak enhancement γ of the JONSWAP spectra characterizes its spectral bandwidth. A value of 3.3 corresponds to broad spectra while a narrow-band spectra is imposed with γ = 20. : Second and third-order short-wave parameters during tests A1, A2 and A3: a) root-mean square wave height Hrms computed as (8 ζ 2 ) 1/2 ; b) corresponding wave amplitude to water depth ratio = Hrms/ √ 2h; c) wave skewness S k ; d) wave asymmetry As and e) Ursell number Ur computed as /µ. Note that these are short-wave parameters, computed using the high-pass filtered surface elevation signal (frequency cutoff at 0.6 fp). The spacing between the adjacent gauges is ∆x = 0.93Lp, with Lp the peak wavelength given by the linear wave dispersion relation (Eq. 14). The phase velocity spectrum is normalised by c L (fp), the phase velocity predicted by linear wave theory at the peak frequency. The red dashed lines in panels a, c and d correspond to the values given by the linear wave dispersion relation. The separation between infragravity and short wave frequencies (0.6 fp) is shown as dashed black line. In panel b, the gray dashed line corresponds to the coherence-squared threshold used for computing ensemble average spectral estimated of κ and c. In panel d, the orange line corresponds to the short-wave envelop propagation speed. Figure 4: Frequency-wavenumber surface elevation power spectra P (ω, k) computed for A2 (left) and A3 (right) at the four stages of propagation described in the paper (I: panels a-b; II: panels c-d; III: panels e-f and IV: panels g-h). For A2, these power spectra were computed at I: x = 9.8 ± 2.8 m; II: x = 38.9 ± 3.6 m; III: x = 54.5 ± 3.5 m and IV: x = 76 ± 2.4 m. For A3, the cross-shore locations were I: x = 9.8 ± 2.8 m; II: x = 53.4 ± 2.4 m; III: x = 62.3 ± 3.5 m and IV: x = 77 ± 2.4 m. Each black contour line correspond to a power of 10. The separation between infragravity and short-wave frequencies (0.6 fp) is shown as the dashed black line. In all panels, the red curves correspond to the linear wave dispersion relation. In panels a-d, the dashed grey curves correspond to the dispersion relation for the second harmonic, and was computed based on the assumption that it travels at a similar speed as its principal component. In panels e-f, the grey curve corresponds to the shallow water wave dispersion relation (c = √ gh) while the dashed green line in panels g-h refer to the modified one (c = gh(1 + )). : From top to bottom: surface elevation energy density spectra, dimensionless wavenumber spectra, phase velocity spectra normalised by c L (fp) and bispectra computed at Ur ∼ 0.3 (stage II) for wave tests A2 (left panels) and A3 (right panels). The cross-shore location corresponding to this stage is indicated in panels a and b (cf figure 2e). The wavenumber and phase velocity spectra were computed with five spacing configurations: each point corresponds to the ensemble-averaged value and the error bars correspond to the standard deviation. For readability, only one on three data points are shown. Red lines in panels c-f correspond to values given by the linear wave dispersion relation. The separation between infragravity and short-wave frequencies (0.6 fp) is shown as the dashed black line. In panels c-d, the red dashed line corresponds to κ(f ) = 2πf /c(fp).
In panels e-f, the orange line corresponds to the short-wave envelop propagation speed denoted cenv.    Figure A1: Energy density spectra imposed at a fictive x 1 position, with different ratio of forced to free second harmonic amplitude. In the absence of forced components, the forcing corresponds to a JONSWAP spectrum with parameter similar to A3 (see Table 1).