Characterizing Wave Shape Evolution on an Ebb-Tidal Shoal

Field measurements of waves and currents were obtained at ten locations on an ebb-tidal shoal seaward of Ameland Inlet for a six-week period. These measurements were used to investigate the evolution of the near-bed velocity skewness and asymmetry, as these are important drivers for wave-induced sediment tranport. Wave shape parameters were compared to traditionally used parameterizations to quantify their performance in a dynamic area with waves and tidal currents coming in from different directions over a highly variable bathymetry. Spatially and temporally averaged, these parameterizations compared very well to observed wave shape. However, significant scatter was observed. The largest deviations from the parameterization were observed at the shallowest locations, where the contribution of wave-induced sediment transport was expected to be the largest. This paper shows that this scatter was caused by differences in wave-breaking, nonlinear energy transfer rate, and spatial gradients in tidal currents. Therefore, it is proposed to include the prior evolution of the wave before reaching a location in future parameterizations in numerical modeling instead of only using local parameters to predict wave shape.


Introduction
Tidal inlets and their deltas are highly dynamic systems [1,2]. Deep channels are formed by the in-and outgoing tidal currents. Sediment is exchanged through a network of channels between the tidal flats, deltas and adjacent barrier island coasts [3,4]. At the seaward end of the ebb channel, a relatively shallow area is located, the ebb-tidal delta. It is shaped as a result of the decelerating tidal flow and the influence of waves from offshore. The ebb-tidal delta serves an important role in these systems. During storm conditions, it provides shelter to the back-barrier basin by dissipating high incoming waves. Furthermore, the ecosystem has a high ecological value containing a large variety of benthic species. According to Finley [5], morphodynamic development of the ebb-tidal delta depends on the tidal flow and waves approaching from offshore. Traditionally, it was considered that waves stir up the sediment, which is then transported by the currents [6]. However, as stated by Nielsen [7], wave-induced sediment transport itself can be of significant influence, even in presence of strong currents. Properly accounting for wave-induced sediment transport in process-based numerical models (e.g., [8,9]) is of high importance to represent bar formation along a straight coast. In a more complex environment, such as a tidal inlet, the evolution of the ebb tidal shoal and delta is formed by the wave-induced sediment transport. A model study by Chen et al. [10] concluded that morphological changes at the ebb-tidal delta are dominated by wave-induced sediment transport, where the other surroundings are dominated by the current-induced sediment transport.
Nonlinearity of the near-bed velocity is an important parameter determining the wave-induced sediment transport [11]. Nonlinear wave shape can be defined as the skewness (vertical asymmetry) and asymmetry (horizontal asymmetry) of a time series [12]. The initial sinusoidal wave shape, as occurs in deep water far away from the coast, becomes nonlinear when moving into shallower water. First, it transforms in a skewed shape with short peaked crests and long flat troughs. Due to the faster propagation of the crest with respect to the trough, waves become pitched forward resulting in an asymmetrical wave shape. Eventually, as the waves become more and more asymmetric and break, skewness decreases. Analogous to gradients in depth, gradients in coastal currents result in wave shape changes. Wave skewness and asymmetry are intra-wave properties, which are not resolved in the wave-averaged models that are used to drive morphodynamic models. Therefore, it is essential to have an accurate parameterization to include their contribution to the total sediment transport. A recent study has shown that properly parameterizing the near-bed velocity substantially improves morphodynamic predictions on the decadal time-scale [13].
First sediment transport models including the effects of velocity skewness were developed by Stive [14], Roelvink and Stive [15]. The importance of adding velocity asymmetry was shown by Drake and Calantoni [16]. Doering and Bowen [17] provided a parameterization of orbital velocity skewness and asymmetry based on field measurements obtained from four beaches. They found a correlation between skewness/asymmetry and the Ursell number, namely Ur = 3 4 a k 2 d 3 , a local nonlinearity parameter based on wave amplitude a, wave number k, and mean water depth d. Although the relationship between bi-amplitude (skewness and asymmetry combined) and the Ursell number is observable in the shoaling and surf zone, they stated that it is less clear in the surf zone. In a similar study, based on many more data points (34.000 vs. 48), Ruessink et al. [18] found a similar relationship as a function of the Ursell number. Differences between the two parameterizations are only significant for Ur > 3. Rocha et al. [19] compared available nonlinearity parameterizations based on field experiments and state that the Ruessink et al. [18] parameterization can predict skewness and asymmetry on low-sloping beaches for specific wave conditions, but that it does not work as well for steeper beaches and longer wave periods. In a follow-up study, Rocha et al. [20] used a numerical model to conclude that a decreasing offshore wave steepness, decreasing bed slope, and decreasing spectral width all lead to a higher maximum nonlinear wave shape (skewness and asymmetry combined) and different development along the beach profile. Although not included in the previously mentioned parameterizations [17,18], the effects of wave steepness [21], spectral bandwidth [22] and bottom slope [22][23][24] on skewness and asymmetry have been shown in previous research.
Although several studies (e.g., [17,18,20]) have investigated the parameterization of the wave shape, most of them focused on (barred) beaches in the surf zone. It is questionable how applicable these parameterizations are at an ebb-tidal delta for a number of reasons. Firstly, at an ebb-tidal delta the bathymetry is much more complex than for an open coastline with milder alongshore variations, leading to a larger range of incoming directions and spreading. Secondly, tidal currents might play an important role. Tidal currents opposing the direction of wave propagation can lead to decreasing wave length, increasing wave height and eventually even wave breaking (e.g., [25,26]). Therefore, it is likely that tidal currents have an influence on the wave shape. Whereas in coastal waters, tidal currents are often (almost) perpendicular to the wave direction, hence barely affecting the waves, at an ebb-tidal delta, the currents and waves can face each other under any arbitrary angle. Thirdly, depths of the ebb-tidal delta are generally larger than those investigated in previous studies.
In this paper, we compare near-bed velocity skewness and asymmetry from a new dataset of measurements at an ebb-tidal delta to existing parameterizations. Data were obtained during the CoastalGenesis2/SEAWAD field campaign in September 2017 at the Ameland Inlet, The Netherlands [27]. In Section 2, the field campaign and data-processing techniques are discussed. Temporal and spatial variations in wave shape, their dependence on the Ursell number and the role of prior wave transformation is addressed in Section 3. Subsequently, in Section 4, we discuss the consequences of the results on sediment transport, and ways to improve future wave shape parameterizations. Concluding remarks are given in Section 5.

