The relationship between sea-swell bound wave height and wave shape

: The nonlinear wave shape, expressed by skewness and asymmetry, can be calculated from surface elevation or pressure time series using bispectral analysis. Here, it is shown that the same analysis technique can be used to calculate the bound superharmonic wave height. Using measured near-bed pressures from three different ﬁeld experiments, it is demonstrated that there is a clear relationship between this bound wave height and the nonlinear wave shape, independent of the measurement time and location. This implies that knowledge on the spatially varying bound wave height can be used to improve wave shape-induced sediment transport predictions. Given the frequency-directional sea-swell wave spectrum, the bound wave height can be predicted using second order wave theory. This paper shows that in relatively deep water, where conditions are not too nonlinear, this theory can accurately predict the bispectrally estimated bound superharmonic wave height. However, in relatively shallow water, the mismatch between observed and predicted bound wave height increases signiﬁcantly due to wave breaking, strong currents, and increased wave nonlinearity. These processes are often included in phase-averaged wind-wave models that predict the evolution of the frequency-directional spectrum over variable bathymetry through source terms in a wave action balance, including the transfer of energy to bound super harmonics. The possibility to calculate and compare with the observed bound super harmonic wave height opens the door to improved model predictions of the bound wave height, nonlinear wave shape and associated sediment transport in large-scale morphodynamic models at low additional computational cost. the difference in the spatial evolution of H m 0 and H b , obs observed already for the single data burst presented above. It shows that for the energetic cases (Group 2), the maximum value of H b , obs is systematically found in the area where H m 0 is already decaying due to breaking.


