Kévin Martins
email: kevin.martins@u-bordeaux.fr

Philippe Bonneton

Hervé Michallet

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

come

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

Introduction

Wind-generated surface gravity waves (hereafter short waves) are the principal driver of nearshore dynamics. Close to shore, short waves eventually break and through this process, they enhance the vertical and horizontal mixing of the water column (e.g., [START_REF] Ting | Dynamics of surf-zone turbulence in a spilling breaker[END_REF][START_REF] Drazen | Turbulence and mixing in unsteady breaking surface waves[END_REF][START_REF] Clark | Vorticity generation by short-crested wave breaking[END_REF], drive a setup near the shoreline (e.g., [START_REF] Longuet-Higgins | Radiation stresses in water waves; a physical discussion, with applications[END_REF][START_REF] Stive | A study of radiation stress and set-up in the nearshore region[END_REF] and control the nearshore circulation at various temporal scales (e.g., [START_REF] Svendsen | Mass flux and undertow in a surf zone[END_REF][START_REF] Peregrine | Vorticity and surf zone currents[END_REF][START_REF] Bühler | Wave-driven currents and vortex dynamics on barred beaches[END_REF][START_REF] Bonneton | Large-scale vorticity generation due to dissipating waves in the surf zone[END_REF][START_REF] Castelle | Rip current types, circulation and hazard[END_REF]. At first order, linear wave theory correctly predicts a number of physical processes associated with the propagation of short waves, such as the refraction or shoaling of directionally spread seas (e.g. [START_REF] Longuet-Higgins | The refraction of sea waves in shallow water[END_REF][START_REF] Guza | Local and shoaled comparisons of sea surface elevations, pressures, and velocities[END_REF]Elgar et al., 1990, and many others). However, as waves shoal and interact with a sloping, mobile bed, non-linear processes become dominant. These are responsible for the changes observed in wave shape, from nearly symmetric to more (horizontally) skewed in the shoaling region and more (vertically) asymmetric prior to breaking and in the surf zone (e.g. Elgar and Guza, 1985a,b;[START_REF] Doering | Parametrization of orbital velocity asymmetries of shoaling and breaking waves using bispectral analysis[END_REF][START_REF] Michallet | Kinematics of waves propagating and breaking on a barred beach[END_REF][START_REF] Rocha | Improving the parameterization of wave nonlinearitiesthe importance of wave steepness, spectral bandwidth and beach slope[END_REF]. Non-linear effects not only affect the sea surface elevation but also the near-bottom wave orbital velocities and thus play a crucial role in short-and long-term beach morphodynamics (Doering andBowen, 1986, 1995;[START_REF] Hoefel | Wave-induced sediment transport and sandbar migration[END_REF][START_REF] Berni | Surf zone cross-shore boundary layer velocity asymmetry and skewness: An experimental study on a mobile bed[END_REF][START_REF] Van Der Zanden | Bedload and suspended load contributions to breaker bar morphodynamics[END_REF].

Weakly non-linear triad interactions (e.g. [START_REF] Phillips | On the dynamics of unsteady gravity waves of finite amplitude. Part 1. The elementary interactions[END_REF][START_REF] Freilich | Nonlinear effects on shoaling surface gravity waves[END_REF]Elgar and Guza, 1985a) occurring as short waves propagate landward over a sloping bottom are responsible for these changes in the wave field. The interaction of two primary components of frequencies f 1 and f 2 excite a secondary component f (either sum f 1 + f 2 or difference f 1 -f 2 ), which is bound to the statistically independent primary components. As such, the bound wave component f does not follow the linear wave dispersion relation [START_REF] Phillips | On the dynamics of unsteady gravity waves of finite amplitude. Part 1. The elementary interactions[END_REF][START_REF] Longuet-Higgins | Radiation stress and mass transport in gravity waves, with application to 'surf beats[END_REF][START_REF] Freilich | Nonlinear effects on shoaling surface gravity waves[END_REF]. As opposed to non-linear resonant interactions between quadruplets in deep water [START_REF] Hasselmann | On the non-linear energy transfer in a gravity-wave spectrum. Part 1. General theory[END_REF], which require very large distances to be effective, non-linear coupling between triads in nearshore areas are non-or near-resonant and can be very efficient in transforming incident wave spectra over just few typical wavelengths (e.g., see [START_REF] Freilich | Nonlinear effects on shoaling surface gravity waves[END_REF], and the references therein). As both forced (or 'bound') and free components of directionally spread seas can co-exist in a wave field, there is no longer a unique relation between a frequency and wavenumber (e.g. Herbers and Guza, 1994). When forced components dominate over a region of the spectrum, large deviations from predictions by the linear wave dispersion relation can be observed in (single-valued) dominant wavenumber and phase velocity spectra [START_REF] Thornton | Energy saturation and phase speeds measured on a natural beach[END_REF][START_REF] Freilich | Nonlinear effects on shoaling surface gravity waves[END_REF]Elgar and Guza, 1985b). In particular near the breaking point or in the surf zone, most wave components of a typical sea-surface spectrum travel at the speed of non-dispersive shallow-water waves (e.g., see [START_REF] Thornton | Energy saturation and phase speeds measured on a natural beach[END_REF]Elgar and Guza, 1985b;[START_REF] Catalán | Remote sensing of breaking wave phase speeds with application to non-linear depth inversions[END_REF][START_REF] Tissier | Field measurements and non-linear prediction of wave celerity in the surf zone[END_REF], which is due to the dominance of amplitude dispersion effects over frequency ones [START_REF] Herbers | Nonlinear dispersion of surface gravity waves in shallow water[END_REF]. As noted by [START_REF] Laing | Nonlinear properties of random gravity waves in water of finite depth[END_REF], the deviations of measured wave phase speed from predictions by the linear dispersion relation, discussed here for nearshore waves, are quite analogous to those observed in growing seas (e.g., see [START_REF] Ramamonjiarisoa | Loi experimentale de dispersion des vagues produites par le vent sur une faible longueur d'action[END_REF][START_REF] Mitsuyasu | On the dispersion relation of random gravity waves. Part 2. An experiment[END_REF][START_REF] Crawford | Effects of nonlinearity and spectral bandwidth on the dispersion relation and component phase speeds of surface gravity waves[END_REF][START_REF] Donelan | Directional spectra of wind-generated ocean waves[END_REF]. In such conditions, growing short-wave fields are dominated by modulated trains of finite amplitude waves to which high-frequency components are bound [START_REF] Lake | A new model for nonlinear wind waves. Part 1. Physical model and experimental evidence[END_REF][START_REF] Coantic | Wind-water tunnel simulation of small-scale ocean-atmosphere interactions[END_REF].

In practice, knowledge on the spatial structure of the wave field is generally lacking and the presence of forced energy is therefore difficult to quantify. As forced components at high harmonics are characterised by lower wavenumbers than free components of the same frequency, large errors from depth-inversion algorithms based on the linear wave dispersion relation can be expected in regions where non-linear effects are important (e.g., see Holland, 2001;[START_REF] Brodie | Evaluation of video-based linear depth inversion performance and applications using altimeters and hydrographic surveys in a wide range of environmental conditions[END_REF]. The over-predictions of the dominant wavenumbers at high frequencies also explain the commonly reported 'blow-up' when reconstructing the free surface elevation from sub-surface pressure measurements with the linear dispersion relation [START_REF] Bonneton | Recovering water wave elevation from pressure measurements[END_REF][START_REF] Bonneton | A nonlinear weakly dispersive method for recovering the elevation of irrotational surface waves from pressure measurements[END_REF][START_REF] Mouragues | Field data-based evaluation of methods for recovering surface wave elevation from pressure measurements[END_REF]Martins et al., 2020b). This is related to the fact that forced high harmonics are much less attenuated across the vertical than free components of the same frequency (e.g. [START_REF] Herbers | Wind-wave nonlinearity observed at the sea floor. Part I: Forced-wave energy[END_REF]Herbers et al., 1992). Nonetheless, most field-based studies on non-linear wave transformation in the shoaling region employed sub-surface hydrodynamic data (whether pressure or orbital wave velocity), not corrected or corrected for depth-attenuation using the linear wave dispersion relation. Field-based studies on wave non-linearity also suffer from other limitations such as a poor spatial resolution and the distance over which waves can be studied. In particular, the cross-shore location where non-linear effects at high harmonics become predominant remains largely unknown.

The present paper uses a high-resolution laboratory dataset (GLOBEX, see [START_REF] Ruessink | GLOBEX: wave dynamics on a gently sloping laboratory beach[END_REF] to study the dispersive properties of irregular waves propagating and breaking over a mildly sloping beach.

Besides confirming past findings, the GLOBEX dataset stands out from previously-published field observations for several reasons. The free surface is directly measured with wave gauges and it is highly-resolved in space (several points per wavelength at any stage of propagation). The former aspect removes uncertainties currently existing on energy levels at high harmonics as measured in the field by sub-surface pressure sensors, and where the choice of the surface elevation reconstruction method has a strong influence [START_REF] Bonneton | A nonlinear weakly dispersive method for recovering the elevation of irrotational surface waves from pressure measurements[END_REF][START_REF] Mouragues | Field data-based evaluation of methods for recovering surface wave elevation from pressure measurements[END_REF]Martins et al., 2020b). Furthermore, the experiments considered unidirectional irregular waves, which removes uncertainty about directional effects. In section 2, the high spatial and temporal resolution experimental dataset collected during GLOBEX is briefly presented. Section 3 introduces the cross-spectral and bispectral analysis techniques and describes the weakly non-linear numerical approach employed here for predicting the cross-shore evolution of energy spectra. In Section 4, the cross-spectral analysis is performed on the surface elevation data from adjacent wave gauges to extract dominant wavenumber spectra κ(f ), phase velocity spectra c(f ) and their evolution across the entire wave flume. From this analysis, we identify four regimes of propagation ranging from a linear up to a surf zone situation, where wavenumber and phase velocity spectra display specific characteristics. In section 5, the bispectral analysis is used to quantify non-linear energy transfers towards harmonics, which play a fundamental role in the patterns observed at high frequencies in wavenumber and phase velocity spectra. The dominant wavenumber is shown to vary with the amount of forced energy at a particular frequency, with larger deviations from the linear wave dispersion expected for higher forced-to-free energy ratios. Finally, section 6 briefly discusses the results and provides the concluding remarks of this study.

Experimental dataset

The Gently sLOping Beach Experiment (GLOBEX) project was performed in a 110-m-long, 1-m-wide, and 1.2-m-high wave flume, located in the Scheldegoot in Delft, the Netherlands [START_REF] Ruessink | GLOBEX: wave dynamics on a gently sloping laboratory beach[END_REF].

The experiments aimed at collecting high-resolution data of free surface elevation and current velocities in order to study infragravity wave dynamics and short-wave propagation and non-linearities (e.g., see [START_REF] De Bakker | Nonlinear infragravity-wave interactions on a gently sloping laboratory beach[END_REF][START_REF] Tissier | Infragravity-wave modulation of short-wave celerity in the surf zone[END_REF][START_REF] Rocha | Improving the parameterization of wave nonlinearitiesthe importance of wave steepness, spectral bandwidth and beach slope[END_REF]. A combination of 21 capacitance and resistance-type of wave gauges sampling at 128 Hz were deployed along the low-sloping 1:80 concrete beach to measure the free surface elevation (see Fig. 1). 18 of these wave gauges were mounted on movable trolleys, which were repositioned after the repetition of each wave test in order to reach the desired spatial resolution (10 repetitions per test). Similarly, five electromagnetic current meters were fixed to the trolleys and allowed the collection of current velocities at numerous locations across the wave flume (the data from these current meters is not used here). This spatial resolution makes the GLOBEX dataset unique as it allows to characterize and quantify non-linearities at various stages of the waves propagation (from a linear situation up to the surf zone). The wave paddle steering signals included second-order wave generation and the paddle was equipped with an Active Reflection Compensation system to absorb long waves radiated from the beach. The water depth at the wave paddle was 0.85 m for all tests.

The present study uses free surface elevation measurements from the 70-min-long irregular wave tests of the A series (A1, A2, A3, see [START_REF] Ruessink | GLOBEX: wave dynamics on a gently sloping laboratory beach[END_REF]. During this series of tests, JONSWAP spectra were imposed, covering moderate to energetic and broad to narrow-banded sea wave conditions, see table 1. The wave conditions along the flume for these runs are displayed in Fig. 2, through a range of second and third-order wave parameters. Fig. 2a first shows the root-mean square wave height H rms computed as (8 ζ 2 ) 1/2 where ζ is the free surface elevation and the overbar denotes the time-averaged operator. Note that ζ was high-pass filtered using a cutoff frequency at 0.6 f p , with f p the peak wave frequency, so that the bulk parameters shown in Fig. 2 are computed on the short-wave frequency band only. As conditions were more energetic during A2, the mean breaking point was located farther seaward than during A1 and A3. In the inner surf zone, waves were found to be depth-limited during all runs, which is also evidenced by the near-constant values reached by the non-linearity parameter = H rms / √ 2h in this region of the wave flume ( ∼ 0.3, see Fig. 2b), where h is the mean water depth. Fig. 2c-d show the short-wave skewness S k and asymmetry A s computed respectively as

S k = (ζ -ζ) 3 (ζ -ζ) 2 3/2 .
(1)

A s = - (H(ζ -ζ)) 3 (ζ -ζ) 2 3/2 , (2)
where H{•} is the Hilbert transform. S k and A s are measures of wave asymmetry along the horizontal and vertical axes respectively and also inform on the energy content at high harmonics [START_REF] Elgar | Relationships involving third moments and bispectra of a harmonic process[END_REF]Martins et al., 2020b;[START_REF] De Wit | The relationship between sea-swell bound wave height and wave shape[END_REF]. For the same imposed wave height, A1 and A3 differed in the shoaling region, with a more intense shoaling process noted during A3 (Fig. 2a and2c), as expected for longer waves. This is also explained by the narrower spectrum imposed during A3, which favoured non-linear energy transfers to higher harmonics as compared to A1 [START_REF] De Bakker | Nonlinear infragravity-wave interactions on a gently sloping laboratory beach[END_REF]. Short waves remain (vertically) symmetric (A s ∼ 0) until the location where the largest waves start breaking. Short-wave asymmetry is maximal in the inner surf zone, where short waves pitch forward and display the commonly observed sawtooth-like shape. Finally, the Ursell number U r is shown in Fig. 2e. Here, U r is defined as the ratio between the non-linearity and dispersion µ = (κ p h) 2 parameters, where κ p is the peak wavenumber given by the linear wave dispersion relation.

Methods

Computation of wavenumber and phase velocity spectra

Cross-spectral analysis between adjacent wave gauges is used to compute the dominant wavenumber and phase velocity spectra across the wave flume. As this approach provides phase differences (or delay) between two signals in the frequency domain (e.g., see [START_REF] Ochi | Spectral analysis[END_REF], it has been successfully used in the past to study, in both the laboratory and field, the dispersive properties of ocean waves propagating in deep (e.g., see [START_REF] Ramamonjiarisoa | Loi experimentale de dispersion des vagues produites par le vent sur une faible longueur d'action[END_REF][START_REF] Mitsuyasu | On the dispersion relation of random gravity waves. Part 2. An experiment[END_REF] and intermediate to shallow waters [START_REF] Thornton | Energy saturation and phase speeds measured on a natural beach[END_REF][START_REF] Freilich | Nonlinear effects on shoaling surface gravity waves[END_REF]Elgar and Guza, 1985b).

Let C x1, x2 denote the cross-spectrum computed with surface elevation timeseries from two gauges located at positions x 1 and x 2 . The coherence coh(f) and phase φ(f ) spectra computed between x 1 and x 2 are then respectively given by

coh x1, x2 (f) = C x1, x2 (f )C * x1, x2 (f ) C x1, x1 (f )C x2, x2 (f ) 1/2 (3) φ x1, x2 (f ) = arctan Im{C x1, x2 (f )} Re{C x1, x2 (f )} , (4)
where Re{•} and Im{•} are the real and imaginary parts of the cross-spectra respectively and * denotes the complex conjugate. As shown in Fig. 3a with an example of cross-spectral analysis performed during A2 in the deepest region of the flume (x ∼ 10 m), the phase spectra φ(f ) provides a phase lag per frequency that is bounded between -π/2 and π/2. The time delay (in sec) per frequency is obtained from the unwrapped phase φ unw which, in the case of progressive waves propagating in one dimension, is easily retrieved from the phase jumps (e.g., see Fig. 3a). The wavenumber κ(f ) and (cross-shore) phase velocity c(f ) spectra are then readily computed as

κ(f ) = φ unw x1, x2 (f ) ∆x (5) c(f ) = 2πf ∆x φ unw x1, x2 (f ), ( 6
)
where ∆x is the spacing between the two wave gauges. As in [START_REF] Herbers | Nonlinear dispersion of surface gravity waves in shallow water[END_REF], κ refers to the single-valued wavenumber modulus, representative of a mixed sea-state composed of both free and forced components (i.e., for a given angular frequency ω = 2πf , energy in the (ω, k) space is spread across several wavenumbers k). In practice, κ and c provide estimates at x = (x 1 + x 2 )/2 of the dominant wavenumber (in an energy-averaged sense) and the corresponding propagation speed respectively, as shown by the analysis on synthetic data in Appendix A. This aspect is also discussed in Section 5. Cross-spectra were computed using Welch's method and 63 Hann-windowed records of 128 seconds, which were overlapping by 50%. This resulted in each spectral estimate having approximately 70 equivalent degrees of freedom and a spectral resolution of 0.008 Hz.

As observed by many authors in the past, the coherence spectra at high harmonics were found to be quite sensitive to the spacing between the two adjacent wave gauges. This is especially true in the deepest parts of the wave flume, where energy levels at these frequencies are quite low and the spacing between the two gauges can represent several wavelengths. An example is provided in Fig. 3b, which shows the coherence spectra computed over the flat section of the flume during A2 with ∆x = 0.93L p , L p being the peak wavelength given by linear wave theory. The coherence remains very high (coh > 0.95) at frequencies between 0.6 and 1.5 f p , which explain over 86% of the short-wave variance at this location. However, valleys in the coherence can be observed around 2.5 f p for this particular spacing configuration and the coherence weakens quickly after 2.8 f p (Fig. 3b). Since using a different spacing slightly shifts the coherence 'valleys', several combinations of wave gauges and spacing were used to obtain spectral estimates at a single location with coherence higher than typically 0.5 at all frequencies (similar idea as that used by [START_REF] Herbers | Nonlinear dispersion of surface gravity waves in shallow water[END_REF], with their field observations). After the removal of data associated with a coherence inferior to 0.5 (coh 2 0.25 in Fig. 3b), spectral estimates were ensemble-averaged. When non-linearities are strong, such as near the breaking point or in the surf zone, the coherence remains high up to 4-5 f p , as long as ∆x is taken sufficiently small (typically 0.2-0.4 L p in those regions). In such cases the ensemble-averaging procedure is not really necessary but it was used all along the flume for consistency in the analysis.

Computation of bispectra

The power bispectrum of discretely sampled data corresponds to a representation in the frequency domain of its third-order cumulant or moment. As it provides information on the strength of the phase coupling between triads of frequencies f 1 , f 2 and f 1 + f 2 , the bispectrum is a useful and powerful tool for studying non-linear phenomena in ocean waves [START_REF] Hasselmann | Bispectra of ocean waves[END_REF]Elgar and Guza, 1985a). Here, we use the definition given by [START_REF] Kim | Digital Bispectral Analysis and Its Applications to Nonlinear Wave Interactions[END_REF]:

B( f 1 , f 2 ) = E [A(f 1 ) A(f 2 ) A * (f 1 + f 2 )] , (7)
where A(f ) are the complex Fourier coefficients and E is an expected, or ensemble-average, value. Similar to the cross-spectrum phase, the biphase is obtained from the bispectrum as

β( f 1 , f 2 ) = arctan Im{B( f 1 , f 2 )} Re{B( f 1 , f 2 )} . (8)
It is particularly insightful to recast the biphase as a function of the Fourier coefficients phases θ(f ): [START_REF] Kim | Bispectrum and nonlinear wave coupling[END_REF]Elgar and Guza, 1985a). The amount of energy transfer between near-resonant components depends on their relative phases, whose information is contained in the biphase [START_REF] Hasselmann | Bispectra of ocean waves[END_REF][START_REF] Kim | Bispectrum and nonlinear wave coupling[END_REF].

β( f 1 , f 2 ) = θ(f 1 ) + θ(f 2 ) -θ(f 1 + f 2 ) (e.g., see
Bispectra were computed on the free surface elevation signals down-sampled to 16 Hz by averaging estimates from 126 Hann-windowed records of 128 seconds, which were overlapping by 75%. Statistical stability was increased by merging estimates over three frequencies (e.g., see Elgar and Guza, 1985a) yielding approximately 205 equivalent degrees of freedom, and a spectral resolution of 0.023 Hz. A validation of the bispectra computation is provided in Appendix B, which compares surface elevation third-order moments computed across the whole wave flume via bispectral (Elgar and Guza, 1985a;[START_REF] Elgar | Relationships involving third moments and bispectra of a harmonic process[END_REF] and statistical definitions (Eq. 1 and 2).

Non-linear energy transfers

In the present study, we are principally interested in the growth of high harmonics due to non-linear energy exchanges between triads and the effect of these forced components on the dispersive properties of shoaling waves. To quantify the relative importance of forced energy at high harmonics, we represent these exchanges with a Boussinesq approximation in a spectral energy balance equation, following a similar approach as developed in [START_REF] Herbers | Spectral energy balance of breaking waves within the surf zone[END_REF]. Such Boussinesq approximation has already demonstrated its potential to predict, in a very computationally efficient manner, the cross-shore evolution of spectral components in the nearshore [START_REF] Herbers | Nonlinear shoaling of directionally spread waves on a beach[END_REF][START_REF] Norheim | Nonlinear evolution of surface wave spectra on a beach[END_REF][START_REF] Herbers | Spectral energy balance of breaking waves within the surf zone[END_REF][START_REF] De Bakker | Nonlinear infragravity-wave interactions on a gently sloping laboratory beach[END_REF][START_REF] Padilla | Transfer and dissipation of energy during wave group propagation on a gentle beach slope[END_REF]. Note that the Boussinesq approximation (detailed below) restricts the present modelling approach to A2 and A3 since conditions during A1 are too dispersive.

For unidirectional waves propagating shoreward on an alongshore-uniform beach and assuming a weak reflection at the shoreline, a balance between the cross-shore gradient of the energy flux spectrum F (f ), a source term S nl (f ) quantifying the non-linear energy exchanges between triads, and a dissipation term S dis (f ) reads (e.g., Eq. 1 of [START_REF] Herbers | Spectral energy balance of breaking waves within the surf zone[END_REF])

∂F (f ) ∂x = S nl (f ) + S dis (f ). ( 9
)
Breaking processes are ignored so that the dissipation term S dis reduces to the energy losses by bottom friction S f r , here simply modelled after [START_REF] Thornton | Transformation of wave height distribution[END_REF]:

S f r (f ) = ρc f 1 6π 2πf sinh |k|h H 3 , ( 10
)
where h is the mean water depth and |k| is the wavenumber modulus obtained from the linear wave dispersion relation (hereafter denoted κ L , see Section 4). The friction coefficient c f was set to 0.0055 after calibration in the deepest section of the wave flume.

Assuming that the wave field is weakly non-linear, weakly dispersive, and that these effects are of similar order (U r ∼ O(1), i.e. Boussinesq regime), the non-linear source term S nl can be approximated with an integral of the bispectrum as follows [START_REF] Herbers | Nonlinear shoaling of directionally spread waves on a beach[END_REF][START_REF] Norheim | Nonlinear evolution of surface wave spectra on a beach[END_REF][START_REF] Herbers | Spectral energy balance of breaking waves within the surf zone[END_REF]:

S nl (f ) = 3πf h ∞ -∞ Im{B * (f , f -f )}df . (11)
Eq. 11 differs from the expression of [START_REF] Herbers | Spectral energy balance of breaking waves within the surf zone[END_REF] (their Eq. 2) in several points: the conjugate of B is taken in order to be consistent with their definition of the bispectrum (conjugate of the present definition); S nl is here defined as a function of f (not ω) and the definition with the full integral is kept (symmetric properties of the bispectrum are not used to decompose it). Previous studies proposed evolution equations for the energy and bispectra (e.g., the stochastic model of [START_REF] Herbers | Nonlinear shoaling of directionally spread waves on a beach[END_REF] to simulate the propagation and transformation of directional seas. As in [START_REF] Herbers | Spectral energy balance of breaking waves within the surf zone[END_REF], we take advantage of the spatial resolution of the present dataset and directly evaluate S nl and S f r at each available cross-shore location using trapezoidal rules for approximating the integrals. Besides removing uncertainties associated with its cross-shore prediction, using measured bispectra has for advantage to relax the limitations of the stochastic Boussinesq model of [START_REF] Herbers | Nonlinear shoaling of directionally spread waves on a beach[END_REF] (such as the distance over which bispectra can be propagated, see [START_REF] Freilich | Nonlinear effects on shoaling surface gravity waves[END_REF]. In order to be consistent with the original equation derived by Herbers and Burton (1997, their Eq. 22a), the energy in Eq. 9 is assumed to propagate at the shallow water wave speed so that F (f ) = ρgE(f ) √ gh, with ρ the water density and g the acceleration of gravity.

After these considerations, Eq. 9 simplifies to

∂E(f ) gh(x) ∂x = 1 ρg (S nl (f ) + S f r (f )) . ( 12
)
Integrating this equation (in space) between the location of the first gauge x 0 and any location x 1 prior to the mean breaking point location yields the following expression for the energy density spectrum at x 1 E(x 1 , f ):

E(x 1 , f ) = E(x 0 , f ) gh(x 0 ) gh(x 1 ) + 1 ρg gh(x 1 ) x1 x0 (S nl + S f r ) dx. (13)
Finally, the ratio x1 x0 S nl dx/F (x 1 ), which represents the energy flux received (or lost) via non-linear coupling between x 0 and x 1 over the total energy flux at x 1 , is used in the following as an approximation of the relative amount of forced energy at x 1 . Outside the surf zone, where dissipative processes dominate, this estimation was found to be more reliable than the bicoherence, which is often used as a proxy for such an estimation but lacks general consensus upon its definition, see Appendix C for more details.

4. Wavenumber and phase velocity spectra of shoaling and surf zone waves

In this section, we present and describe the main results from the cross-spectral and bispectral analyses.

The computation of dominant wavenumber κ and phase velocity c spectra for varying degrees of nonlinearity resulted in the identification of four different regimes of propagation, which broadly consist of: a linear regime (stage I), a shoaling regime relatively far from the mean breaking point (stage II), a shoaling regime near the mean breaking point (stage III) and a surf zone regime (stage IV). The Ursell number U r (Fig. 2e) was used as a reference to define these regimes of propagation, which are characterized by similar patterns -between all wave tests considered here -in wavenumber and phase velocity spectra. For all stages, κ(f ) and c(f ) spectra are compared to predictions from the linear wave dispersion relation, which links the spatial and temporal information of a linear wave field:

ω 2 = g|k| tanh (|k|h). (14)
In the following, the subscript 'L' is used throughout the manuscript to refer to κ or c values that are solution to Eq. 14 (i.e. κ L and c L = ω/κ L ). For conciseness, the propagation as linear waves is discussed using the results from the most non-linear case (A2), already shown in Fig. 3. The other stages of propagation focus on the differences between broad and narrow spectra, i.e. between A2 and A3. Note that κ and c characteristics during A1 are very similar to those obtained during A2 for similar U r numbers. Finally, given the relatively high spatial resolution, frequency-wavenumber power spectra P (ω, k) (e.g., see [START_REF] Redor | Experimental evidence of a hydrodynamic soliton gas[END_REF] were computed using 2D Fourier analysis at four different cross-shore sections corresponding to stages I, II, III and IV, by averaging over 63 Hann-windowed time-records of 128 seconds overlapping by 50% (Fig. 4). The chosen length of the cross-shore section (width of the data window for computing the 2D Fast-Fourier transform) is a best compromise for reaching a proper resolution in wavenumber k and still assuming little change in wave type over the section. The 2D Fourier analysis shown in Fig. 4 provides a qualitative information on energy spreading in (ω, k) space and is particularly useful to illustrate what the single-valued dominant wavenumber represents in a mixed sea-state. Fig. 4 will be used in the present and following sections to analyse the link between non-linear energy exchanges and the dispersive properties of the wave field.

4.1. Propagation as free, linear waves (I)

Fig. 3 shows for the most non-linear case A2 an almost perfect match for f > 0.6 f p between the measured wave phases, wavenumbers and phase velocity spectral estimates with the predictions from the linear wave dispersion relation (Fig. 3a, 3c and 3d respectively). The decrease of coherence at high frequencies (Fig. 3b) can be explained by the relatively low energy content in the spectrum tail (> 86% of the variance is contained between 0.6 f p and 1.5 f p during A2). It also provides an explanation for the slight deviations of measured κ and c from predictions by the linear wave dispersion relation (Fig. 3c-d). However, the 2D Fourier analysis performed over the flat section of the flume (Fig. 4a-b) reveals the presence of energy along the dispersion relation for bound high-frequency components (e.g., see at f > 2.5f p for A2 and around 2f p for A3), suggesting that energy transfers towards high harmonics already occur in this region. Nonetheless, the energy is mainly distributed along the linear dispersion relation, which indicates that short-wave components were mostly propagating as free waves in the deepest regions of the wave flume

(U r ∼ O(0.01 -0.1)).
Data points in the infragravity frequency band (f < 0.6 f p ) correspond to wave components generated via difference interactions and that are bound to short-wave groups. As such, bound infragravity waves propagate at a speed close to the corresponding group velocity, which is lower than their intrinsic phase velocity. This is confirmed in Fig. 3d for A3, which explains the significant growth of the third harmonic 3f p observed for that test (Fig. 5b).

Consistent with field observations (e.g. Elgar and Guza, 1985a;[START_REF] Norheim | Nonlinear evolution of surface wave spectra on a beach[END_REF], relatively strong nonlinear exchanges by sum interactions are observed during A2, despite the broader spectrum conditions (note the scale difference between the A2 and A3 cases). Short-wave frequencies distant from f p (e.g. ∼ 1.5f p )

show a strong coupling with f p , which is the consequence of a broader spectrum for A2 compared to A3.

A striking result in this regime of propagation is the fact that, despite non-linearities becoming important ( ∼ 0.12 -0.2, S k ∼ 1), the principal wave components follow the linear wave dispersion relation. The spatial structure of the wave field for f ∈ [0.6 f p ; 1.5 f p ] is therefore well-described by the linear wave theory (Fig. 5c and5d). However, this is not the case for f > 1.5 f p , where the overestimation of wavenumbers by the linear dispersion relation increases with f (up to a factor 2.5 at 3.5 f p ). As the biphase at both (f p , f p ) and (f p , 2f p ) are close to 0, this regime of propagation is consistent with Stokes-like non-linearities (e.g., see Elgar and Guza, 1985a).

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

Energy transfers towards the infragravity band concentrate at frequencies around 0.1-0.2f p (Fig. 5a and5b) and principally originate from strong difference interactions which transfer energy from f p towards that infragravity frequency and components at frequency slightly lower than f p (see the negative imaginary part of B * along the f p anti-diagonal, Fig. 5g and5h).

Shoaling of non-linear waves in the vicinity of the breaking point (III)

The next regime analysed in Fig. 6 corresponds to non-linear waves approaching the mean breaking point location (U r ∼ 0.7). At this Ursell number, wavenumber and phase velocity spectra become less frequency-dependent at high frequencies. This can be explained by the fact that most of the energy at those frequencies is bound to primary components in the range [0.6 f p ; 1.5 f p ], and that c values display little variations at these frequencies (c/c L (f p ) is within 5% of c(f p )/c L (f p )). As a result, wavenumbers at high frequencies are relatively well described by a simple dispersion relation given by κ(f ) = 2πf /c(f p ) (Fig. 6c-d). Note that at these depths (κ p h ∼ 0.57 and 0.47 for A2 and A3 respectively), this approximately corresponds to the shallow water dispersion relation κ sw (f ) = 2πf / √ gh, which explains why most of the energy is spread around this dispersion relation in Fig 4e-f. Non-linear amplitude effects become increasingly important ( = 0.25 and S k = 1.4-1.5) and deviations between c(f p ) and predictions by the linear dispersion relation are observed (Fig. 6e-f). As noted and observed by [START_REF] Herbers | Nonlinear dispersion of surface gravity waves in shallow water[END_REF], this is the result of nonlinear interactions, which alter the wavenumber of all three wave components involved. Indeed, bound high harmonics contribute to an increase in skewness and height of principal wave crests, which has for effect to enhance their propagation speed. Variations in the strength of non-linear interactions thus explain the varying magnitude of the deviations of c(f p ) from predictions by the linear dispersion relation as observed in both field and laboratory settings [START_REF] Thornton | Energy saturation and phase speeds measured on a natural beach[END_REF][START_REF] Herbers | Nonlinear dispersion of surface gravity waves in shallow water[END_REF][START_REF] Tissier | Field measurements and non-linear prediction of wave celerity in the surf zone[END_REF].

For this range of Ursell number, overestimations of κ by the linear dispersion relation remain close to a factor 2 at f = 3.5 f p .

In these regions of the wave flume, surface elevation bispectra are still dominated by real values (note the x2 in imaginary values, Fig. 6g-h), typical of (horizontally) skewed and nearly symmetric waves (As ∼ 0, see Fig. 2c-d and also [START_REF] Masuda | A note on the imaginary part of bispectra. Deep Sea Research Part A[END_REF]Elgar and Guza, 1985a). The strong coupling between first and second harmonics, which for instance explains the energy peak at 3f p in A3 (Fig. 6b) is here evident in the imaginary part of the bispectra at (f p ,2f p ) (see also the bicoherence in Fig. C1). Infragravity frequencies between 0.1-0.2f p keep receiving energy via relatively strong difference interactions. As a consequence, prior to short-wave breaking infragravity waves are still bound to short-wave groups during both tests (see Fig. 6e-f, and also [START_REF] De Bakker | Nonlinear infragravity-wave interactions on a gently sloping laboratory beach[END_REF].

Surf zone waves (IV)

In the surf zone (U r taken at 3.2), wavenumber and phase velocity spectra are frequency-independent (Fig. 7e-f). The large differences observed between the measured phase velocity at peak frequency and that predicted with the linear dispersion relation demonstrate the dominance of amplitude effects over dispersive ones [START_REF] Herbers | Nonlinear dispersion of surface gravity waves in shallow water[END_REF]. All components travel at a speed slightly larger than that of shallow water waves (∼ 10 -20% depending on the wave test), which is consistent with many past field and laboratory observations of wave dispersion in the surf zone [START_REF] Inman | Mixing in the surf zone[END_REF][START_REF] Svendsen | Wave characteristics in the surf zone[END_REF][START_REF] Stive | Velocity and pressure field of spilling breakers[END_REF][START_REF] Thornton | Energy saturation and phase speeds measured on a natural beach[END_REF], and many others). For both A2 and A3 tests, the relation c = gh(1 + ) is shown in Fig. 7e-f (see also Fig 4g-h, for the corresponding frequency-wavenumber surface elevation power spectra). This formula is asymptotic in shallow water to the non-linear dispersion relation given by [START_REF] Hedges | An empirical modification to linear wave theory[END_REF] with the adjustment of [START_REF] Booij | Gravity waves on water with non-uniform depth and current[END_REF] (see also [START_REF] Catalán | Remote sensing of breaking wave phase speeds with application to non-linear depth inversions[END_REF], for an assessment of these relations for depth-inversion applications). In the inner surf zone, where ∼ 0.3 during all tests (Fig. 2b), this yields c ∼ 1.14 √ gh, a relation also found in the field by [START_REF] Tissier | Field measurements and non-linear prediction of wave celerity in the surf zone[END_REF]. This relation appears as an upper bound, in an average sense (e.g. over several wave groups), for the propagation speed of wave components in the inner surf zone.

In the short-wave frequency range, bispectra display relatively weak and nearly equal real and imaginary parts, characteristics of asymmetric (pitched forward) broken waves, with biphases approaching -π/2 (Elgar and Guza, 1985a;[START_REF] Doering | Skewness in the nearshore zone: A comparison of estimates from marshmcbirney current meters and colocated pressure sensors[END_REF]: for instance during A3, β(f p , f p ) ∼ β(f p , 2f p ) ∼ -π/3. Since these values were still close to 0 during stage III, this indicates that bound high-harmonic components (here only shown for 2f p and 3f p ) slowly drift out of phase from the principal component starting from the mean breaking point location. This is in contrast with short-wave group-forced bound infragravity waves that slowly drift out of phase with the short-wave envelope in the shoaling region (e.g., see Elgar and Guza, 1985a). Since this process occurs over long distances, the differences in phase velocities between f p , 2f p and 3f p are not sufficiently large to be observed in the present dataset. Finally, in the inner surf zone, intense and complex energy transfers occur within the infragravity band and also between infragravity and short-wave components. For more information on these processes during the GLOBEX experiments, the reader is referred to the study of [START_REF] De Bakker | Nonlinear infragravity-wave interactions on a gently sloping laboratory beach[END_REF].

Role of non-linear energy transfers on κ and c

Most differences observed in wavenumber and phase velocity spectra between broad and narrow-banded wave conditions concentrate in the shoaling region (stages II and III). In stage II, c values at 2f p and 3f p during A2 lie between c(f p ) and the values predicted by linear wave theory while c(f p ), c(2f p ) and c(3f p ) are all equal during A3 (Fig. 5e-f). As these differences between broad versus narrow-band spectra are likely explained by the relative importance of forced energy at those frequencies (see the cross-spectral analysis performed on synthetic data and discussed in Appendix A and also Herbers and Guza, 1992;[START_REF] Herbers | Spectral energy balance of breaking waves within the surf zone[END_REF], we analyse non-linear energy transfers to 2f p and 3f p in more details here.

Wave amplitudes at f p , 2f p and 3f p (a fp , a 2fp and a 3fp respectively) computed from energy spectra modelled with Eq. 13 are compared against observations in Fig. 8. Overall, the Boussinesq approach of [START_REF] Herbers | Nonlinear shoaling of directionally spread waves on a beach[END_REF] for the non-linear energy transfers between triads accurately predicts the growth of second and third harmonics across the shoaling zone for both broad and narrow-band wave tests.

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

A closer look at the cross-shore evolution of the source terms (Fig. 9) shows that friction effects are negligible at high frequencies but not around the peak frequency where the energy dissipated through bottom friction is of similar order than the energy lost via non-linear coupling. During A3 (Fig. 9b), the steady increase of S nl (2f p ) indicates a gradual growth of a 2fp in the shoaling region (see also Fig. 8b). This growth is most intense right before the mean breaking point location, around the location corresponding to stage III. Computations of the ratio respectively [START_REF] Kim | Digital Bispectral Analysis and Its Applications to Nonlinear Wave Interactions[END_REF].

The effect of varying relative amounts of forced energy not only explain the differences observed in κ and c spectra at high frequencies between A2 and A3 but also their variation across frequencies. During A3 for instance, non-linear energy transfers in the short-wave frequency band were predominantly towards 2f p and 3f p (Fig. 5h and6h). At these harmonics, energy is predominantly forced (> 70% around 2f p and > 90% around 3f p ), thus, the dominant κ and corresponding c values are directly related to the values at the peak frequency. In contrast, frequencies located in valleys between harmonics (i.e. 1.5f p and 2.5f p ) receive very little energy through the coupling of triads. This eventually leads to the patterns in c observed in Fig. 5f, which are very similar to those obtained by [START_REF] Crawford | Effects of nonlinearity and spectral bandwidth on the dispersion relation and component phase speeds of surface gravity waves[END_REF] with their model for modulated wave trains. During A2, lower relative amounts of forced energy are found at high harmonics and dominant wavenumbers lie between κ(f p ) and the intrinsic value at that frequency predicted by the linear dispersion relation. Furthermore, due to the broader spectrum conditions, this forced energy is more spread across frequencies, which explains the weaker variations of κ and c across frequencies as compared to A3 (e.g., see Fig. 5e-f).

The present study considers unidirectional wave fields transforming across a mildly sloping beach. A typical field situation where incident wave spectra exhibit varying degrees of directional spreading is unclear as far as non-linear energy transfers are concerned. According to Boussinesq theory (e.g., [START_REF] Herbers | Nonlinear shoaling of directionally spread waves on a beach[END_REF], directional spreading has only a weak influence on the efficiency of non-linear energy transfers by sum interactions, as opposed to difference interactions. For instance, these authors predicted a weaker growth of high-harmonic bound wave by 10-20% for very large directional spreads (60 • ). For more realistic spreading angles in the nearshore, de Wit et al. ( 2020) also recently noted variations of approximately 10%.

Although small, such a decrease in efficiency of non-linear energy transfer towards high harmonics due to an increasing directional spreading might lead to lower relative amounts of forced energy, subsequently leading to slightly less important deviations of κ and c spectral estimates from predictions by the linear wave dispersion relation (see also Appendix A). The larger spectral bandwidth generally associated with wave spectra exhibiting a large directional spreading is, however, likely to be the principal reason explaining weaker couplings between triads and thus weaker relative amounts of forced energy at high harmonics.

Discussion and concluding remarks

Cross-spectral and bispectral analyses were employed on a highly-resolved surface elevation dataset to study the dispersive properties of waves shoaling and breaking over a mildly sloping beach. For all wave tests considered here, four regimes of propagation (I to IV) with specific characteristics in dominant wavenumber and phase velocity spectra could be defined using a local Ursell number. Stage II (U r ∼ 0.3) is particularly interesting as it simultaneously shows significant non-linear effects at high harmonics (typically f > 1.5f p ), which are evidenced by the large deviations of κ and c spectra from predictions by the linear wave dispersion relation, and a clearly linear spatial field for the primary components. Although less energy was imposed during A3 as compared to A2, the narrow-banded conditions promoted relatively more important non-linear energy transfers towards high harmonics such that phase velocity at 2f p and 3f p were found equal to those at the peak frequency. For A2, with a broader spectrum imposed, forced energy at high frequencies was found in relatively smaller proportion and c spectral estimates lied between values at the peak frequency and that given by the linear dispersion relation. Closer to the mean breaking point (stage III), non-linear energy transfers were found more intense, which, together with the less dispersive conditions compared to stage II explain that κ and c spectra become less frequency dependent. In the surf zone (stage IV), wavenumber and phase velocity spectra essentially showed the same characteristics for all tests and a modified shallow water dispersion that accounts for non-linearity (c = gh(1 + )) accurately describes the dispersive properties of the wave field. Some of the results presented here for shoaling and breaking waves confirm previous findings obtained in the field [START_REF] Thornton | Energy saturation and phase speeds measured on a natural beach[END_REF]Elgar and Guza, 1985b;Herbers et al., 2002, among others), yet the unprecedented level of accuracy obtained with the GLOBEX dataset allowed to highlight and quantify the effect of varying levels of forced energy at high harmonics on the dispersive properties of the wave field.

The detailed analysis of the different regimes of propagation highlighted here indirectly provides an assessment of the validity of the linear wave dispersion relation at several locations in the nearshore and for a range of wave conditions. Two specific regions in the short-wave frequency band can be defined: primary components (f ∈ [0.6f p , 1.5f p ]) and high harmonics (f > 1.5f p ). For the primary components, deviations from the linear wave dispersion relation are due to non-linear amplitude effects and concentrate in a region near the mean breaking point and in the surf zone. These effects can cause deviations of κ and c from predictions by the linear wave dispersion of the order of 10-20%. At high frequencies, deviations are much larger (O(100%)) and also occur much farther seaward. For instance, stage II corresponds to κ p h = 0.71, 0.70 and 0.62 for wave tests A1, A2 and A3 respectively. The fact that such deviations from predictions by the linear dispersion relation are found in κ and c spectra so far seaward of the surf zone for different wave conditions is expected to have strong implications for a wide range of applications, two of which are briefly discussed below.

Depth-inversion algorithms such as cBathy [START_REF] Holman | cBathy: A robust algorithm for estimating nearshore bathymetry[END_REF] estimate the water depth by applying the linear wave dispersion relation to the most coherent pairs of wave frequency and wavenumber extracted from timeseries of video images [START_REF] Stockdon | Estimation of wave phase speed and nearshore bathymetry from video imagery[END_REF][START_REF] Plant | Ocean wavenumber estimation from wave-resolving time series imagery[END_REF]. However, due to the working principles of video cameras, video data can be more coherent at high harmonics and frequencies greater than 1.5 f p are often picked to invert the water depth [START_REF] Stockdon | Estimation of wave phase speed and nearshore bathymetry from video imagery[END_REF][START_REF] Bergsma | Video-based depth inversion techniques, a method comparison with synthetic cases[END_REF]. In such cases, large deviations of the dominant wave phase velocity from linear wave theory are expected, which likely explains why [START_REF] Brodie | Evaluation of video-based linear depth inversion performance and applications using altimeters and hydrographic surveys in a wide range of environmental conditions[END_REF] noted an overestimation in the detected phase velocities by 20-100% up to 50 m seaward of the mean breaking point location. In the absence of knowledge on the relative amounts of forced energy at high harmonics, pairs of frequencies and wavenumbers around the peak frequency should be preferred in regions where non-linearities are expected. Similarly, the large overestimation of κ by the linear dispersion relation causes the widely reported blow-up at high harmonics when correcting sub-surface pressure signals for depth attenuation [START_REF] Bonneton | Recovering water wave elevation from pressure measurements[END_REF][START_REF] Mouragues | Field data-based evaluation of methods for recovering surface wave elevation from pressure measurements[END_REF]Martins et al., 2020b). The use of a cutoff frequency prevents for instance the accurate description of third-order parameters and wave height distributions (e.g., see Martins et al., 2020a). In shallow water (µ 0.3), a weakly dispersive formula recently developed allows for an accurate correction of the energy levels at high harmonics both seaward of [START_REF] Bonneton | A nonlinear weakly dispersive method for recovering the elevation of irrotational surface waves from pressure measurements[END_REF][START_REF] Mouragues | Field data-based evaluation of methods for recovering surface wave elevation from pressure measurements[END_REF] and in the surf zone (Martins et al., 2020b). However in deeper water, a non-linear fully dispersive reconstruction requires knowledge on the spatial structure of the wave field, which is generally lacking.

thus validating the present computations of the bispectrum. Note that for these comparisons, blocks were not tapered and rectangular windows were used for computing the Fast Fourier Transforms. Using any other types of windows resulted in small differences between the two definitions of both third-order wave parameters.

Appendix C. Normalised bispectra (bicoherence)

We generally seek a normalisation of the bispectrum so that it takes 1 as value when there is a full coupling between the three components involved and 0 when there is none. In the case of surface gravity waves, this would correspond to a situation with the only presence of forced or free energy respectively.

At present, there does not seem to be a consensus in the literature on what definition of the bicoherence should be used to robustly quantify or at least characterize the relative dominance of forced energy (e.g. [START_REF] Haubrich | Earth noise, 5 to 500 millicycles per second: 1. spectral stationarity, normality, and nonlinearity[END_REF][START_REF] Kim | Digital Bispectral Analysis and Its Applications to Nonlinear Wave Interactions[END_REF]Elgar and Guza, 1985a;[START_REF] Hinich | Normalizing bispectra[END_REF]de Bakker et al., 2015, and many others). The normalisation proposed by [START_REF] Haubrich | Earth noise, 5 to 500 millicycles per second: 1. spectral stationarity, normality, and nonlinearity[END_REF] is frequently used for studying non-linearities in the nearshore area but was shown to lead to values greater than 2 (e.g. Elgar and Guza, 1985a). In the present study, this was also the case, but only when merging across frequencies was performed beforehand. Nonetheless, a slightly different normalisation that was proposed in [START_REF] Hagihira | Practical issues in bispectral analysis of electroencephalographic signals[END_REF] is used in the present case as it seems the most appropriate normalisation (see their Appendix):

b( f 1 , f 2 ) = |B(f 1 , f 2 )| E [|A(f 1 )| • |A(f 2 )| • |A * (f 1 + f 2 )|] (C.1)
As explained in [START_REF] Hagihira | Practical issues in bispectral analysis of electroencephalographic signals[END_REF] In the absence of forced components, the forcing corresponds to a JONSWAP spectrum with parameter similar to A3 (see Table 1).

  with a good correspondence between measured c values with shortwave groups propagation speed (c env , shown as the orange line), estimated via cross-spectral analysis of short-wave groups envelop. The latter was computed following[START_REF] Janssen | Long waves induced by short-wave groups over a sloping bottom[END_REF], as the absolute value of the analytical signal of ζ: |ζ + iH{ζ}| (no low-pass filtering is needed for the present application). Note that the discrepancies observed at f ∈ [0.3 f p ; 0.6 f p ] can be explained by the lower levels of energy at these frequencies:E(f ) ∼ O(10 -5 ) m 2 /Hz as opposed to E(f ) ∼ O(10 -4 ) m 2 /Hz around f = 0.1 f p .4.2. Shoaling of weakly non-linear waves (II)At higher Ursell number (U r ∼ 0.3), wave non-linearities become significant and non-linear energy transfers to high harmonics promote the dominance of bound high harmonics (Fig.4c-d). This is observed in the considerable deviation of dominant wavenumbers and phase velocity spectral estimates from the predictions of the linear dispersion relation for f > 1.5 f p (see Fig.5c-f). During A3, the phase velocity values at 2f p and 3f p are equal to that of the peak frequency (difference < 1%, c(f p ) = 1.95 m/s), while during A2 they lie between the value at f p and the predictions of the linear dispersion relation (∼ 7% difference from c(f p ) = 2.18 m/s). At the second harmonic 2f p , most of the non-linear energy transfers occur via sum interactions of components around f p (self-self interactions) as evidenced by the positive imaginary part of B * around (f p , f p ), see Fig.5gand 5h for tests A2 and A3 respectively. Although less evident in the bispectra, the normalised bispectra shown in Fig.C1indicate that strong coupling involving f p and higher harmonics also exist in the shoaling region. These are particularly strong around (2f p , f p )

Figure 1 :Figure 2 :

 12 Figure 1: Elevation z of the 1:80 concrete beach against the cross-shore distance x in the Scheldegoot flume during the GLOBEX project. The wave paddle is located at x = 0 m. The grey '+' symbols show the position of the wave gauges across the wave flume.

Figure 3 :

 3 Figure3: Phase (a), coherence-squared (b), dimensionless wavenumber (c) and phase velocity (d) spectra computed at x = 10.65 m for case A2. 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 :

 4 Figure4: 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 + )).

Figure 5 :

 5 Figure 5: 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 6 :

 6 Figure 6: Similar to figure 5, but with Ur ∼ 0.7 (stage III).

Figure 7 :

 7 Figure 7: Similar to figure 5, but with Ur ∼ 3.2 (stage IV). Additionally, the green dashed lines in panels e) and f) correspond to the modified shallow water dispersion relation c = gh(1 + ).

Figure 8 :Figure 9 :

 89 Figure8: Comparison between observed and predicted (Eq. 13) wave amplitudes at fp, 2fp and 3fp (a fp , a 2fp and a 3fp respectively). The cross-shore locations corresponding to the regimes of propagation discussed in section 4 are indicated by the vertical dashed lines for both wave tests. Since the present modelling approach neglects wave breaking-induced energy dissipation, modelled amplitudes are shown only until the approximate mean breaking point location (slightly after stage III). The amplitude growth predicted with Green's law is also shown.

Figure A1 :

 A1 FigureA1: 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 Table1).

Figure A2 :Figure

 A2 Figure A2: Wavenumber (a) and phase velocity (b) spectra computed on the synthetic surface elevation timeseries for varying ratio of forced to free second harmonic amplitudes.

Figure C1 :

 C1 Figure C1: Squared bicoherence b 2 computed using the definition of Hagihira et al. (2001) for wave test A2 (left) and A3 (right). The octant above the 1:1 diagonal corresponds to the bicoherence for stage II while stage III is that under the diagonal. Due to the symmetric properties of the bispectrum, only one octant is shown for each stage. The separation between infragravity and short-wave frequencies (0.6 fp) is shown as the dashed black line. Only bicoherence values greater than the 95% significance level as defined by Haubrich (1965) (b 2 95% ≥ 6/d.o.f., with d.o.f. the equivalent number of degrees of freedom) are shown.

  x1x0 S nl dx/F (x 1 ) at the location x 1 corresponding to stage II (resp. III) suggest that approximately 70% (resp. 80%) of the energy at 2f p is forced during A3. At 3f p , these estimations range between 90% and 100% during stage II and III respectively. Although the resolution in wavenumbers does not allow a clear separation of forced and free energy at 2f p , the frequency-wavenumber power spectra shown in Fig.4d(II) and 4f (III) corroborate these numbers and clearly indicate that most of the energy at 3f p is forced. In contrast, non-linear energy exchanges between triads display a more complex picture for the broader conditions of A2. S nl (2f p ) oscillates around 0 so that when spatially integrated, non-linear energy exchanges explain only 20% (30%) and 30% (50%) of the total energy at 2f p and 3f

p respectively during stage II (III). These numbers are quite consistent with the bicoherence values shown in Appendix C. For both wave tests, b 2 values are found to vary little between stage II and III: for A2, b 2 (f p , f p ) ∼ 0.3 and b 2 (f p , 2f p ) ∼ 0.2, while these values oscillate between 0.55 and 0.65 for A3. b 2 (f p , f p ) and b 2 (f p , 2f p ) can be considered as crude estimates of the relative amount of forced energy at 2f p and 3f p

Table 1 :

 1 , the bicoherence is 1 only when phase angles of all triple products 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.

	are equal.
	Squared bicoherence values b 2 for both II and III regimes of propagation are shown in Fig. C1. As
	frequency merging is inappropriate with this definition of the bicoherence (due to the expected value in the
	denominator), b values were averaged over a 3x3 square in the frequency space.

Acknowledgments

Kévin Martins greatly acknowledges the financial support from the University of Bordeaux, through an International Postdoctoral Grant (Idex, nb. 1024R-5030). The GLOBEX project was supported by the European Community's Seventh Framework Programme through the Hydralab IV project, EC Contract 261520. The GLOBEX data used in this research can be accessed on Zenodo at https://zenodo.org/record/4009405 or from the authors, and can be used under the Creative Commons Attribution 4.0 International licence.

The bispectral analysis tools developed and used in this study are accessible from the first author github repository at https://github.com/ke-martins/bispectral-analysis. Fruitful exchanges with Steve Elgar regarding the bispectral analysis were greatly appreciated. The authors thank two anonymous reviewers for their helpful comments and suggestions.

Appendix A. Sensitivity analysis of the cross-spectral technique to the presence of forced energy

The situation where both free and forced components exist in a wave field is not entirely clear as far as the cross-spectral analysis is concerned. Some studies consider the cross-spectral estimates to be biased towards bound harmonics (e.g., see [START_REF] Lake | A new model for nonlinear wind waves. Part 1. Physical model and experimental evidence[END_REF][START_REF] Thornton | Energy saturation and phase speeds measured on a natural beach[END_REF], arguing that the travelling distance between two gauges is reduced for bound components compared to free ones. This is evidently exacerbated when using sub-surface pressure sensors or current meters since bound harmonics are much less attenuated and will dominate the spectrum at some depth.

To analyse the sensitivity of the cross-spectral analysis to different relative levels of forced energy, some tests are performed here on synthetic data. Starting at a position x 1 , a surface elevation timeseries following a JONSWAP spectrum with similar parameters to the A3 test is generated (Fig. A1). A Gaussian-shaped perturbation around 2 f p is added to represent forced harmonics: for any f around 2 f p , the forced harmonic wavenumber is set at κ f orced = 2κ(f /2), so that we have c f orced (f ) = c(f /2). For the phases of forced harmonics, we choose to impose θ f orced (f ) = θ(f /2) for simplicity. The maximum ratio of forced to free harmonic amplitudes a 2fp, f orced /a 2fp, f ree is set to 6, which corresponds to the maximum ratio of energy E(f p )/E(2f p ) of approximately 0.4 that was observed in the flume outside the surf zone. Each wave component is then propagated using its corresponding wavenumber to a second location x 2 .

The results from the cross-spectral analysis performed on the synthetic timeseries at x 1 and x 2 are provided 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 .

In the absence of forced components (a 2fp, f orced /a 2fp, f ree = 0), components around 2 f p follow the dispersion relation given by linear wave theory, i.e. κ(2 f p ) = κ L (2 f p ) and c(2 f p ) = c L (2f p ). Interestingly, when forced and free components are found in equal proportion (a 2fp, f orced /a 2fp, f ree = 1), the relation

2 is verified. For higher ratio of forced to free energy, κ(2 f p ) and c(2 f p ) values rapidly converge towards the values corresponding to the peak frequency. Overall, these results suggest that the cross-spectral analysis is not biased towards forced components, but simply provides a dominant wavenumber (or averaged in terms of energy) and its corresponding phase velocity.

Appendix B. Validation of the bispectrum computations using the skewness and asymmetry test

The wave skewness and asymmetry are third-order moments characterizing the wave shape. These parameters are directly related to the energy content at high frequencies of the spectra and bispectra [START_REF] Elgar | Relationships involving third moments and bispectra of a harmonic process[END_REF]. The fact that the statistical (see Eq. 1 and 2) and the bispectrum-based definitions are theoretically equivalent provide a means to validate the bispectrum's calculations.

Due to its symmetry properties, the bispectrum can be uniquely defined in a single octant in the frequency-space [START_REF] Hasselmann | Bispectra of ocean waves[END_REF]Elgar and Guza, 1985a). If we denote by f N the Nyquist frequency, the octant with positive frequencies is bounded by the vertices at (0,0), (f N/2 , f N/2 ) and (f N ,0).

The wave skewness and asymmetry are defined by the sum of the real and imaginary parts of the bispectrum over this octant respectively, normalized by the cube of the free surface elevation standard deviation [START_REF] Elgar | Relationships involving third moments and bispectra of a harmonic process[END_REF]Elgar and Guza, 1985a). This reads: with n > l and n + l < N .

The cross-comparison between statistical and bispectrum-based definitions of the wave skewness and asymmetry computed for A3 is provided in Fig. B1 and shows a perfect match between both definitions,