Field Campaign
To obtain the system knowledge necessary to plan future sand nourishments in ebb-tidal deltas, a large field campaign was conducted in September 2017 at the Ameland Inlet, the CoastalGenesis2/SEAWAD field campaign [27]. The Ameland Inlet is located in the north of The Netherlands, between barrier Islands Ameland and Terschelling. The area is highly dynamic, with combined action of waves and tidal currents. The main tidal channel has depth-averaged currents up to 2 m/s in a maximum depth of 20 m. North of the channel, a shallow sandy area is formed by the deceleration of the diverging ebb tide, the ebb-tidal delta (see Figure 1). This area is shaped by the combination of waves and tide. Waves, locally generated in the North Sea, approach the ebb delta from North to Northwestern directions. They are short (5-10 s) and very variable in height (0.5-5 m). Bathymetric map of Ameland ebb-tidal delta with instrument frames (black triangles), pressure sensors (pink circles), and offshore wave buoy (green circle). In the bottom left corner is a zoomed-in map of the ebb-tidal shoal. where most of the instruments were located.
Five instrument frames were deployed for a six-week period, measuring velocities, waves, water levels, sediment concentrations, bed forms, conductivity, temperature and density. Eight pressure sensors surrounded the frames to get more spatial information on waves. Frame and pressure sensor locations are indicated in Figure 1. Additionally, drifter measurements and 13-h tidal-inlet discharge measurements were performed. Furthermore, grab samples, box cores, and water samples were obtained to characterize the sediment composition in the bed and water column. For a detailed description of the field campaign, see van Prooijen et al. [27].
For this study, the focus was on near-bed, high-frequency measurements of velocity and pressure. At the frames, this was measured by a downward-looking High-Resolution Acoustic Doppler Current Profilers (ADCP), measuring the velocity profile in the lowest 0.5 m of the water column above the bed and pressure 0.5 m above the bed at a frequency of 4 Hz. The standalone pressure sensors (PS) only measured pressure at a frequency of 10 Hz approximately 0.5 m above the bed. To compute the depth-averaged currents, the above-mentioned downward looking ADCPs were used in combination with upward looking ADCPs measuring the velocity profile from 2.3 m above the bed to the free surface at a frequency of 1.25 Hz. Mean water depth, bed slope, instrumentation, and sampling frequency are given in Table 1. The bed slope (indicated in Table 1) was computed over a distance of 200 m in the mean direction of wave incidence (north-northwest), in order to reflect the slope the waves encountered before reaching a certain sensor. Figure 2 shows the main bathymetric transect including the positions of the measurement locations. The nine sensors available were grouped into four clusters based on their position on the ebb-tidal delta: • Shelf: PS1, PS8, F4, and PS3. These sensors were at the milder shelf in mean depths between 10.5 and 8 m. • Seaward slope: PS5, PS7, and F5. These sensors were at the steeper seaward edge of the shoal in mean depths between 8 and 5 m. • Shoal: PS2 was at the transition between the seaward slope and the flat top of the shoal. • Landward slope: PS4 was at the basin end slope of the shoal. It was 300 m more landward than PS2 in a region where the depth starts increasing after the shallowest point, which was located between PS4 and PS2.  This standard deviation was computed over 91 transects with the same length as the main transect rotated through a fixed point at PS2 between −45 • and +45 • from the main transect to visualize the three-dimensionality of the area. Sensor locations are indicated in a red circle if they were on the main transect and in a white circle if they were not on the main transect. The pressure sensors that were not on the main transect were placed such that they were at the depth at which they were located.
During the field campaign, a wide range of conditions were encountered with two storms alternated by calm periods. In Figure 3, the conditions throughout the campaign are presented. The offshore wave conditions (Figures 3a-c) were obtained from a wave-buoy (green circle in Figure 1), and tidal elevation and depth-averaged velocities (Figures 3d-f) from the most offshore measurement frame F4. Waves, coming in from the North Sea (northern to western directions), were between 0.5 and 5 m. Two storms occurred, with significant wave heights up to 5 m. Due to limited fetch, wave periods were generally short, 6-8 s, with peaks to 10 s during the second storm. Water levels varied with a tidal range of 2.6 m. Additionally, a storm surge of up to 1 m was observed. At the seaward slope of the ebb-tidal delta, the tidal velocities ranged from 0.5 m/s towards southeast during flood to 0.5 m/s northwest during ebb. During the first storm, an outflow exceeding 1 m/s was observed, probably related to wind-driven outflow from the Wadden Sea basin [28].

