Characterizing Wave Shape Evolution on an Ebb-Tidal Shoal

Abstract

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.

1. 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=\frac{3}{4}\frac{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.

2. Methodology

2.1. 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).

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. Overview of instruments, their location, the mean depth and bed slope they were in, and the sampling frequency they operated in.

Cluster	Location	Depth (m)	Bed Slope (-)	Instrument	Frequency (Hz)
Shelf	PS1	10.4	0.004	PS	10
	PS8	9.5	0.004	PS	10
	F4	8.5	0.004	ADCP down	4
				ADCP up	1.25
	PS3	8.2	0.002	PS	10
Seaward slope	PS5	7.9	0.005	PS	10
	F5	6.6	0.006	ADCP down	4
				ADCP up	1.25
	PS7	5.3	0.009	PS	10
Shoal	PS2	4.3	0.011	PS	10
Landward slope	PS4	4.6	−0.002	PS	10
Other	F1	6.6	0.009	ADCP down	4
				ADCP up	1.25

. 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.

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 (Figure 3a–c) were obtained from a wave-buoy (green circle in Figure 1), and tidal elevation and depth-averaged velocities (Figure 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].

2.2. 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.

2.2.1. 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 ${\theta }_U$.

2.2.2. 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+\langle p\rangle /\rho g$, in which $z_p$ is the height of the pressure sensor above the bed, $\langle p\rangle$ is the burst averaged pressure, $\rho$ is the water density and g is the gravitational acceleration. The surface elevation $\eta$ 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\sqrt{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 ${\theta }_w$ and spreading ${\sigma }_{{\theta }_w}$. Wave number k was computed by iteratively solving the linear dispersion relationship:

$$\begin{array}{c}\hfill \omega =\sigma +k{U}_n=\sqrt{gktanh\left(kd\right)}+k{U}_n\end{array}$$

(1)

with absolute frequency $\omega =2\pi /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, ${\theta }_U$, and ${\theta }_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.

2.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:

$$\begin{array}{cc}\hfill S{k}_{dim}& =\langle q^3\rangle \hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill A{s}_{dim}& =\langle H{\left\{q\right\}}^3\rangle \hfill \end{array}$$

(3)

where q is the filtered time series of a quantity, which can be either velocity or pressure, $H\{\dots \}$ denotes the imaginary part of the Hilbert transform, and $\langle \dots \rangle$ 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:

$$\begin{array}{cc}\hfill Sk& =\frac{\langle q^3\rangle }{{\langle q^2\rangle }^{3/2}}\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill As& =\frac{\langle H{\left\{q\right\}}^3\rangle }{{\langle q^2\rangle }^{3/2}}.\hfill \end{array}$$

(5)

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.,

$$\begin{array}{c}\hfill Ur=\frac{3}{4}\frac{a}{k^2{d}^3},\end{array}$$

(6)

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.

2.2.4. 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 $\sigma$ 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:

$$u_p=\sum _i{p}_i\frac{{\sigma }_i}{\rho gtanh\left(k_{i}d\right)},$$

(7)

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.

2.3. 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 $k{U}_n$ compared to the magnitude of $\sigma$. 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 $\sigma$, 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, $\sigma$, and $Ur$ changed significantly when including the currents.

3. Results

3.1. Data Overview and Selection

Significant wave height $H_{m0}$, dimensional skewness $S{k}_{dim}$ and asymmetry $A{s}_{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 $S{k}_{dim}$ and $A{s}_{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 $S{k}_{dim}$ increases while moving to the shallower sensor, and a clear tidal pattern can be seen in the $S{k}_{dim}$ during the second storm. Large values of $|A{s}_{dim}|$ occur when conditions are close to wave breaking or once the waves are breaking. This explains why near-zero $A{s}_{dim}$ was observed at PS1, but that waves get asymmetric at PS2.

For the following analyses, the focus was on cases with considerable $S{k}_{dim}$ and $A{s}_{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 ${\langle u^2\rangle }^{3/2}$ is really small. A combined nonlinear parameter $B_{dim}=\sqrt{S{k}_{dim}^2+A{s}_{dim}^2}$, 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.

3.2. 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. Correlation coefficient $R^2$ between measured $Sk$ and $As$ and $Ur$.

Location	$\mathit{Sk}$${\mathit{R}}^2$	$\mathit{As}$${\mathit{R}}^2$
PS1	0.89	0.03
PS8	0.88	0.05
F4	0.87	0.10
PS3	0.89	0.07
PS5	0.91	0.29
F5	0.87	0.47
F1	0.83	0.36
PS7	0.76	0.74
PS4	0.28	0.01
PS2	0.13	0.05

). 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$.

$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 $\left|As\right|$ 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.

3.3. 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 $\left|As\right|$-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.

3.3.1. Variability Caused by Wave Breaking

High $\left|As\right|$-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=\frac{H_{m0,1}\sqrt{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:

$$c_{g_i}=(\frac{1}{2}+\frac{k_i{d}_i}{sinh\left(2k_i{d}_i\right)})\frac{{\sigma }_i}{k_i}.$$

(9)

If G is higher than 0.5, which is within the common range of breaking criteria for irregular waves over sloping bathymetry for $\gamma =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 $\left|As\right|$ stays near-zero for non-breaking cases and slightly increases for breaking cases (Figure 9c). The increase in $\left|As\right|$ 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.

3.3.2. 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 $\left|As\right|$ 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).

3.3.3. 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$.

4. Discussion

4.1. 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.

4.2. 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, $U{r}^{*}$ (defined at depth $d^{*}=d+\Delta d$, with $\Delta 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 $\Delta 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 $\Delta 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 $U{r}^{*}$ is computed as:

$$U{r}^{*}=\frac{3}{4}\frac{a^{*}}{k^{*2}{d}^{*3}},$$

(10)

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. Depth by which $Ur$ needs to be delayed to get $U{r}^{*}$ per sensor accompanied by the mean depth and bed slope the sensors are in.

Location	PS1	PS8	F4	PS3	PS5	F5	PS7	PS2
$\Delta d$ (m)	0	0.25	0.2	0.2	0	1.1	0.8	0.85
d (m)	10.4	9.5	8.5	8.2	7.9	6.6	5.3	4.3
bed slope (-)	0.004	0.004	0.004	0.002	0.005	0.006	0.009	0.011

shows per location the $\Delta d$ by which the Ursell number needs to be corrected in order to obtain the optimized $U{r}^{*}$. It is seen that the delay is negligible ($\Delta 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 $U{r}^{*}$. For $Sk$, it is seen that variability reduces when using $U{r}^{*}$ (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 $U{r}^{*}$ is optimized for $Sk$, $As$ also seems to be better predicted by $U{r}^{*}$ since the individual sensors for high $\left|As\right|$ follow similar trends. The averaged standard deviation for $As$ does not decrease (0.07 for $Ur$ and $U{r}^{*}$), 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.

4.3. 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 $k{a}_{off}$ 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.

5. 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 $\left|As\right|$ 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.