Introduction
Coastal management decisions, such as nourishment strategies and sea level rise scenarios, rely more and more on morphodynamic model simulations. Within these simulations, fluxes in sediment transport, caused by hydrodynamic forcing mechanisms, result in changes in the bathymetry. An important contribution to the sediment transport fluxes is the wave shape-induced sediment transport driven by the skewness and asymmetry of the individual waves [1][2][3][4]. Although its instantaneous magnitude is often smaller than other contributions, it can have a considerable net effect on the bathymetric evolution as the contribution is typically in the dominant wave direction [5,6]. As such, it is important for beach recovery after storm impact [7,8], onshore bar motion [7,[9][10][11] and the evolution of ebb-tidal shoals (see, e.g., in [12,13]).
Current large-scale morphodynamic modeling approaches generally combine a spectral wave transformation model [14][15][16][17][18] and a flow model (see, e.g., in [19][20][21][22][23]) to predict the local wave, flow and sediment transport conditions (see, e.g., in [19,24]). Using a local parameterization based on the wave height, wave period and water depth the wave skewness, asymmetry, and associated sediment transport are obtained (see, e.g., in [25][26][27][28][29]). However, as was shown by Rocha et al. [30], Rocha et al. [31], and De Wit et al. [32], predicting the wave shape using a local approach has its limitations, related to the fact that the prior evolution of wave shape is not taken into account. As a result, the wave shape can be different although the local wave height, period, and water depth are exactly equal, if, for instance, the bed slope is different [33,34], the conditions are rapidly changing [32], or the offshore wave steepness is different [31]. Thus, there is a need for a better way to predict the wave shape that includes the history of the waves before reaching a certain location.
The wave skewness and asymmetry can be computed with a bispectral analysis corresponding to the sum of the real and imaginary parts of the bispectrum, respectively (see, e.g., in [35,36]). The bispectrum is a reflection of the coupling between the primary waves and the bound super and sub harmonics [36,37]. This implies that there is a close connection between the nonlinear wave shape and the proportion of bound wave energy.
The bound portion of energy in the super harmonics within a directionally spread sea-swell wave field can be predicted with the second order theory of Hasselmann [38] based on a local equilibrium over a horizontal bed. However, in the presence of a variable bathymetry and thus spatially evolving sea-swell wave field this may lead to an erroneous estimate as demonstrated by Herbers and Burton [39]. On the other hand, spectral wind-wave models often include a source term to describe the transfer of wave energy from the primary wind-waves to bound super harmonics through triad sum interactions over variable bathymetry [37,40,41]. The modeled bound fraction of superharmonic wave energy is an integration of the source term in the down-wave direction showing up as an additional spectral peak at twice the primary frequencies (see, e.g., in [37,42]). However, to speed up the computations to enable morphodynamic computations at realistic time scales, the phase information is ignored and even though the spatially evolving fraction of bound energy is implicitly predicted, the accompanying skewness and asymmetry are not known.
Examining the three-dimensional (3D) wavenumber-frequency spectrum is a relatively straightforward way to discriminate between the bound and free wave energy as these follow different dispersion relations (see, e.g., in [43]). However, estimating the full 3D wave spectrum requires high-resolution spatial information that is rarely available in the field (see, e.g., in [44] for one of the exceptions). Alternatively, bispectra can be used to characterize the portion of bound energy in a given frequency range. Most efforts to quantify and analyze bound harmonic energy have focussed on the sub-harmonic range (see, e.g., in [45][46][47][48][49][50][51]) following the work of Herbers et al. [52] who demonstrated that the bound fraction of subharmonic (i.e., infragravity) energy could be obtained from the difference interactions in the bispectrum. Significantly less attention has been devoted to quantifying the bound energy in the super harmonic range, with the most notable contributions being the work by Herbers and Guza [53,54] and Herbers et al. [55] who examined bound wave energy in intermediate water depths.
The aforementioned studies [53][54][55] showed that triad sum interactions between wave components with large difference angle of propagation can contribute significantly to the bound near-bed pressure variance at these depths. Interestingly, these are typically associated to negative interaction coefficients according to the theory of Hasselmann [38], while the more classical sum interactions between wave components with small difference in angle of propagation have a positive contribution. Thus, for a given sum frequency in the superharmonic range, both positive and negative contributions from primary wave pairs can occur such as the bispectrum is expected to yield a lower limit of the bound super harmonic energy. Several authors [36,56,57] additionally mention, based on the work of McComas and Briscoe [58], that estimating the bound super harmonics from the bispectrum in a broad-banded spectrum is not straightforward. This inhibits a direct comparison with the predictions of the wind-wave spectral models. Notably, Herbers et al. [55] did find a good match between predictions by the theory of Hasselmann [38] and observations in a case of narrow-banded energetic swell conditions (their Figure 10d) in contrast to conditions with crossing sea states (their Figure 10a-c). This raises the question to what extent the bispectral estimate of the bound super harmonic fraction can work for sea-swell conditions.
In the following, we therefore construct a method to first estimate the bound portion of the energy in the super harmonics in a realistic directionally spread wave field and secondly to use this as a predictor of the nonlinear wave shape controlling wave skewness and asymmetry. To that end, the velocity and pressure data obtained at nine locations on the Ameland ebb-tidal delta from the CoastalGenesis2/SEAWAD field campaign in September 2017 are examined using bispectral analysis. The bound superharmonic fraction is expressed as an equivalent observed bound wave height that is compared with the predicted bound wave height obtained from the equilibrium theory of Hasselmann [38] to explore its spatial evolution. Next, the correspondence between the bound wave height and nonlinear wave shape is examined to explore the potential of using a wave shape parameterization based on the predicted bound wave height instead of a local parameterization. This is followed up with a discussion on the general applicability of such an approach and the necessary steps in spectral wave modeling to enable these predictions.

The Spectrum
The surface elevation is represented as a summation of discrete frequencies as in which A( f m ) is the complex amplitude at discrete frequency f m = m∆ f with ∆ f being the frequency resolution, A * ( f m ) indicates the complex conjugate of A( f m ), i is the imaginary number, and t is time.
The number of discrete spectral estimates is 2N + 1, which are bound by the Nyquist frequencies: ± f N = ± f s /2, in which f s is the discrete sampling frequency of the surface elevation time series. The complex amplitudes are obtained by applying a discrete fast Fourier transformation on the surface elevation. Subsequently, the (double-sided) variance spectrum 1 is defined as in which E[...] denotes the expected value. If besides pressure, also collocated x-and y-velocities are present, the 2D frequency-directional variance spectrum can be computed: in which Θ( f m , θ) is the directional distribution over the discrete number of directional bins N θ resulting in a directional resolution ∆θ = 2π/N θ . Second-order statistics are derived from the variance spectrum, such as the sea-swell significant wave height: in which m min corresponds to the index of the first discrete frequency in the sea-swell range. The factor 2 in Equation (4) arises from the fact that the double-sided variance spectrum is only summed over the positive frequencies. 1 Note that different conventions are found in literature. In this paper, the double-sided spectrum is presented, because the bispectrum (Equation (5)) is commonly also presented in a double-sided form and we want to prevent mixing up singleand double-sided spectra. Furthermore, for readability purposes, we present all equations as a function of the variance instead of the variance density. The variance density spectrum can be obtained by dividing the variance spectrum by ∆ f .

The Bispectrum
The bispectrum is a spectral representation of third-order statistics that can be used to analyze the nonlinear interactions between a triad of frequencies f m , f n and f p that satisfies f m + f n = f p . The discrete variance bispectrum is defined as If the three components are statistically independent, there is no phase correlation and B( f m , f n ) = 0. In that case, the third component f p is not bound to f m and f n but freely propagating. On the other hand, a non-zero bispectrum B( f m , f n ) indicates that (part of) the variance at f p is bound to the energies at f m and f n .
In contrast to the variance spectrum, the bispectrum is complex. The normalized magnitude and phase of the bispectrum are the bicoherence b 2 and the biphase β, given by (7) in which and denote the real and imaginary parts, respectively. According to Kim and Powers [59], the bicoherence characterizes the relative degree of coupling between three waves at f m , f n , and f p , which can be used to determine the bound variance at f p . Different equations to calculate the bicoherence have been presented in literature, all slightly differing in the way the bispectrum is normalized. Here, the equation of Haubrich [60] is presented, as was later also applied by Herbers et al. [52]. Furthermore, Elgar and Guza [61] showed that the statistical reliability of the bicoherence is insensitive to the normalization method. The intensity of the imaginary part of the bispectrum is indicative for the strength of the nonlinear energy transfers, which result in temporal or spatial changes in the spectrum (see, e.g., in [39,62]). The real and imaginary parts of the bispectrum are also closely related to the wave shape, as will be described in Section 2.3.
Every triad interaction appears in the bispectrum multiple times. Due to symmetry in the bispectrum, it is redundant to calculate and analyze the full bispectrum, but all sum and difference interactions are present in the triangle in ( f m , f n )-space bounded by ( f m = 0, f n = 0), ( f m = f N/2 , f n = f N/2 ), and ( f m = f N , f n = 0). For a detailed description of the symmetry regions in the bispectrum the reader is referred to Kim and Powers [59].

Wave Shape
The nonlinear shape of a wave can be described by its skewness (asymmetry w.r.t. the vertical axis) and asymmetry (asymmetry w.r.t. the horizontal axis). Skewness Sk and asymmetry As are third-order statistics [36], which are proportional to the real and imaginary parts of the bispectrum, normalized by the variance. The sea-swell Sk and As are computed as The factors 6 and 12 arise from the fact that the previously introduced bispectrum triangle covers 1/6 of the bispectral diagonals and 1/12 of the remaining bispectral area. The factor 2 in the denominator arises from the fact that the double-sided variance spectrum (Equation (2)) is only summed over the positive frequency range in order to determine the variance. Sk and As can be combined in the nonlinear wave shape parameter S [28]: S = Sk 2 + As 2 (10)

Bound Variance
This section outlines the method to obtain both the predicted as well as the observed bound wave variances, and the equivalent bound wave heights, from measurements.

Predicted Bound Variance for Equilibrium Conditions
Using second-order finite depth theory, the variance associated with the bound super harmonics can be predicted [38]. Based on this theory, the bound variance at a given frequency f p resulting from all sum interactions between primary sea-swell components which contribute to the variance at f p is calculated as in which r and s are discrete indices defining the directional bins such that θ r = r∆θ and θ s = s∆θ. D is the nonlinear coupling coefficient for seafloor pressure given by Herbers et al. [52] and d is the mean water depth. The corresponding bound super harmonic wave height is computed as As only sum interactions between primary sea-swell frequencies (of indices ≥ m min ) are considered in the calculation of P b,pred , the index of the first frequency where bound variance is present is 2m min .

Observed Bound Variance from the Bispectrum
As was mentioned in Section 2.2, Kim and Powers [59] 2 pointed out that the proportion of bound variance is related to the bicoherence. Integrating the bicoherence over all triad sum interactions contributing to frequency f p gives the proportion of variance at f p which is bound: The observed bound variance in the super harmonics can be written as 2 KP79 refers to the methods and equations from Kim and Powers [59].
The observed bound wave height is subsequently computed as Herbers et al. [52] 3 proposed an expression for the observed bound proportion of variance in the infragravity wave range due to difference interactions between primary sea-swell waves using the bispectrum. In a similar way, we define the observed bound proportion of variance at frequency f p associated with the sum interactions (i.e., bound super harmonics) as (16) in which α( f p ) is a weighting factor to account for differences in interaction strength between all triads contributing to frequency f p . As discussed by Herbers et al. [52], however, this effect is small for sea-swell waves and α can be assumed to be 1. Subsequently, the bound wave variance is expressed as Finally, the resulting bound wave height is computed as Comparing Equations (14) and (17) illustrates the differences between the two methods to compute the bound wave variance. While Kim and Powers [59] first calculate all individual bicoherences before summing them, Herbers et al. [52] first sum over the bispectrum and the cross products of the spectrum, subsequently square the absolute value of the summed bispectrum, and finally calculate the ratio. In Appendix A, the performance of both methods is investigated as a function of the statistical reliability of the expected values of the spectrum and bispectrum (Equations (2) and (5)). From this, it is decided to use the HEG94 formulation to compute the observed bound wave height, because it provides the most reliable estimate for a low number of degrees of freedom. This is of key importance in order to be applicable to field data in which the time series duration is usually limited in order to satisfy stationary conditions. Essentially, by first summing the bispectrum and spectral cross terms individually, additional averaging is applied which results in a statistically more reliable spectrum for the same duration. It should be noted that this additional way of averaging is only valid if the part of the bispectrum over which is summed is dominated by positive sum interactions, because otherwise the contributions of the sum and difference interactions cancel each other out. Thus, in the following H b,obs refers to the observed bound wave height calculated following the HEG94 method (Equations (17) and (18)). 3 HEG94 refers to the methods and equations from Herbers et al. [52].

CoastalGenesis2/SEAWAD Field Campaign
The 6-week CoastalGenesis2/SEAWAD field campaign was conducted in September and October 2017 by a consortium of universities and research institutes in order to improve the understanding of physical processes at the Ameland inlet [63,64]. The Ameland Inlet is a tidal inlet between the barrier islands Terschelling and Ameland at the north of the Netherlands. The barrier islands are located between the North Sea and the Wadden Sea. Seaward of the Ameland inlet, an ebb-tidal shoal has formed due to the deceleration of the ebb-tidal flow. During storm conditions with waves incoming from the north, it is on this shoal that the waves start feeling the bottom, reshape, and eventually break.
This paper focuses on pressure and velocity measurements obtained at nine locations clustered together on the seaward side of the ebb-tidal shoal (see also [32]). At two measurement frames (F4 and F5, see red dots in Figure 1), collocated pressure and near-bed velocity measurements were obtained. This was done using high-resolution, downward looking, Acoustical Doppler Velocity Profilers (ADCP) that measured the velocity profile over the bottom 50 cm of the water column and concurrent pressure with a sampling frequency of 4 Hz. At seven other locations surrounding the two frames, standalone pressure sensors (P1-P8, black dots in Figure 1) were deployed measuring the pressure continuously with a sampling frequency of 10 Hz. The two frames and three pressure sensors were aligned along a main transect, while the four other pressure sensors were deployed on side transects to investigate two-dimensional spatial variability. The mean water depth and the sensor height above the bed is given for all locations in Table 1.  Table 1. Overview of measurement locations, the mean depth and sensor height above the bed as deployed in the field, and the quantities which were measured (p stands for pressure, and u and v are the horizontal components of the velocity field).

Data Processing
Measured pressure is expressed in meters of water column by dividing by ρg, in which ρ is the density of sea water (=1025 kg/m 3 ) and g the gravitational acceleration (=9.81 m/s 2 ). For two reasons it is chosen not to reconstruct the surface elevation accounting for wave-induced pressure attenuation with depth: First, wave shape-induced sediment transport is driven by the near-bed wave shape. Second, reconstructing the surface elevation from the near-bed pressure requires the use of a transfer function that generally relies on linear wave theory. This introduces uncertainty, particularly in the high frequency range, where the bound higher harmonics are found. As a result, wave heights referred to in this paper are actually pressure-head derived wave heights, meaning that the actual wave height as would be observed at the surface is slightly larger.
Hourly time series are subdivided in 71 semi-overlapping blocks of 100 s. Subsequently, the spectrum, bispectrum, and bicoherence are estimated with 142 degrees of freedom (Equations (2), (5) and (6)). The wave height (Equation (4)) and peak period T peak = 1/ f peak are computed from the spectrum in which f peak is the frequency at which the variance is maximum. The observed bound wave height is computed following the method of HEG94 (Equations (17) and (18)). The focus is on the bound variance present in the superharmonic frequency range that originates from sum interactions between frequencies in the sea-swell range. Therefore, interactions involving components in the infragravity range should be excluded. This is done by using a frequency cut-off separating sea-swell and infragravity wave frequencies, here defined as f peak /2 following Roelvink et al. [49]. Hence, the index m min is defined such that f m min = f peak /2.
At F4 and F5, where collocated pressure and velocity signals are present, the 2D frequency-directional spectrum is computed using the Maximum Entropy Method [65] with ∆θ = 5 • , from which the energy-weighted mean direction θ and directional spreading σ θ are computed. In this paper, directions are presented in a Cartesian convention, thus the direction in which the wave is propagating is measured counterclockwise from the east. At the other locations, the normalized directional distribution of the closest frame is used in order to construct the 2D spectrum from the measured pressure spectrum using Equation (3). Subsequently, at all locations the predicted bound wave height is obtained using Equations (11) and (12). In order to make a fair comparison between the predicted and observed bound wave heights, they need to be calculated over the same frequency range, such that the same triad interactions are included in both estimates. This is done by using the same index m min as described in the previous paragraph. The uncertainty introduced by using a nearby directional distribution Θ on the calculation of H b,pred is discussed in Section 6.1.1.
In a few rare occasions, some minimal variance (at most 10% of the total) is present in directional bins opposing the peak incoming wave direction. These could not be explained by the concurrent wind conditions nor crossing sea states, and are most likely an artifact of the method used to construct the frequency-directional spectra. As the nonlinear interaction coefficient can be orders of magnitude stronger for opposing wave components [66], this minimal amount of variance can adversely affect the bound wave height prediction. Furthermore, these opposing components do not contribute to the sea-swell wave shape of interest because their bound wave length is much longer than the primary waves. Therefore, contributions to the predicted bound wave height for interactions with D < 0 are not taken into account.
The tidal current is a ubiquitous feature on an ebb-tidal delta. Its presence causes a shift between the absolute frequency ω and relative frequency σ, the Doppler shift. Formally wave theories, as used by Hasselmann [38], are valid in a moving frame of reference [67,68], thus using the relative frequency. This requires, however, that the current magnitude in direction of wave propagation is known, which is not the case at the seven standalone pressure sensors. Estimating the current in the direction of wave propagation U n at these locations would require a number of assumptions on wave refraction and the evolution of the current over the complex bathymetry. Therefore, it is decided to use ω at all sensors to compute k using the linear dispersion relationship. Subsequently, ω and k are used to determine D. The mismatch in frequency and wavenumber at the frames where velocities are measured is on average 3% and at most 10%. This mismatch is expected to be larger at the shallower sensors where k and U n are both expected to increase. The uncertainty on H b,pred introduced by this mismatch is discussed in Section 6.1.2.
Wave shape parameters (Sk, As, and S) are calculated from the bispectrum (using Equations (8)-(10)). In order to obtain the wave shape associated with bound superharmonics and to be consistent with the bound wave height formulations, only interactions in the bispectrum with f > f peak /2 are included.

Data Selection and Overview
For this study, it is chosen to only present cases in which the wave height at P1 exceeds 0.5 m. Cases with a lower wave height showed a negligible amount of bound variance as well as near-zero Sk and As, and are therefore disregarded in this paper. After this data selection, a total of 347 one-hour cases are included in this study. The majority of these cases (287) occurred during two storm events that coincided with the field campaign. To give an idea of the conditions, the wave height, period, direction, and mean water level at F4 during these two storms are shown in Figure 2. Another 60 cases outside of these storms are included, which occurred during four smaller events with 0.5 < H m0 < 1 m. For analyses involving all sensors, the total number of data points is 3123.

Spatial Evolution of the Spectrum and Bispectrum
This section discusses the spatial evolution of the spectrum, bispectrum, and associated wave statistics. First, results are shown for a single burst on 4th October, 12:00 ( Figure 3) representing storm conditions with a H m0 of 1.9 m, a T p of 10 s, and a mean direction θ of −72 • around slack tide. Subsequently, all bursts are combined in order to analyze general trends and the variability from these trends for all locations (Figure 4). Figure 3 shows the wave spectra (panels (a-d)) as well as the real (panels (e-h)), and imaginary (panels (i-l)) parts of the bispectra at four locations along the main transect (see Figure 3p for the position of the selected sensors). As expected in these energetic conditions, waves break when they reach the outer slope of the ebb-tidal shoal resulting in significantly lower variance levels at the shallowest sensors P2 and P4 (Figure 3c,d) than at F4 and F5 (Figure 3a,b). In addition to the primary sea-swell peak at f peak = 0.1 Hz, two secondary peaks can be clearly identified at the deepest selected locations, one in the infragravity range ( f < f peak /2), reaching maximum variance below 0.02 Hz, and another one around 0.2 Hz, i.e., at 2 f peak . The large (absolute) values of the real and imaginary part of the bispectra around ( f m , f n ) = (0.1 Hz,0.02 Hz) and (0.1 Hz,0.1 Hz) indicate that variance at these secondary peaks is, at least partly, nonlinearly coupled to the sea-swell primary peak. In the following, we focus on interactions between sea-swell frequencies as they are responsible for the buildup of the bound super harmonic variance. This means that interactions involving infragravity waves are excluded from further analysis. The importance of properly separating infragravity and sea-swell variance can be understood by looking at the bispectrum at for instance F4. Including interactions containing one frequency in the infragravity wave range and one in the sea-swell range (blue part in Figure 3e and red part in Figure 3i) would lead in this case to a lower (higher) integrated value for {B} ( {B}), modifying not only skewness and asymmetry but also the bound wave height estimate (see Equations (8), (9), (17) and (18)).
The magnitude of the real and imaginary parts of the bispectrum varies significantly along the transect (see changes in color scale in Figure 3e-l). At F4 and F5, the real part of the bispectrum is positive over the entire sea-swell range and of considerably larger magnitude than the imaginary part. This suggests that sea-swell waves are skewed, but not asymmetric, which is common for waves in deeper water that are not close to the breaking limit. At P2, the real and imaginary parts of the bispectrum are of the same order of magnitude, with {B} positive over the full sea-swell frequency range and {B} mostly negative. This means that waves are both skewed and asymmetric (saw-tooth shaped) at this location. Moreover, the consistently negative value of {B} around ( f peak , f peak ) at F4, F5, and P2 indicates that variance is being transferred from f peak to the sum frequency 2 f peak along most of the transect 4 , including at P2 where waves are already breaking. This nonlinear energy transfer contributes to the observed growth in variance of the first higher harmonic (2 f peak ) from F4 to P2, although other processes, such as (linear) shoaling and the changes in wave-induced pressure attenuation with depth also play a role. Finally, at P4, i.e., landward of the shallowest point of the ebb-tidal shoal, the magnitude of both {B} and {B} has significantly decreased while the total variance stays close to the one observed at P2. This suggests a weaker nonlinear coupling, and thus a decrease in bound wave variance, as well as more linear wave shapes.
These trends are confirmed when looking at the evolution of the integrated wave statistics for all nine sensors (Figure 3m-o). Initially, the total sea-swell wave height H m0 very gradually increases from P1 to F5 before significantly decreasing from F5 to P2 (Figure 3m). The bound superharmonic wave height, H b,obs , increases at a much higher rate than H m0 as waves propagate over decreasing water depth (Figure 3n), which is consistent with the variance increase observed at 2 f peak in Figure 3a-c. Interestingly, H b,obs keeps increasing beyond F5, while the total sea-swell wave height H m0 is already decreasing due to breaking. H b,obs finally decreases from P2 to P4 while H m0 stays almost constant. This suggests a release of bound higher harmonics over the shoal, as was observed under laboratory conditions by Beji and Battjes [69].
Finally, the evolution of the dimensionless wave shape parameters is visualized in Figure 3o. The skewness (Sk) gradually increases while moving into shallower water towards a maximum at P2 before it starts decreasing. In contrast, the asymmetry (As) is near-zero for most locations and is only of significance at P7 and P2, the two locations where the wave height is significantly decreasing and where the portion of breaking waves is expected to be the largest. As a consequence, the combined nonlinearity parameter S is close to the Sk except for P7 and P2, where the contribution of the As makes S slightly higher than Sk. 4 Different sign conventions can be found in the literature for the imaginary part of the bispectrum. In the present paper, we adopt the same representation as Norheim et al. [33], in which a negative value of {B( f 1 , f 2 )} is indicative of an energy transfer from f 1 and f 2 to f 1 + f 2 . Note that the opposite convention is used in, e.g., Herbers and Burton [39] and De Bakker et al. [51]. and dimensionless sea-swell wave shape parameters Sk, As, and S (o) at all locations (see panel (p) for information on the deployment depths). The dashed lines in the spectra and bispectra indicate f peak , the dotted lines f peak /2 and 2 f peak . The thick black diagonal lines in the bispectra are the symmetry lines. Note that the limits of the color scales for the bispectral plots are not all the same. The bathymetry of the main transects is shown in panel p, with ± one standard deviation indicated by black dashed lines representing bathymetric variability for cases with non-oblique incoming wave directions (−45 • < θ < 45 • ). Sensor locations are indicated by the red (if sensor is on main transect) and white boxes (if sensor is not on main transect; sensor is placed such that the mean depth is the same as the depth on the main transect).
The data analysis as presented above for the single case on 4th October is performed for all cases and the wave statistics and their variability is shown in Figure 4. The data is divided in two groups based on the significant wave height at the most offshore location (P1). Group 1 (black lines in Figure 4) contains all cases with 1 ≤ H m0,P1 ≤ 2 m and Group 2 (red lines in Figure 4) contains all cases with H m0,P1 ≥ 2 m. As expected, larger offshore waves (Figure 4a) lead to larger bound waves (Figure 4b) and more nonlinear wave shapes (Figure 4d). This figure also confirms the difference in the spatial evolution of H m0 and H b,obs observed already for the single data burst presented above. It shows that for the energetic cases (Group 2), the maximum value of H b,obs is systematically found in the area where H m0 is already decaying due to breaking.

Predicted and Observed Bound Wave Height
The evolution of H b,pred across the measurement transect is displayed in Figure 3n for the earlier selected burst and in Figure 4c for the entire dataset. From these figures, it can be seen that H b,pred increases at a similar rate as H b,obs up to sensor P7, where they both reach their maximum values but that the mismatch between predicted and observed bound wave heights increases when the depth decreases. After P7, the trends exhibited by H b,obs and H b,pred differ more significantly, with in particular a much stronger decrease in H b,pred than in H b,obs when waves propagate from P7 to P2.
Overall, these first comparisons suggest that the ability of Hasselmann's [38] theory to predict the bound wave height, as implemented in Section 3.1, varies spatially. To examine this in more detail, the data is clustered in three regions: the shelf (P1, P8, F4, P3, and P5), the seaward slope (F5 and P7), and the ebb-tidal shoal (P2 and P4). Figure 5 shows the observed bound wave height as a function of the predicted bound wave height for these different regions.
At the shelf (Figure 5a), the predicted bound wave height is very similar to the observed bound wave height with a strong correlation coefficient (R 2 = 0.94) and a linear regression slope of 0.96 suggesting that the observed bound wave response is in equilibrium with the local sea-swell forcing. At the steeper seaward slope (Figure 5b), the correlation between the predicted and observed bound wave height is still high (R 2 = 0.93), but the slope of 1.11 reveals a slight overestimation of the predicted bound wave height. This over prediction is consistent with nonequilibrium conditions where at the steeper part of the slope waves experience rapid changes in depth inhibiting the higher harmonics to fully develop.
On top of the shoal and right behind it ( Figure 5c) the predicted bound wave height deviates significantly from the observed bound wave height (lower correlation coefficient R 2 = 0.72 and a slope of 0.80). Although this linear regression slope indicates an underestimation of the predicted bound wave height on average, Figure 5c shows that the predicted bound wave height is both underand overestimated, depending on the conditions. At these relatively shallow locations the changes in sea-swell conditions are controlled by wave breaking, rapid refraction and wave current interaction and as such the equilibrium theory of Hasselmann [38] is not expected to hold. The errors introduced by currents, wave breaking, and directionality, and their effects on the predicted bound wave height are further discussed in Section 6.1.

Wave Shape as a Function of Observed Bound Wave Height
The relationship between the wave shape and bound wave height is examined next. Figure 6 shows the dimensionless combined wave shape parameter S as a function of the dimensionless predicted and observed bound wave height over sea-swell wave height ratios H b,obs /H m0 and H b,pred /H m0 , respectively. There is a very strong correlation, R 2 = 0.99, between the wave shape and the observed bound to total wave height ratio (Figure 6a). This strong relation between S and H b,obs /H m0 was expected, as S and H b,obs /H m0 are both computed by summing over the bispectrum and subsequently normalizing by the variance to some power. More specifically, S is obtained after dividing by the variance to the power 3/2 (see Equations (8) and (9)), while the bound wave height ratio is obtained after division by the significant wave height, i.e., the variance to the power 1/2. The strong relation between S and H b /H m0 suggests a mathematical connection between the two variables that still needs to be established. The small scatter around the best fit line can be partly explained by the reliability of the estimated spectrum and bispectrum. Additional tests (not shown) have indeed revealed that the scatter decreases for increasing number of degrees of freedom (DOFs). Although this relationship between dimensionless wave shape and bound wave height may seem trivial, to the authors knowledge, it has not been presented before.
The correlation between wave shape and bound wave height deteriorates significantly using the predicted bound wave height, H b,pred (right panel of Figure 6) with R 2 = 0.80. The observed scatter is related to the mismatch in predicted bound wave height with respect to the observed bound wave height (presented in Figure 5). Although scatter is present, the predicted value is reasonable for cases with a low proportion of predicted bound wave height H b,pred /H m0 < 0.15. For cases with a higher proportion of predicted bound wave height, the scatter increases, and thus the predictive skill decreases accordingly. It is for those cases where other commonly used wave shape parametrizations based on equilibrium conditions [25,26,28] also struggle to accurately predict the wave shape. Here, improvements in wave shape predictions can be readily obtained with better model predictions of the bound wave height ratio, as is evident from the comparison of the panels in Figure 6. Once the error in the predicted bound wave ratio is understood, it opens up the avenue for future modeling perspectives as discussed in the next section.

Errors in Determining the Predicted Bound Wave Height
Despite the fact that there is an acceptable agreement between the predicted and observed bound wave height in deeper water, a significant mismatch between the two is observed once it gets shallower (see Figure 5). In this latter part, the effects of directional spreading, ambient currents, and wave breaking on the prediction of the bound wave height are prevalent and discussed next.

Directional Spreading
A source of uncertainty is introduced by the assumptions required to obtain the frequency directional spectrum at the locations where velocities were not measured. A marginal spatial difference in directional spreading is observed between F4 and F5 (root mean squared difference is 4 • ), from which it can safely be assumed that the directional distribution at the other deeper located pressure sensors (P1, P8, P3, and P5) will not be too different. At the shallower locations (P7, P2, and P4), however, significant changes in depth and currents could lead to both an increase or decrease of the directional spreading. Easy ways to compute the refraction commonly rely on parallel depth contours are not applicable in this study due to the complex bathymetry. Therefore, this section discusses the effect a larger or smaller directional spreading has on the subsequent computation of H b,pred (Equation (12)).
A study by Herbers et al. [70] observed a maximum difference in directional spreading of 10 • along a cross-shore transect in the nearshore zone. This study was conducted in shallower and more nonlinear environment than our field campaign. Therefore, it is assumed that a mismatch of 10 • in directional spreading is an upper limit. The normalized directional distribution Θ is adapted such that the observed directional spreading is varied by plus or minus 10 • . This is achieved by taking the observed directional distribution to the power p and subsequently normalize it again to ensure that the sum of Θ is one for each frequency bin and hence the frequency distribution of the variance is not being affected: p > 1 gives more weight to the energetic directional bins, making the directional distribution narrower and thus decreasing the directional spreading. Conversely, p < 1 results in an increased directional spreading. For each burst, p is obtained using an optimization routine to obtain the desired 10 degree increase or decrease in the observed directional spreading. Figure 7 shows that an increased or decreased directional spreading of 10 • leads to an underestimation of 8% and an overestimation of 9% of H b,pred , respectively. Because a difference in directional spreading of 10 degrees is expected to be an upper limit in this location, the error introduced by using a wrong directional distribution is less than 9%.

Current
Tidal inlets are characterized by strong ebb and flood currents. It is known that these currents affect the wave dynamics due to wave current interaction and lead to a Doppler shift. This Doppler shift describes the difference between the absolute frequency ω (as observed at one fixed location) and the relative frequency σ (as observed in a frame of reference moving with the current): (20) in which U n is the current magnitude in direction of wave propagation. In absence of an ambient current σ = ω, whereas in presence of following (U n > 0) and opposing (U n < 0) currents σ < ω and σ > ω, respectively. A given observed absolute frequency thus results in a smaller wavenumber k for following than for opposing currents. The nonlinear interaction coefficient D, as used in Equation (12), is higher for longer waves (lower σ and k) than for shorter waves. However, as no reliable estimates of the current magnitude and direction were available at most measurement locations, ω and the corresponding k were used to obtain H b,pred without taking into account the tidal current. This could lead to an over or underprediction of H b,pred depending on the current direction.
The influence of the current on the over and under prediction of H b,pred is visualized in Figure 8, where the color of the data points indicates following currents (red dots: U n > 0.1 m/s), no currents (gray dots: −0.1 m/s < U n < 0.1 m/s), and opposing currents (blue dots: U n < −0.1 m/s). The current direction, magnitude, and wave direction at F4 are used to distinguish the different current conditions because the most reliable directional estimates were obtained at that location.
At the deeper sensor F4, no clear correlation is seen between the over and under prediction of H b,pred and U n (Figure 8a). However, the shallower it gets, the more evident the effect of the current on the over-and underprediction of H b,pred is (Figure 8b-d). The reason why the influence of the ambient current is more significant in shallower water, is related to the relative importance of the Doppler shift. This importance is described by the ratio kU n /σ or as U n /c in which c is the wave celerity. In deeper water, the current magnitude is lower and the celerity is higher, resulting in a marginal influence of the current. However, in shallower water, the current magnitude increases whereas the celerity decreases, making the importance of including the current more and more important.
To separate the possible effect of wave breaking from current effects in shallower water, sensor P4, located in a deshoaling zone and as such only marginally affected by wave breaking, is examined in more detail. At P4 (Figure 8d), for cases with a following current, on average the ratio of H b,pred /H b,obs is 0.78, whereas for cases with an opposing current this ratio is 1.38. Therefore, a significant part of the scatter observed in Figure 5c is related to the under-and overestimation of H b,pred for opposing and following conditions, respectively. H b,pred is under predicted in following current conditions, as the waves are actually longer than measured at the fixed measurement location. These longer waves should have a higher D and thus a higher predicted bound wave height. Vice versa, in opposing current conditions, the waves are shorter than measured, so H b,pred is overpredicted.

Wave Breaking
In the following, (near-)breaking conditions are defined as bursts for which As < −0.2. Using this criterion most wave breaking is observed at P7 and P2. To eliminate the previously discussed effect of the current at these two locations, a subset of 175 bursts is considered where the current in direction of wave propagation is negligible (|U n | < 0.1m/s), of which 53 and 50 meet the breaking criteria at P7 and P2, respectively. At P7, H b,pred is higher than H b,obs for 43 out of 53 breaking cases with an average over estimation of 13% (see Figure 9a). In contrast, at P2 H b,pred is lower than H b,obs for 47 out of 50 cases with an average under estimation of 23% (see Figure 9b). Figure 9. H b,pred as a function of H b,obs at P7 (a) and P2 (b) for the subset of data with small currents (|U N | < 0.1 m/s) with black markers for (near-)breaking conditions (As < −0.2) and gray markers for non-breaking conditions. A probable cause for the overestimation of the breaking waves at P7 was already discussed in Section 5.2, where it was explained that the wave shape can not instantly change in the rapidly changing bathymetry. The underestimation at P2 also seems to be related to the rapidly changing conditions. Due to wave breaking, significant energy loss is seen in the primary wave components. As H b,pred is proportional to the primary wave energy (Equation (12)), it decays at approximately the same rate as H m0 . It is known, however, that breaking conditions coincide with negative asymmetry, thus a negative imaginary part of the bispectrum, and as a result nonlinear energy transfers towards the bound higher harmonics. Therefore, although H m0 is decaying due to wave breaking, this positive energy transfer can lead to a growth of H b,obs . Additional support for this explanation can be found in Figure 4. Here, for the highly energetic cases (red lines) only a minor decrease in H b,obs is observed from P7 towards P2, whereas the decay in H m0 and H b,pred is much stronger. The fact that H b,obs decreases does not conflict with the earlier mentioned positive nonlinear energy transfer because the spatial change in variance is also affected by wave breaking, shoaling, and depth attenuation of the pressure signal.

Overall Validity of Equilibrium Bound Wave Height Theory
In the previous sections it is shown that a small amount of error in H b,pred is introduced by the unknown directional distribution at the shallow sensor but that more error is introduced when strong currents are present or when conditions are (near-)breaking. However, even when excluding cases with a strong current or (near-)breaking conditions, still a significant under or over estimation of H b,pred is observed. This is related to the validity region of the wave theory used by Hasselmann [38].
Many different definitions can be found in literature for the applicability region of the second order wave theory. Most of these describe the validity as a function of the Ursell number (see, e.g., in [71][72][73]), a dimensionless nonlinearity parameter. Here, we follow the definition as provided by Ruessink et al. [28]: Le Méhauté [74] showed that second order theory is only valid for Ur < 0.1. Madsen [75], Guza and Thornton [76] found a different value of Ur = 0.25, based on the argument that second order theory performs well as long as the secondary amplitude is at least 4 times smaller than the primary amplitude. According to Hedges [68], however, second-order Stokes theory can be used until Ur = 0.38. To investigate the consistency between the second order wave theory validity regions and the findings from our study, the RMSE between H b,pred /H m0 and H b,obs /H m0 is calculated for cases with low nonlinearity (Ur < 0.25) and high nonlinearity (Ur > 0.25). The threshold of 0.25 is used as it is an average of values found in literature. It is chosen to present the RMSE between the dimensionless wave height as this is used as a predictor for the wave shape. Cases with ambient currents and (near-)breaking conditions are excluded for the calculation of the RMSE. The RMSE between H b,pred /Hs and H b,obs /Hs is only 0.012 for cases with Ur < 0.25 whereas it is 0.050 for cases with Ur > 0.25, showing that the validity of second order wave theory also restricts the predictive skills of the bound wave shape.
In conclusion, it can be said that the wave shape can be accurately predicted using second order wave theory as long as waves are not breaking and nonlinearity is not too high. Furthermore, a proper prediction of the bound wave height requires directional information and, in presence of strong currents as in this study, information about the current field to properly account for the effect of the Doppler shift on the wavenumber.

Applicability in Different Areas
The disadvantage of many wave shape parametrizations is that they are site-specific, which means that they are only applicable for certain locations or conditions. In order to test the applicability of the findings of this paper, two additional datasets are analyzed. Although the SEAWAD dataset covers a wide range of conditions, it is limited in two ways: First, the shallowest sensors were in a mean water depth of 4 m, which limits kd to a minimum value of 0.4. Second, the wave periods observed during the campaign were rather short due to the fetch-limited conditions which are typically found in the North Sea. The two datasets that are used to validate the relationship between bound wave height and wave shape are the COAST3D [77][78][79] and SandyDuck [80,81] campaigns.
The COAST3D data was collected during a 6-week field campaign that took place at the beach of Egmond aan Zee in the Netherlands. It is thus mostly dominated, as SEAWAD, by wind-sea waves (T p = 6-12 s) but includes data in shallower depths. The SandyDuck dataset includes on the other hand 5 weeks of measurement along the North Carolina coast, and is therefore characterized by longer and more regular swell waves (T p = 12-20 s).
For both datasets, all bursts corresponding to a mean water depth below 0.5 m are discarded to exclude data that could be intermittently dry. For COAST3D, this selection results in a total of 5015 one-hour long pressure timeseries collected at nine locations in a mean depth ranging between 2.1 and 5.2 m. For SandyDuck, this means that a total of 1115-hour-long time series from seven pressure sensors are analyzed, covering a water depth range from 0.7 to 4.0 m.
The relationship between the dimensionless observed bound wave height and the dimensionless wave shape is presented in Figure 10 for the three campaigns. Note that observations at multiple locations are combined per campaign. For all campaigns the correlation coefficients are larger than 0.98, confirming the relationship between the bound wave height and wave shape. Furthermore, the slope of the linear regression line is fairly similar, between 2.7 and 2.8, indicating that the relationship between these two parameters is not specific to a certain area or conditions. This shows that the wave shape is known if it is possible to predict the bound wave height, independently of the area or conditions. The best linear fits between S and H b,obs (red lines in Figure 10) are presented for each field campaign. It can be seen that slightly more scatter is present for data points from the shallower field campaigns. The points that deviate the most from the relationship are shallow cases, in which significant infragravity variance was present in the spectrum (not shown).

Future Modeling Perspective
This paper shows that there is a direct relationship between the bound wave height and the wave shape regardless of location or type of conditions (see Figure 6). Therefore, with a properly predicted spatial and temporal evolution of the bound wave height the wave shape is known, which would in turn be instrumental for accurate calculations of the wave shape induced sediment transport [1,82,83]. Predictions of the equilibrium-bound wave height using second-order finite depth wave theory [38] proved to be accurate in relatively linear conditions (Ur 0.25). In more nonlinear conditions, and where significant changes in bathymetry and wave conditions are observed, there is a clear mismatch between the predicted and observed bound wave height. Although using higher order wave theories to predict the bound wave height might stretch the applicability region, it will not be able to capture non-local aspects that influence the bound wave height and associated wave shape. This omission can be overcome by using an evolution equation for the bound wave height taking into account bathymetric variability. In the following we discriminate between large scale phase-averaged models and more detailed phase-resolving models.
The effects of waves in large-scale morphodynamic models are commonly accounted for by phase-averaged spectral models (see, e.g., in [14,16]), which calculate the evolution of the variance spectrum in time and space. The energy transfers to the higher harmonics due to nonlinear triad interactions can be included using a source term function S nl3 . This source term has been derived from the bispectral evolution equations and is a function of the variance spectrum and a parametrized form of the biphase (see, e.g., in [26,84]). The S nl3 source term is extensively studied in recent years and has led to more and more reliable predictions of the nonlinear energy transfers in spectral wave models [40,41]. Therefore, essentially, the variance that is transferred to the higher harmonics (bound variance) is known and the bound wave height can be estimated by integrating this energy transfer in the down-wave direction. In this way, the evolution of the bound wave height is taken into account with the expectation that this will lead to a significant improvement of wave shape predictions in the shoaling zone, where the proportion of bound harmonics is consistently growing. When the bound higher harmonics are released and/or their variance decays due to wave breaking, this simple integrative approach is expected to yield less accurate results. Further improvements in estimating the wave shape can be achieved by modeling the effects of these processes on both the spectrum and bispectrum in order to be able to more accurately predict the bound wave height in all areas and conditions.
Alternatively, the evolution of the bound wave height can be obtained from more detailed phase-resolving models. These can be divided in time-domain and bispectral models. The former resolve the spatial evolution of the surface elevation, pressure, and velocities in time (see, e.g., in [85][86][87][88][89]). The latter resolve the spatial evolution of the spectrum and bispectrum (see, e.g., in [33,37,39]). For these type of models, the computation of the bound wave height as a wave shape predictor is unnecessary as the skewness and asymmetry can directly be obtained from the time series or bispectra that are provided by these models. Unfortunately, high computational times of such models prevent their usage to drive large-scale morphodynamic models. However, these detailed models can be used to study in detail the release of higher harmonics and (bi)spectral decay due to wave breaking and their effect on the evolution of the bound wave height. Subsequently, these effects can be included in the phase-averaged spectral models by modifying the S nl3 source term to improve bound wave height and wave shape predictions in complex systems where equilibrium conditions do not hold.

Conclusions
This paper shows that bispectral analysis of time series can be used to calculate the observed bound superharmonic wave height. In this study the method is applied to near-bed pressure time series, but it is also applicable for surface elevation or velocity time series. Despite several references in literature that such a methodology is not straightforward [36,[56][57][58], we found that summing over the bispectrum in a similar way as is done by Herbers et al. [52] for the bound subharmonic wave height, provides sufficient statistical reliability to obtain the bound superharmonic wave height from one hour time series. This does require that the methodology is restricted to the sea-swell frequencies and that this part of the bispectrum is dominated by positive sum interactions determining the bound wave height. In case difference interactions or negative sum interactions of crossing sea states contribute significantly to the bispectrum, the methodology should be treated with more care because the positive and negative contributions to the bispectrum might cancel each other out leading to an underestimation of the observed bound wave variance.
By using measurements at nine locations in the vicinity of an ebb-tidal shoal, a clear relationship (R 2 = 0.99) is found between the normalized observed bound wave height and the dimensionless sea-swell wave shape: S ≈ 2.75 H b,obs H m0 . As the same relationship is found for two other datasets that were collected along sandy beaches respectively dominated by wind-sea and sea-swell wave conditions, we conclude that it is insensitive to the environmental conditions. Thus, the wave shape is known at locations where we know the bound wave height. Knowing the wave shape at any given location would significantly increase morphological modeling capabilities because the wave shape induced sediment transport is resolved more accurately. However, as time series are not available in commonly used phase-averaged models, the bound wave height can not directly be computed. As an alternative, the bound wave height can be predicted using second-order wave theory by assuming equilibrium conditions [38]. From the analysis of the field data, it is concluded that the method of Hasselmann [38] is accurate in deeper water, but fails to accurately predict the bound wave height in (near-)breaking conditions or when nonlinearity is so high that second order wave theory is invalid. The accuracy of the predicted bound wave height is significantly improved if besides pressure also velocity measurements are known as the estimates are strongly dependent on the directional spread of the incoming sea-swell wave field. To improve future model capabilities, a next step is to add the evolution equation of the bound wave variance to spectral wave models. By including this evolution equation the bound wave height will be better predicted than using the method of Hasselmann [38] because it allows deviations from equilibrium conditions, which is key when considering wave transformation over rapidly changing bathymetry. The triad source term S nl3 , which is already included in such models, can be used as a source term for the bound wave variance. One of the challenges ahead is how to take into account the decrease of bound wave variance due to wave breaking and the release of bound superharmonics.
where f m = m∆ f with ∆ f = 0.01 Hz and N = 500. The super harmonics are then generated using the second-order theory of [38]: where D is the interaction coefficient introduced in Equation (11), calculated for a depth d = 5 m. The target bound wave height of the time series is computed as It should be noted that the bound wave height directly computed from the variance of the constructed time series (Equation (A2)) slightly deviates from the target bound wave height, because it is only a single realization of finite duration. However, when averaged over a sufficient amount of independently constructed time series, the ensemble-averaged bound wave height from those time series is the same as the target bound wave height.

Appendix A.2. Formulations and Degrees of Freedom
When working with observations, the number of degrees of freedom is typically increased by subdividing the measured timeseries in blocks before performing the spectral analysis or applying frequency merging [61,90]. Here, the use of synthetic data allows us to generate several realizations of our sea-state to estimate the expected values in the spectrum and bispectrum. The number of degrees of freedom is therefore increased by increasing the number of synthetic time series used in the calculations. The bound wave heights using Equations (15) and (18) for different DOFs are compared to the a priori known bound wave height from the time series in Figure A1.
A common way to characterize the statistical reliability of the bispectrum and bicoherence spectrum is to define the 95% bicoherence confidence interval, calculated as b 2 95% = 6/d.o. f .. If b 2 ( f m , f n ) > b 2 95% , B( f m , f n ) and b 2 ( f m , f n ) are considered statistically reliable (see, e.g., in [50,51,57,91]). To see what the effect of the bicoherence confidence interval is, Equations (15) and (18) are additionally applied with and without discarding bispectral estimates with b 2 < b 2 95% (solid vs dashed lines in Figure A1). It can be seen that the KP79 method (blue lines) needs many more DOFs to converge to the correct solution than the HEG94 method (red lines). This is because by first summing the bispectrum and energies individually instead of summing the bicoherences, additional averaging is added, leading to more reliable estimates for less degrees of freedom. If the bicoherence 95%-confidence interval is used to discard values, both KP79 and HEG94 underestimate H b . The underestimation is larger for few DOFs, because of the large b 2 95% -value that leads to the exclusion of many interactions where actual bound variance is present (all of them for the smallest values). Reasonable estimates of the bound wave height are only obtained with these methods for extremely large DOF values (>10 4 ).
The opposite behavior is seen for KP97 without a confidence interval. In that case the bound wave height is overestimated for all considered DOF values. This overestimation originates from the fact that, in this formulation, all interactions are considered, even the non-statistically significant ones, and that all of these interactions contribute positively to the bound wave height estimate as it depends on b 2 (Equation (13)).