Data Processing
All instruments operated in burst mode, with bursts of 30 min. The way the signals were processed in time-averaged quantities is described in this subsection.

Velocity
High-resolution (4 Hz) profile measurements of velocity were obtained 0.5 m above the bed by the downward-looking ADCPs. Additionally, lower resolution (1.25 Hz) profile measurements were obtained above the frame, from 2.3 m above the bed to the water surface. All velocity signals were filtered based on correlation [29], despiked [30,31], and converted to an East-North-Up reference frame. Depth-averaged currents were computed from the upward-and downward looking ADCP profiles (with the gap between the two instruments linearly interpolated), and characterized by a mean current magnitude U and direction θ U .

Pressure
Water pressure p was obtained by subtracting the air-pressure p air from pressure signals p tot measured by the ADCPs (locations F4 and F5) and pressure sensors (PS1-PS8). Mean water depth was calculated using d = z p + p /ρg, in which z p is the height of the pressure sensor above the bed, p is the burst averaged pressure, ρ is the water density and g is the gravitational acceleration. The surface elevation η was reconstructed by using a frequency dependent correction factor, based on linear wave theory, accounting for pressure attenuation over depth. Significant wave height H m0 = 4 √ m 0 and period T m−1,0 = m −1 /m 0 were obtained by applying spectral analysis on the surface elevation in the short-wave frequency range (0.05-0.3 Hz), in which m n is the spectral moment of order n. Additionally, 2D spectral analysis was performed on surface elevation, eastward and northward velocity time series in the same frequency range to obtain the mean wave direction θ w and spreading σ θ w . Wave number k was computed by iteratively solving the linear dispersion relationship: with absolute frequency ω = 2π/T m−1,0 , relative frequency sigma, and U n the depth-averaged current in direction of wave propagation, which can be computed from U, θ U , and θ w . The estimation of depth-averaged currents, which are required to compute k for locations PS1-PS8, where velocity measurements were not available, is described in Section 2.3.

Calculation of Nonlinear Parameters
Wave skewness, asymmetry and the Ursell number were derived from time series of pressure and velocity. To only consider the short-wave contribution, a frequency filter of 0.05-0.3 Hz was applied to the time series. For the first analysis (Section 3.1), the dimensional skewness and asymmetry was used, because it is reflective of the amount of wave-induced sediment transport. They were computed as: where q is the filtered time series of a quantity, which can be either velocity or pressure, H{...} denotes the imaginary part of the Hilbert transform, and ... indicates time averaging over the burst. To compare our dataset to previous field campaigns and parameterizations, skewness and asymmetry were also presented in the conventionally used dimensionless form: As the measured Sk and As were compared to the Ur-based parameterization introduced by Ruessink et al. [18], the Ursell number was estimated exactly as in this paper, i.e., with wave amplitude a = H m0 /2, k calculated with Equation (1) using T m−1,0 as a period estimate, and d the mean depth.

Estimation of Near-Bed Velocity Time Series from Pressure Signals
Near-bed velocity data were only available at three locations, in mean water depths between 7 and 10 m. To extend our dataset to a larger depth range, between 4 and 12 m, near-bed orbital velocities were also estimated from the pressure time series measured by the standalone pressure sensors PS1-PS8. The near-bed orbital velocity signal was computed from the pressure signal using linear wave theory. Since the correction depends on σ and k, the correction factor was frequency-dependant and thus should be applied in the frequency domain. The resulting near-bed, orbital velocity time series are given by: where subscript i refers to the wave component of frequency f i . To verify this method, reconstructed velocity skewness and asymmetry were compared to directly measured velocity skewness and asymmetry at the measurement frames ( Figure 4). For F4, the method systematically over predicts skewness (red data points in Figure 4), which is also reflected in a Root-Mean-Squared-Error (RMSE) of 0.12. For F1 (RMSE = 0.06) and F5 (RMSE = 0.07), skewness is very well predicted. Not a lot of asymmetry was observed at the measurement frames, and a similar error (RMSE = 0.04) was found for all three frames. . Correlation between (a) skewness and (b) asymmetry directly measured from the near bed velocity (Sk vel and As v el) and from reconstructed near-bed velocity (Sk p2v and As p 2v) at measurement Frames 1 (blue), 4 (red), and 5 (yellow). The black dashed line is unity.

Depth-Averaged Currents
At locations PS1-PS8, the currents were not measured. They were needed, however, for solving the dispersion relationship (Equation (1)). Therefore, an estimate of burst averaged currents at these locations was provided by a Delft3D simulation of the campaign, including the effects of wind, waves and tide for the duration of the field campaign. The setup and details of the numerical model are described in Nederhoff et al. [32], and the validity in predicting the correct magnitude and direction of the currents is shown in Appendix A by comparing predicted and measured depth-averaged currents at F1, F4, and F5.
The importance of including time-averaged currents in the dispersion relationship is related to the magnitude of kU n compared to the magnitude of σ. The importance is thus higher for shorter waves and larger tidal currents in or against direction of wave propagation. Because of mass continuity, largest currents were observed at the shallowest location, PS2, which was therefore the location where the effect of the current on waves was expected to be most noticeable. At this location, when not including currents, the average and maximum error for k was 6% and 33%, respectively. Since k was also included in the conversion from pressure to surface elevation (Section 2.2.2) and near-bed orbital velocity (Section 2.2.4), the aforementioned error in k influenced σ, a, Ur, Sk, and As with maximum errors of 20%, 12%, 38%, 3%, and 3%, respectively. The error introduced by not including currents is visualized in Figure 5, which clearly shows that a, Sk, and As were barely influenced, whereas k, σ, and Ur changed significantly when including the currents.

Data Overview and Selection
Significant wave height H m0 , dimensional skewness Sk dim and asymmetry As dim at the deepest (PS1, d = 10.4 m) and shallowest (PS2, d = 4.3 m) locations, are presented in Figure 6. The magnitudes of Sk dim and As dim are highly variable in time (Figure 6b,c). The largest absolute values were mainly observed in stormy periods at 13 September and 3-7 October, where H m0 > 2 m. Depth-induced shoaling generally leads to slightly increasing wave heights between PS1 and PS2, apart from the two storms in which wave breaking leads to a decrease in wave height (Figure 6a). It can be seen that Sk dim increases while moving to the shallower sensor, and a clear tidal pattern can be seen in the Sk dim during the second storm. Large values of |As dim | occur when conditions are close to wave breaking or once the waves are breaking. This explains why near-zero As dim was observed at PS1, but that waves get asymmetric at PS2.
For the following analyses, the focus was on cases with considerable Sk dim and As dim for two reasons: first because this reflects when wave-induced sediment transport occurs, and second to prevent nonphysically high values of dimensionless skewness and asymmetry when u 2 3/2 is really small. A combined nonlinear parameter B dim = Sk 2 dim + As 2 dim , calculated at the most nonlinear location PS2, was used to determine which cases to include. The threshold was set such that 99% of cumulative nonlinearity B dim was taken into account in the following, which was caused by only 42% of the bursts. Eventually, 813 half-hourly bursts were selected (see shaded areas in Figure 6), i.e., 17 days of unique data in which significant Sk and As are present.

Dependance of Wave Shape Parameters on the Local Ursell Number
Traditionally, the Ursell number Ur is used to explain the variations in skewness Sk and asymmetry As [17,18]. For our dataset, Ur ranges from 0 to 0.7, i.e., values similar to those typically found in the shoaling and outer surf zone. In this range, both Sk and As are expected to increase for increasing Ur, with the increase in As starting at higher Ur-values (see parameterization by Ruessink et al. [18], dotted blue lines in Figure 7).
When Sk and As of all locations were combined and binned as a function of their local Ursell number, the results compare remarkably well with the parameterization of Ruessink et al. [18], as can be seen in Figure 7a,b. An almost perfect match was observed between the bin-meaned Sk of all locations and the parameterization. A similar match was found for the As bins at low Ur, but for Ur > 0.3 observed As is significantly lower than the parameterized value. Overall, this suggests that, although based on measurements on barred beaches, the parameterization is applicable in a more complex environment such as an ebb-tidal delta. However, it should be noted that the standard deviation per bin is significant with respect to the mean of the bin, which is also reflected by the extensive clouds of data points in Figure 7a,b. As shown in the following, although the wave shape parameters are well predicted by parameterizations depending on Ur for the entire dataset, when considering individual locations or specific wave conditions, the parameterization can systematically under-or overestimate Sk and As.
At the shelf, Sk is fully controlled by Ur (Figure 7c), with correlation coefficients R 2 between 0.87 and 0.89 (see Table 2). When it gets shallower, two trends can be observed. A larger scatter in Sk was observed at a given Ur, and the slope of the best linear fit between Sk and Ur decreases from the shelf to the seaward slope (Figure 7e). On top of the shoal and at the landward slope (Figure 7g,i), the correlation between Sk and Ur decreases to 0.28 and 0.13, thus variability in Sk cannot be well described by the variability in Ur.  Table 1): (a,b) all sensors combined; (c,d) the shelf; (e,f) seaward slope; (g,h) ebb-tidal shoal; and (i,j) landward slope. Individual data points, visualized by the dots, are in red when local wave breaking is expected (G > 0.5) and grey if no wave-breaking is expected (G < 0.5), based on Equation (8). Black error bars indicate mean per bin ± one standard deviation. Five equally sized bins are defined between the minimum and maximum observed Ur within each subset of the data. The blue dashed line is the parameterization by Ruessink et al. [18]. As is near-zero at the shelf (Figure 7d) although Sk is already present (Figure 7c), which is as expected for Ur < 0.3 (see blue lines in Figure 7). The near-zero As at this location explains the low correlation between As and Ur ( Table 2). At the seaward slope, As starts to develop indicating pitched forward waves, which are breaking or close to breaking (Figure 7f). On top of the shoal, high variability in As was observed (Figure 7h). There is no clear trend as high |As| occurs for a wide range in Ur between 0.25 and 0.4. Moreover, for this range, cases were also observed where As is near-zero although Ur is up to 0.45. For the bursts with non-zero As, there is a clear difference between the slope of As as a function of Ur at the seaward slope ( Figure 7f) and the shoal (Figure 7h), which leads to an increased scatter when all sensors are combined (Figure 7b). At the landward slope, where the depth has slightly increased As is near-zero for the full Ur-range. This behavior was also observed and modeled by Elgar et al. [33], who found that highly skewed and asymmetric waves can become symmetrical again in increasing water depths.
In summary, although it is slightly underestimating Sk at the shelf and overestimating As at the seaward slope, the parameterization by Ruessink et al. [18] is capable of predicting Sk and As reasonably well at the shelf and the seaward slope (see dashed blue lines in Figure 7c-f). On top of the shoal, the data cloud suggests that no parameterization based on Ur only can properly predict Sk and As. The same holds for Sk at the landwards slope (see dashed blue lines in Figure 7g-

i). No
As was observed behind the ebb-delta, although for higher Ur, some As would be expected based on the parameterization (see dashed blue lines in Figure 7j). The fact that observed Sk and As for a given Ur-band (e.g., 0.3-0.4) are very different depending on the location and that large scatter was observed at the shallow locations suggest that more physical processes need to be accounted for than local Ur only.

Role of Wave Transformation
Deviations with respect to the parameterization by Ruessink et al. [18] occur for high Sk in deeper water ( Figure 7c) and high |As|-values (Figure 7f,h). Besides that, a significant increase in variability was observed on top of the shoal (Figure 7g,h). All these differences are explored here by examining wave transformation effects associated with wave breaking, non-linear interactions and currents.

Variability Caused by Wave Breaking
High |As|-values typically correspond to breaking waves. Therefore, the high values at PS2 and PS7 (Figure 7f,h) suggest that differences in breaking patterns could describe part of the variability. The relative wave height, H/d, is a common parameter to discriminate between breaking and non-breaking waves. To quantify whether waves are expected to break between PS1 and PS2, a new parameter, G, was defined, which is the predicted wave height at PS2 divided by the water depth at PS2. The predicted wave height at PS2 was obtained by shoaling the wave height from PS1 towards PS2 assuming energy conservation and unidirectional wave propagation. In equation form, G, was computed as: G = H m0,1 c g,1 /c g,2 d 2 (8) with group celerity c g i at location i computed according to linear wave theory as: If G is higher than 0.5, which is within the common range of breaking criteria for irregular waves over sloping bathymetry for γ = H m0,br /d [34], it is assumed that wave breaking took place. The added value of using G over the more traditional local estimate of H/d is that it not only indicates whether waves are expected to break or not, but also for how long waves have been breaking and thus what proportion of a random wave field is expected to break. Thus, G essentially describes how far in the surf zone the data have been collected. Figure 8a,b shows how the wave shape parameters, Sk and As, vary as a function of G. For G > 0.5, the absolute values of Sk and As are also high. For low values of G, a wide range in Sk and near-zero As was observed. This difference is consistent with the fact that shoaling waves start developing Sk earlier than As. If we use a threshold of G = 0.5, we clearly see that the two distinct trends in As at PS2 are explained by the fact whether waves have broken or not before reaching PS2 (see red markers in Figure 7h). In addition, we see that most variability in Sk is associated with low G values, i.e., non-breaking cases. Figure 8c finally shows that Ur, although commonly used to predict Sk and As, is not able to discriminate between breaking and non-breaking waves (see for instance the wide range of G for Ur = 0.3 in Figure 8c). To further investigate the role of wave transformation on the observed variability in Sk and As, we examined the changes in wave characteristics over the bathymetric transect ( Figure 2). More specifically, we examined the changes in Sk, As and Ur during wave propagation from PS7 (d = 5.3 m) to PS2 (d = 4.3 m), as PS7 was the last sensor before scatter in Sk and As as a function of Ur significantly increased. Sk is remarkably constant between the two sensors regardless of the fact whether waves are breaking or not (Figure 9b) while |As| stays near-zero for non-breaking cases and slightly increases for breaking cases (Figure 9c). The increase in |As| for breaking cases is explained by the fact that depth has decreased, thus a larger proportion of the waves will be breaking waves at PS2. This is consistent with the strong dependence of As on G described above. In contrast to the marginal changes in Sk and As, Ur significantly changes between PS7 and PS2, with a clear decrease for breaking waves and an increase for non-breaking waves (Figure 9a). The observed decrease in Ur is caused by increasing k, and decreasing a due to breaking. The observed increase in Ur for non-breaking cases, on the other hand, is a result of an increase in a and k and a decrease in d. This decrease of Ur for breaking waves (high Sk) in contrast to the increase for non-breaking waves (low Sk) explains the increased scatter at PS2 compared to PS7 (black vs. grey and red markers in Figure 9d, in which the red and grey arrows represent the general trend for breaking and non-breaking cases, respectively). In a similar manner, scatter in asymmetry is increased from PS7 to PS2 (see Figure 9e). It also explains why the most asymmetric waves were observed at lower Ur numbers at PS2 than at PS7 (compare Figure 7f,h). Thus, apparently, when wave-breaking occurs, the Ur-number is not the best parameter to predict the wave shape, since large scatter was observed. An approach to reduce the scatter in parameterizations predicting the wave shape is discussed in Section 4.2. As for non-breaking (grey) and breaking (colored) cases, from PS7 to PS2. Ur, Sk and As at PS2 are visualized as a function of Ur, Sk and As at PS7, respectively. Colors for the breaking cases indicate U N (m/s). The black dashed line has slope 1, thus indicates no change. The effect on (d) Sk and (e) As as a function of Ur is shown by visualizing data-points at PS7 in grey (non-breaking) and red (breaking) and at PS2 in black (all). The general evolution of data clouds from PS7 to PS2 for non-breaking and breaking conditions is visualized by grey and red arrows, respectively.

Variability Caused by Nonlinear Energy Transfer Rate
Nonlinear interactions lead to energy transfer from the primary harmonics to the higher harmonics and are thus essential for the development of Sk and As. From Hasselmann [35], it is known that the wave length, and thus period, of the primary wave field is a key variable influencing the nonlinear energy transfer strength. As a result, longer waves are expected to develop Sk sooner and at a higher rate than shorter waves. This was confirmed by the modeling study of Rocha et al. [20], who showed that waves with higher period have faster development of Sk and As.
In line with findings by Hasselmann and Rocha et al. [20,35], in this dataset, the incoming wave period T m−1,0 (at the deepest sensor PS1) explains a significant part of the variability in Sk as a function of Ur at the shallow sensors. Figure 10a shows that, within each Ur-bin, Sk varies as a function of T m−1,0 at PS1. Although the wave period is included in Ur through k, a long low wave can have the same Ur as a shorter higher wave when in the same depth, whereas a higher amount of Sk and |As| is expected for the longer wave based on the different nonlinear energy transfer rate.
Besides variability at the shallow sensors, the higher Sk in deeper water is also related to high incoming wave periods. The longer waves, which feel the bottom earlier, have already transformed into a more skewed wave shape than predicted by the parameterization by Ruessink et al. [18] based on local Ur (Figure 10b).

Variability Caused by Currents
At PS2, a large range of As and Ur is seen for the breaking waves (red dots in Figure 7h), which cannot be explained by a difference in incoming wave period because all those waves are long (9-10 s, not shown). It can be seen in Figure 11c that opposing and following currents lead to two distinct branches. This results in a different As observed for opposing and following waves for the same Ur. The two different branches can be explained by the fact that, in breaking conditions, Ur decreases more when the waves propagate from PS7 to PS7 for opposing currents than for following currents (compare blue and red colors in Figure 11a). This stronger decrease is due to the shortening of wave length if the current magnitude increases in shallow water during the opposing currents, and subsequently the change in k influences Ur.

Relevance in Morphodynamic Modeling
Long-term simulations using process-based models do not model waves to such extent that Sk and As are resolved. Therefore, parameterizations are needed to predict Sk and As, usually as a function of local Ur.
A recent numerical study for this specific tidal inlet by Reniers et al. [36] investigated the relevance of wave-shape induced transport compared to the traditional transport where sediment is stirred up by waves and transported by the currents. They found that the contribution of wave-shape transport is the dominating mode of transport in stormy conditions at the shallow locations, and can counteract the tide-induced sediment transport leading to a net onshore sediment transport.
This dataset has shown that it is in those conditions, when wave-shape transport is dominant, that most deviations from the parameterizations are observed (e.g., for breaking waves and at the shallow sensors). Since wave-shape transport is always onshore directed, small errors will accumulate over long-term simulations and can thus significantly influence the development of the shallow ebb-tidal shoal.

Correcting for Delayed Response in Wave Shape
Changes in Sk and As when waves propagate from PS7 to PS2 are much less significant than changes in Ur (see Figure 9). This suggests a delayed response of the wave shape, and in particular of the skewness, to changes in the Ursell number when waves propagate towards shallower water depths. In the following, we want to check if using an Ursell number corresponding to an earlier stage of wave transformation, Ur * (defined at depth d * = d + ∆d, with ∆d > 0) has better predictive skills than using the local Ursell value Ur (at depth d). Since PS4 (at the landward slope) is in a different regime, where ∆d < 0, it is excluded for this analysis. We assume that, at the deepest sensor, the skewness is close to being fully developed and thus has an optimal fit with the Ursell number (see also the high R 2 value in Table 2 for PS1), which we take as a reference for the optimization procedure. More specifically, we iteratively estimate, for each sensor, the depth correction ∆d needed to get the best fit with the observation at PS1 in terms of both slope and intersect of the least square linear fit. The modified Ur * is computed as: with k * obtained by deshoaling k from depth d to d * using linear wave theory, and amplitude a * , at depth d * obtained by linear interpolation between the location under consideration and the location before it in deeper water. Table 3 shows per location the ∆d by which the Ursell number needs to be corrected in order to obtain the optimized Ur * . It is seen that the delay is negligible (∆d < 0.25 m) at the deeper sensors (PS1, PS8, F4, PS3, and PS5), but is significant at the shallower sensors on the steeper seaward slope (F5, PS7, and PS2). Figure 12 shows Sk and As as a function of local Ur and delayed Ur * . For Sk, it is seen that variability reduces when using Ur * (compare Figure 12a,c), which is also reflected by a small decrease in standard deviation (averaged of the standard deviations of the five bins) from 0.09 to 0.08. Although Ur * is optimized for Sk, As also seems to be better predicted by Ur * since the individual sensors for high |As| follow similar trends. The averaged standard deviation for As does not decrease (0.07 for Ur and Ur * ), which is probably because most points have close to zero As, thus will not be affected by the delayed parameterization. Visually, it can be seen, however, that it decreases the number of outliers.

Future Work
This paper characterizes conditions in which the wave shape deviates from commonly used parameterizations. It is proposed in a follow-up study to include offshore wave steepness ka o f f as this determines the shoaling characteristics and breakpoint of the waves [20]. In addition, the strength of the nonlinear interactions should be included since this proved to explain most of the variability at the shallow sensors for non-breaking waves. Once waves break, the parameterization can be improved by taking into account the proportion of waves which are breaking, because this was shown to be the dominant parameter defining the amount of As in breaking conditions. Finally, it should be investigated whether the delayed response of the wave shape to changes in the Ursell number under rapidly changing conditions should be accounted for, and if so how it varies with bed slope, wave length and directional spreading.
The efficiency of the above-mentioned improvements cannot be properly quantified with our field data due to their limited spatial resolution. Although it is an exceptionally comprehensive dataset, wave characteristics were only measured at ten locations over a relatively large area, with only one measurement point in the breaking zone and one in the de-shoaling zone. Future work will include the use of detailed wave-resolving simulations, from which Sk and As can be estimated with a high spatial resolution. This will allow us to study the development of wave shape at many more locations and under controllable conditions in which only a single input parameter is varied. This will allow us to further improve our understanding of the physical processes relevant to predict wave shape evolution and improve future parameterizations.
As mentioned above, parameterizations for Sk and As are usually used to parameterize wave-induced sediment transport in long term, process-based models. In such a modeling approach, wave predictions are usually based on a phase-averaged spectral wave model. It is interesting to note that such models provide much more information than the wave parameters currently used to calculate Ur and hence to parameterize Sk and As. They for instance also provide data on directional spreading, frequency spreading and the proportion of breaking waves and associated dissipation patterns, which could also be included in an improved parameterization.

Conclusions
This paper presents a new dataset of waves and currents obtained at 10 locations on the ebb-tidal shoal seaward of the Ameland Inlet. Pressure signals were transformed in near-bed orbital velocity signals in order to investigate the occurrence of Sk and As, as these are the main drivers of wave-induced sediment transport. Significant Sk and As were observed during the stormy periods with H m0 > 2 m, and it was seen that Sk and |As| increased moving to shallower water.
Wave shape parameterizations are usually a function of local Ur. For this dataset it was found that variability in wave shape was well-explained by variability in local Ur in deeper water. However, for shallower locations, variability in Ur could only partially explain the variability in Sk and As, indicating that some physical processes are not properly accounted for. The increased scatter at the shallow locations could be attributed to a combination of three main physical processes: wave-breaking, spatial variability in tidal currents, and the nonlinear energy transfer rate.
Towards the shallow sensors, a rapid decrease in Ur for breaking conditions and increase for non-breaking conditions is observed, although Sk and As only marginally change.
This leads to significant scatter in Sk and As as a function of Ur at these shallow locations (see Figure 9). For opposing currents, this decrease in Ur for breaking conditions is even bigger than for following currents, because the wave length decreases and increases, respectively. Another influence contributing to the scatter of Sk as a function of Ur is the nonlinear energy transfer rate, which is not properly accounted for in Ur.
It is shown that a potential improvement for future parameterizations would be to correct for the fact that the wave shape does not instantly respond to changes in depth and currents. For this dataset, this led to a reduced variance from 0.09 to 0.08 for Sk. It is expected that the added value will be higher if relatively more cases with high nonlinearity are included in the dataset.
Author Contributions: All co-authors contributed to the used methodology. F.d.W. performed post-processing and analysis of field data and wrote the manuscript. M.T. and A.R. contributed to supervising, editing, and improving this paper. water levels and depth-averaged velocities at Frames 1, 4, and 5. Table A1 shows that the order of magnitude of water levels and velocities is well-predicted by the model, with RMSE not exceeding 0.25 m and 0.13 m/s for water levels and depth-averaged velocities, respectively. Figure A1. Comparison between measured and modeled η, v east , and v north at Frames 1, 4, and 5 where black dots are individual data points and red dashed lines are unity. Table A1. RMSE between half-hourly-averaged modeled and measured water levels and depth-averaged velocities at Frames 1, 4, and 5.