Shoaling Wave Shape Estimates from Field Observations and Derived Bedload Sediment Rates

: The shoaling transformation from generally linear deep-water waves to asymmetric shallow-water waves modiﬁes wave shapes and causes near-bed orbital velocities to become asymmetrical, contributing to net sediment transport. In this work, we used two methods to estimate the asymmetric wave shape from data at three sites. The ﬁrst method converted wave measurements made at the surface to idealized near-bottom wave-orbital velocities using a set of empirical equations: the “parameterized” waveforms. The second method involved direct measurements of velocities and pressure made near the seabed: the “direct” waveforms. Estimates from the two methods were well correlated at all three sites (Pearson’s correlation coefﬁcient greater than 0.85). Both methods were used to drive bedload-transport calculations that accounted for asymmetric waves, and the results were compared with a traditional excess-stress formulation and ﬁeld estimates of bedload transport derived from ripple migration rates based on sonar imagery. The cumulative bedload transport from the parameterized waveform was 25% greater than the direct waveform, mainly because the parameterized waveform did not account for negative skewness. Calculated transport rates were comparable to rates estimated from ripple migration except during the largest event, when calculated rates were as much as 100 times greater, which occurred during high period waves. cumulative bedload calculated using the parameterized method was about 25% greater than that calculated using the direct method. These results demonstrate the efﬁcacy of using the parameterized method to obtain representative waveforms for use in asymmetric bedload formulations. The direct measurements revealed instances of negative skewness at all three sites, but the equations used for the parameterized method cannot generate negative skewness. Episodes of negative skewness drive offshore transport, which explains why the parameterized approach predicts a larger onshore wave-driven bedload than the direct method. While this work shows the efﬁcacy of Ruessink’s method in estimating asymmetric waveform characteristics, that approach relies on local measurements and does not account for shoaling history Other studies have included non-local and offshore wave parameters to estimate wave asymmetry. They studied the effects of offshore wave steepness, offshore spectral bandwidth, and beach slope and showed a 50% better prediction of asymmetric parameters nearing breaking regimes. Incorporating the effects of negative skewness in parameterizations and testing the effect of adding non-local wave conditions to estimate asymmetric waveforms could potentially improve the pre-dictions of bedload. In addition, more observations at shallower depths where breaking waves lead to higher acceleration asymmetry could improve the empiricisms involved in the asymmetric bedload formulation, which is currently limited to non-breaking waves. Author Contributions: Conceptualization, C.R.S. and J.C.W.; methodology, C.R.S., T.S.K. and S.E.S.; software, T.S.K., S.E.S., C.R.S. and G.R.L.; validation, T.S.K. and S.E.S.; Analysis, T.S.K., S.E.S., A.L.A. and C.R.S.; writing—original draft preparation, T.S.K. S.E.S.; writing—review and T.S.K., A.L.A.,


Introduction
Waves propagating towards the coast undergo transformations from their sinusoidal wave shape in deep water to non-linear wave shape in shallow water. "Skewed", as used here, describes wave shapes, bottom-orbital velocities, and sediment transport generated by waves that cannot be adequately represented by linear (Airy) wave theory (e.g., [1]). The fundamental equations for continuity of mass and momentum used as the starting point for most wave theories are linear, but the surface boundary conditions are not. Linear wave theory neglects the non-linear components by assuming that the ratio of wave height/wavelength (wave steepness) is small. As waves shoal and steepen, this assumption no longer holds, and non-linear equations are needed to accurately represent the waves [1]. Outside the surf zone, this shoaling transformation causes the near-bottom wave orbital velocity to differ between the wave crest and trough cycles, referred to as velocity skewness [2]. As waves propagate further into the surf zone, the differential duration in crest and trough half cycles can lead to a differential acceleration in each half cycle, referred to as acceleration asymmetry [3].
The effect of shoaling waves on resulting sediment dynamics has been quantified in various studies. In previous flume experiments [4][5][6], skewed waves led to onshore sediment transport near the bed. Under mild wave conditions, velocity skewness has been found to be the main driver of onshore sediment transport, generating onshore sandbar migration [7][8][9]. The combined effect of different asymmetric wave shapes in the presence of currents acting against the direction of waves and showed that under decreasing skewed-asymmetric wave shapes, net offshore sediment transport increased [10]. Onshore sediment transport has been also related to ripple migration associated with velocity skewness in several field observations [11][12][13]. Other than flume and field data, numerical models [14][15][16] have shown that acceleration asymmetry caused high onshore sediment transport and contributed towards beach recovery. In addition, process-based morphodynamic modeling studies have shown the impact of wave-induced cross-shore sediment transport that caused inlet closure during weak inlet current conditions [17][18][19]. In another study [20], a numerical model in an idealized inlet domain that accounted for asymmetric waves showed that larger waves contributed to the migration and growth of shoals on ebb-tidal deltas. Using a realistic domain [21], it was shown that the skewness of the waves was responsible for the enhancement of sediment transport into the inlet and the evolution of the ebb-tidal delta during storm conditions. Recent modeling [22][23][24][25][26] studies have shown the importance of wave-driven transport under combined action during wave-current interaction.
These studies highlight the importance of estimating the asymmetric waveforms that affect nearshore sediment transport. Other efforts [27,28] emphasized that understanding wave-driven sediment transport is fundamental to the improvement of morphodynamic models. In the presence of waves, morphodynamic models need to account for both suspended and bedload transport. Of the two modes of sediment transport, bedload transport is most affected by near-bottom wave orbital wave motions [29], thus necessitating the need to estimate asymmetric waveforms in modeling frameworks. Asymmetric waveforms can be computed by wave-resolving models [30,31], but these models are computationally expensive and impractical for regional-scale morphodynamic models. One approach to estimating asymmetric waveforms used a statistical method [32] that applied a genetic algorithm to obtain velocity skewness based on field observations made under shoaling and breaking wave conditions. Another approach to quantify asymmetry from shoaling waves is to use the bulk wave statistics and local depth [33,34]. The approach mentioned in [34] would be referred to as Ruessink's method. A recent study [35] compared the modified version of the approach mentioned in [33] to Ruessink's method [34] to simulate long-term morphodynamic change. They [35] found that Ruessink's method provided a better representation of wave and near-bed orbital velocities, required less calibration and, therefore, resulted in more realistic morphodynamic predictions. A modification to Ruessink's approach was used in a study [36] to accelerate the model calculations to capture the effects of nourishment on cross-shore bedload fluxes at different locations on an idealized planar beach. Ruessink's parameterization also showed a good skill predicting shoreward evolution of wave skewness and asymmetry under high-energy waves ( H s ∼ 2 m) during offshore sandbar migration events [37].
In the present work, we analyzed the efficacy of the widely used Ruessink's parameterization by taking observations from three sites located outside the surf zone. At these sites, conditions (including depth and wave characteristics) only generate velocity skewness in the waveform shape; acceleration skewness becomes more important in shallower water. The three sites were located at Fire Island, New York, Martha's Vineyard Coastal Observatory, Massachusetts, and Matanzas Inlet, Florida. The efficacy of Ruessink's method was determined by parameterized estimates computed from bulk surface-wave statistics with direct observations of skewed waveforms made near the bottom. While the parameterizations depended on bulk wave statistics (significant wave height, representative wave period) and water depth, and were translated to the bottom using linear wave theory, the direct observations were obtained from measurements of near-bottom velocities, as described later in this work. The characteristics of the skewed waveforms from both methods include: crest cycle velocity, (u c ), trough cycle velocity, (u t ), crest cycle period, (T c ), trough cycle period (T t ), and accelerative periods of crest and trough cycle, (T cu ) and (T tu ), respectively. Following this comparison, the skewed waveforms obtained from the parameterized and direct methods were used to drive a bedload-transport formulation that accounted for asymmetric waveforms [38]. We compared the bedload transport results with the widely used Meyer-Peter and Müller (MPM) equation [39] and bedform-transport rates from sonar data that tracked the ripple migration on the seafloor. The bedload-transport rate comparisons provide an insight into the implementation of asymmetric bedload formulation in morphodynamic models that would rely on Ruessink's method to first obtain the asymmetric waveform in a way that is analogous to the "parameterized" calculations described in the present observational work. The paper is organized as follows: The dataset, instrumentation, and analytical methods used to obtain the skewed waveforms from parameterized and direct methods along with the methods used to obtain bedload are described in Section 2. The results are presented in Section 3, followed by a discussion in Section 4. Section 5 presents the conclusions.

Field Study Sites
The datasets used in the present work were obtained at three field study sites: Fire Island, New York (FI), Martha's Vineyard Coastal Observatory, MA (MVCO), and Matanzas Inlet, FL (MI). The FI dataset was part of a field campaign [40,41] to study the along-shore variations in waves and sediment processes. FI data from the site shown in Figure 1a correspond to a mean depth of 16 m located approximately 2 km offshore from the southern shoreline of western Fire Island in New York. The MVCO dataset [42] was collected at a site with an average depth of 14 m offshore located approximately 2 km from the southern shoreline of Martha's Vineyard (Figure 1b). The Matanzas Inlet dataset was collected [43] to study the effects of cross-shore currents on waves, and measurements were made at a site with an average depth of 9 m located approximately 0.8 km offshore from the eastern shoreline ( Figure 1c). Table 1 mentions the filenames containing the datasets used to obtain near bottom and surface wave measurements.  The parameterized method to compute the near-bottom wave forms uses surface bulk wave parameters. The parameterized method to find skewed waveform characteristics (u crest , u trough , T c , T t , T cu , T tu ) is described in Appendix A and requires local values of significant wave height H S , a representative wave period T, depth h, and wave orbital velocity amplitude u w . Surface wave spectra were measured with an acoustic Doppler current profiler (ADCP; Table 1). The ADCP beam orbital velocities measured at bins near the surface were used to calculate wave spectra at the three sites using the WavesMon software. An exception to this methodology was the Matanzas Inlet site where, after the first 31 days of the 80-day deployment, measurements from a near-bed pressure sensor (SeaBird Seagauge) were substituted because the ADCP malfunctioned. In these pressure-sensor data, the mean H S was 7% smaller and the mean T p was 8% larger than those measured with the ADCP for the 31 days when both the sensors worked. The methods described by [44] were used to obtain near-bed estimates of u w and T from the surface wave spectra. Note that for the application of Appendix A, we used a representative wave period (T) that was defined as T = m 0 m 1 , where m 0 is the variance of water surface elevation and m 1 is the first moment of the spectral density of surface elevation. Table 2 provides sensor specifications used from each of the three sites. The characteristics of the parameterized waveform obtained using this method are described later with the use of subscript "(p)". Table 2. ADCP sensor specifications for the three sites. The sensor height for the velocity measurements changed during the deployment period but the median elevation remained around the nominal elevation (latitude and longitude in decimal degrees).

Direct Method: Skewed Waveform Estimation from Near-Bed Data
The direct method uses near-bed observations to determine a representative skewed waveform. Waveforms were observed directly from near-bed velocity measurements made with acoustic Doppler velocimeters (ADVs; Table 3). The measured time series of wave velocity directly obtained from ADV recordings are in two components, u raw and v raw close to the seafloor in bursts (ref. Table 3 for site-specific burst length) at an elevation close to the seafloor but above the wave boundary layer. The process is described here for one component, that is, u raw for the sake of brevity; a similar analysis was performed for v raw . u raw was decomposed into a slowly varying mean tidal component u , an incident-wave band u incident , an infragravity wave band u in f ragravity and noise associated with turbulence and instrument error u . Velocities from each burst time series were first detrended to remove u . Incidentwave velocities u incident with periods of 4-20 s were isolated using a Butterworth bandpass filter. The lower limit of 4 s removed higher-frequency motions that were assumed to not affect the calculations of bedload at the bottom, thus reducing the noise in the signal. The upper limit of 20 s excluded waves in the infragravity band. The infragravity band was removed so that only the sea-swell frequency band remained, to match the frequency range in Ruessink's data [34] and to apply an asymmetric bedload-transport equation that was formulated with the exclusion of the infragravity band. Figure 2 compares the raw and incident wave spectrum from a storm event in the Fire Island dataset that corresponds with an H S of 3.7 m associated with a corresponding T P of 11.6 s.
(a)  A principal components analysis (PCA) was applied to the u incident and v incident wave signals to determine the axis of dominant wave motions. PCA determined the dominant wave direction and the corresponding velocity signal (referred to as u) for each burst. Representative waveforms were found within burst time series by (a) identifying individual waves using zero-up crossings, (b) determining characteristics (u crest , u trough , T crest , T trough , T cu , and T tu ) of the individual waveforms, (c) calculating parameters that describe waveform according to the formula of [3], and computing a weighted average of those parameters. Zero up-crossings times were used to find short time series of velocity u(t) with period T i for each wave in the time series, where subscript "i" refers to the individual waveforms in a wave burst time series. From these individual wave motions, the crest and trough velocities u c,i , u t,i and crest and trough periods T c,i and T t,i were determined. In the next step, three parameters corresponding to the individual waveforms were calculated that included the amplitude of wave orbital velocity as defined by: and their velocity skewness (R i ) according to [45]: The inputs of non-dimensional u w,i , T i and R i were used to calculate the third parameter b 1,i a fitting parameter that results in a free-stream velocity description [3] by finding the maximum and minimum positions of velocity in the individual waveforms. The equations to calculate b 1,i are provided in the Matlab script "abreu_fit" provided by [46]. The weighted means of individual waveform u w,i , T i , R i , and b 1,i were used in the Matlab script "abreu_fstream" [46] to obtain the representative waveform for each burst time series (Figure 2b as an example). Weighting was based on the product of the squared wave-orbital amplitude and wave period, as: The characteristics of the parameterized waveform obtained using this method are described later with the use of subscript "d".

Calculation of Bedload for Skewed Waveforms
We derived three estimates of wave-driven cross-shore bedload transport using the Matanzas Inlet dataset because that site had the highest proportion of skewed waveforms (Section 3.4). Both the parameterized and direct skewed waveforms were used to drive a bedload transport formulation that accounted for skewed waves [38], referred to as asymmetric bedload formulation. Results were also calculated using a traditional excess-stress bedload formulation, the Meyer-Peter and Müller (MPM) equation [39]. Appendices B and C detail the equations used for the asymmetric and MPM formulations. The grain size used in the calculations was d 50 = 0.23 mm [43]. In contrast to the asymmetric bedload formulations, which represent waveforms with wave asymmetry indices, the MPM formulation (Equation (A43)) uses the entire time series of instantaneous wave motions ( u incident ) in the bedload calculations.

Sonar Imagery for Bedform Transport Rate from Ripple Migration
Measurements of bedform transport, a proxy for bedload transport, were determined from sonar imagery that tracked ripple migration on the seafloor. The analysis was performed for the shallowest site (Matanzas Inlet) where good images of active ripple migration were available. Images of the seabed were obtained every 30 min from a 2.25 MHz Imagenex model 881A fan-beam rotating imaging sonar mounted 0.45 m (nominal) above the bed. The raw images of backscatter intensity from the sonar were corrected for slant range and rotated into geographical coordinates (positive x = east). This resulted in images with a horizontal resolution of 0.01 m and covering a circular region with a 5 m radius ( Figure 3a). In these images, the light and dark pixel values corresponded to acoustic backscatter intensity (light indicates high backscatter that revealed the presence of ripples). Manual inspection was used to find a region with ripples that were well defined for most of the time series. This corresponded to a 1-m × 1-m square (Figure 3a) orientated in the onshore direction located 1-2 m from the center of the image. The ripple wavelength (λ r ) was determined in the sub-patch region using the two-dimensional fast Fourier transform method described in these studies [47][48][49][50].
We used a single cross-shore transect (shown in Figure 3a) located at 0.05 m distance from the southern edge of the square to form a time-stacked image (Figure 3b). Manual trials showed that the quality of the time stacked image remained unchanged for ripple tracking purposes by choosing any other current transect within the square region. The cross-shore position of pixel intensity was determined using edge detection [51,52]. The cross-shore migration rate (d r ) was determined as the change in cross-shore position at half-hour intervals. The resulting volumetric bedform transport rates ( → q b , m 3 /m) estimated from the sonar data were computed [11] as: assuming a bed porosity factor ε = 0.4, dimensionless ripple shape factor f r = 0.5, and a constant ripple height to ripple wavelength ratio η r λ r = 0.16 [11].

Surface Wave Statistics
MVCO had the least-energetic waves of the three sites ( Figure 4). Fire Island was the deepest site and had the highest proportion of high H S . Matanzas Inlet had the highest proportion of waves with high T P . Energetic waves (high H S and T P ) and shallower water depths would be expected to generate larger wave skewness.

Comparison of Representative Wave Period and Orbital Velocity from Surface Wave Data and Near-Bed Wave Data
Peak wave orbital velocity (u w,d and u w,p ) and representative wave periods (T w,d and T w,p ) are compared in Figure 5 for the three sites, based on both (parameterized and direct) methods described in Sections 2.2 and 2.3 where subscript "w" refers to the representative waveforms obtained from the two methods (parameterized and direct). The correlation (Pearson's correlation coefficient, r) between u w,d and u w,p was above 90% at all three sites and the best linear fit (identified by the slope value) was obtained at the shallowest site (Matanzas Inlet; Figure 5a,c,e). At Fire Island and Matanzas Inlet, u w,p was underpredicted compared to u w,d as significant wave heights (H s ) increased (Figure 5a,e). The correlation between T w,d and T w,p was found to be above 90% at all three sites (Figure 5b,d,f). T w,p was consistently underpredicted compared to T w,d at all three sites for almost all wave heights. This pattern can be attributed to the damping of high-frequency motions in the near-bottom measurements, leading to relatively higher values of T w,d . The linear fits to T w,d and T w,p were better in the shallower depths at MVCO (Figure 5d) and Matanzas Inlet (Figure 5f), compared to those at Fire Island (Figure 5b).

Skewed Waveform Characteristics: Parameterized and Direct Methods
We used the equations described in Appendix A and methods described in Sections 2.2 and 2.3 to determine the parameterized and direct representative skewed waveform characteristics which included u crest , u trough , T c , T t , T cu , and T tu . These characteristics were used to calculate the wave-driven bedload calculations (Appendix B). The comparison of u crest,d and u crest,p ( Figure 6a) and u trough,d and u trough,p (Figure 6b) showed a r > 0.9 and that linear fit improved for the Matanzas Inlet site compared to the other two sites. The correlations between T c,d and T c,p , T t,d and T t,p , T cu,d and T cd,p and T tu,d and T tu,p remain high (r > 0.85) (Figure 6c-f). The pattern of linear fit for half periods was consistent with the results of the representative wave period (i.e., Figure 5b,d,f).
The E 1 error was highest for the shallowest site of Matanzas Inlet for all waveform characteristics (Table 4). E 2 error that accounted for scaling the error with the RMS values showed that the Matanzas Inlet site had the least error corresponding to the velocity characteristics of the waveform (u w , u crest , u trough ). On the contrary, E 2 error was highest for period characteristics (T, T c , T t , T cu , and T tu ) at the Matanzas Inlet site.

Wave Skewness
Wave-orbital velocity skewness drives bedload transport and can be used as a diagnostic tool to indicate transport potential. Non-dimensional wave-orbital velocity skewness (hereafter referred to as skewness, S u ) was calculated at the bottom using the ADV-based direct method as: where u is wave-orbital velocity, angle brackets indicate time-averaging over the burst, and subscript rms indicates root mean square over the burst. S u tends to increase with Ursell number (U r ) above a threshold of about U r = 0.04 (Figure 7). This is consistent with the observations of [34] who mentioned that U r ≤ 0.04 have near-zero skewness. Most of the observations at the two deeper sites (93% at FI and 75% at MVCO) fell below this threshold.
In contrast, only 25% of the MI data had a U r ≤ 0.04. For U r > 0.04, about 22% of the waveforms at MI were associated with negative skewness (similar to the results in [34]) (Figure 7). We examined the directional wave spectrum, D spec = f (ω, θ) of bottom skewness, where ω is the frequency in Hz and θ is the direction in degrees to determine whether negative skewness was related to spectral characteristics. The events associated with positive and negatively skewed waves (S u ) in the directional wave spectrum were binned separately. We first averaged D spec along all directions to get mean frequency spectra ( D spec−ω , Figure 8). The spectra associated with positive S u (green shaded region) frequencies are shifted toward lower frequencies compared to spectra associated with negative S u (red shaded region), which are shifted towards higher frequencies at all three sites (Figure 8a-c). Positively skewed waves also had narrower frequency distributions. The differences in mean frequency and standard deviation about the mean for spectra with positive versus negative skewness were significant with (p < 0.01). We also found that the statistical skewness of the spectral distributions was higher for positively skewed waves (significantly different with p < 0.001). This pattern confirmed a shift of the spectral distribution towards lower frequencies associated with positively skewed waves. Next, we averaged D spec across all frequencies to get spectra D spec−θ varying over direction and binned based on the sign of skewness (Figure 9a-c). Positively skewed waves (green lines) were associated with a smaller directional spread of energy while negatively skewed waves (red) corresponded to a larger directional spread of energy (Figure 9a-c). This was confirmed by the kurtosis of the direction distributions, which was significantly different (p < 0.001) with positively skewed waves exhibiting more peakiness (leptokurtic), and negatively skewed waves exhibiting broader directional spread, especially at MVCO. In summary, data from all three sites (Figures 8 and 9) show that spectra associated with positively skewed waves were distributed over a lower peak frequency with a narrower directional spread. Conversely, spectra for negatively skewed waves were distributed over higher frequencies and with a larger directional spread of energy at all three sites. A similar pattern was noticed in data [12] from an 18-h storm event. They concluded that the sign of the skewed waves originated from wave-wave interactions, and the spectral shift towards higher frequencies originated from mixed sea and swell conditions, whereas unimodal swell conditions with lower frequencies led to positive skewness, consistent with Figure 8.

Estimating Wave-Driven Cross-Shore Bedload Using Skewed Waveform Parameters
We compared our three estimates of wave-driven cross-shore bedload transport (Section 2.4) using the Matanzas Inlet dataset because that site had the highest proportion of skewed waveforms (Section 3.4). The first two estimates were calculated using a slightly modified version of the SANTOSS formula ([38]; Appendix B). This formula requires estimates of wave asymmetry and calculates transport separately for each part of the wave cycle. We derived asymmetry values from surface waves (our "parameterized" method using the [34] approach; Appendix A) and from near-bottom measurements (our "direct" method; Appendix B) and compared it with a traditional excess stress-based bedload formulation (i.e., the MPM method detailed in Appendix C). Figure 10a,b compare instantaneous (q b ) and time-integrated net ( t 0 q b dt) bedload transport rates where positive values in Figure 10a indicated bedload along the wave aligned axis and negative values are opposite to that. The bedload values reach local peaks during similar time intervals between all three methods and are associated with high wave energy events (refer to Figure 4) that result in larger skewness (Figure 10a). The parameterized method's q b exceeded the direct one during 51% of time instances because it always predicted a positive skewness and as a result, the parameterized q b showed an increase of 25% in the ( t 0 q b dt) (Figure 10b). In both the direct and MPM method, there is a change in the sign of q b during certain periods although the magnitude of negative q b is much smaller over the entire time series for both methods. For all the events where q b > 0, the ( t 0 q b dt) with MPM method was 700% smaller when compared to the direct skewed waveform based (  Analysis of the time-stacked ripple images identified seven periods that were continuous, which allowed tracking of the ripple front. Manual interpretation (red lines, Figure 11a) was used to detect the starting and ending point of each continuous period. The ripple wavelength in these periods ranged from 0.065 to 45 cm. Using the distance traversed by ripple wave fronts over the duration of each segment, we found ripple migration rates ranging from 0.3 cm/day to 0.48 cm/day. The corresponding significant wave height (H s ) and representative wave period (T) at the bottom ranged between 0.4 to 1.5 m and 5.8 to 14.4 s respectively (Figure 11b,c). The Ursell (U r ) number during these periods ranged from 0.024 to 0.23 (Figure 11d). The relatively low values of these parameters (Figure 11b-d) for the periods when ripples could be tracked suggest the sonar was only capable of providing bedform transport estimates during weak events. During other times, ripples were not clearly visible in the images.
A comparison of instantaneous q b from sonar imagery and empirical methods (Figure 11e; note logarithmic y-axis) shows that there was substantial scatter within events, and that the four methods often overlapped. The q b estimates are typically within 2 orders of magnitude with each other. The exception is the strongest event, between 06-09 March, when the sonar q b estimates were lower than q b calculated from both the asymmetric waveform methods and the MPM method. There are periods (02-03 February, 06-07 February, and 03-04 April) when the parameterized method results in zero q b (not visible on the logarithmic axis) because the calculated stress did not exceed the threshold critical stress in crest or trough cycles. The RMSE of the cumulative bedload calculated from the asymmetric bedload method is lower than MPM for two events (1 and 7) but in all other events, it is greater (Table 5). For the asymmetric bedload, events 1 and 2 are significantly (positively) correlated (r 1 ) with the sonar data explaining between 56-62% of variance. Events 4 and 5 also have a positive correlation but are not significant. For MPM, event 2 is the only one for which the explaining variance is above 50%. Table 5. Event-based cumulative error (RMSE) and correlation coefficient (r 1 ) between the calculated transports (both asymmetric direct and MPM) and the sonar estimates. Significant correlations (p < 0.05) are indicated in bold.

Event Period
Asymmetric (

Skewness Parameterization for Operational Morphodynamic Models
In observational datasets, the near-bottom velocity skewness of the wave shape is often calculated with Equation (6) [53] and u 3 , can be used to estimate a net sediment transport rate driven by velocity-skewed oscillatory flows [4]. At the three study sites, parameterized S u is always positive while the direct skewness changes sign due to the factor of u 3 (also seen in Section 3.4) (Figure 12). Observations [12] showed the dependence of bedform migration direction on the sign of skewness where positively skewed waves resulted in an onshore migration of ripples and negatively skewed waves resulted in the offshore migration of ripples. However, in this study, the parameterized method that depends on Ruessink's approach always predicts positively skewed waveform shapes, indicating onshore wave-driven sediment transport. Therefore, the application of this methodology in models may overpredict onshore wave-driven transport and require one to calibrate velocity skewness based on field data and/or modify the current and wave contributions [35].

Bedload Comparison from Various Methods
The large differences between the cumulative bedload ( Figure 10b) from asymmetric bedload formulations and the MPM method may be attributed to several factors, all of which enhance transport in the direction of wave propagation. The asymmetric bedload formulation was designed to account for differential crest and trough cycle stresses and for processes associated with free surface wave propagation. These processes included the effect of modifying crest and trough periods due to surface waves that led to a longer crest period and shorter trough period experienced by sediments at the bed ( [54], Equations (38) and (39) in Appendix B of this paper). Another important process is the increase in crest cycle stress and a decrease in trough cycle stress due to surface wave propagation (Equation (40)), also referred to as "wave streaming" or "boundary layer streaming" [38, 55,56]. Using a high-resolution CFD model for waves and sediment that included a free surface [57] showed that progressive wave streaming produced a 60% increase in onshore sediment transport under nonbreaking wave conditions (only velocity skewed). Another mechanism that affects the asymmetric bedload calculation under surface waves is the adjustment of phase lag parameters under crest and trough to account for the differential horizontal sediment advection [38,58]. This mechanism leads to a decrease in effective settling velocity of sediments during crest (decreasing phase lag during crest) and an increase during trough (increasing phase lag during trough). However, the biggest difference between the parameterized method for wave shape used in conjunction with the non-linear transport formulation is the lack of negative skewness.
There was only limited sonar imagery suitable for extracting bedform transport rates for comparison with bedload calculations. The low pixel intensities and loss of bottom reflectors during storms due to plane-bed conditions or suspended material compromised most of the record. We explored the sensitivity of our migration-rate estimates to the starting point in the first two events (02-03 February and 06-07 February) in the time stack (Figure 11e). The resulting variations in migration rates varied as much as a factor of three, which is small compared to the order-of-magnitude changes in our bedform-transport estimates.
The choice of parameters (bed porosity factor, ripple shape factor, and ripple height to ripple wavelength mentioned in Section 2.5) would modify the results but not impact the order of magnitude changes in the sonar-based bedform transport. In future experiments, the ripple imagery could be improved by averaging over several sweeps to obtain clearer images of the ripples. This could reduce uncertainty in crest tracking and might also reduce the effort of manual interpretation and identification of the start and end points of migrating ripples in the time-stacked image. Despite the shortcomings, the dataset provides seven instances for comparison during weak to mild wave conditions and shows that transport rates from the asymmetric bedload formulations fall within an order of magnitude of rates from sonar data. This study is the first to compare and show that the asymmetric bedload formulations fall in the range of sonar data for all events (except event 5, likely due to the larger period waves). Overall, these results highlight the challenge in assessing the relative performance of bedload formulations and comparing it with field observations.

Conclusions
The use of skewed waveforms in sediment-transport calculations represents an advance in computing bedload transport that accounts for wave-driven processes in morphodynamic models. The present work compares parameterizations of skewed waveforms from bulk surface-wave statistics with waveforms obtained from direct measurements near the seafloor, using observations from three different sites under relatively calm conditions (u W < 0.76 m/s). Overall, the correlations for all the waveform characteristics (velocity and period in each half cycle) were significant and remained above 0.85. The two methods produced skewed waveform characteristics that matched well at the three sites. The cumulative bedload calculated using the parameterized method was about 25% greater than that calculated using the direct method. These results demonstrate the efficacy of using the parameterized method to obtain representative waveforms for use in asymmetric bedload formulations. The direct measurements revealed instances of negative skewness at all three sites, but the equations used for the parameterized method cannot generate negative skewness. Episodes of negative skewness drive offshore transport, which explains why the parameterized approach predicts a larger onshore wave-driven bedload than the direct method. While this work shows the efficacy of Ruessink's method in estimating asymmetric waveform characteristics, that approach relies on local measurements and does not account for shoaling history [29]. Other studies [59,60] have included non-local and offshore wave parameters to estimate wave asymmetry. They [60] studied the effects of offshore wave steepness, offshore spectral bandwidth, and beach slope and showed a 50% better prediction of asymmetric parameters nearing breaking regimes. Incorporating the effects of negative skewness in parameterizations and testing the effect of adding non-local wave conditions to estimate asymmetric waveforms could potentially improve the predictions of bedload. In addition, more observations at shallower depths where breaking waves lead to higher acceleration asymmetry could improve the empiricisms involved in the asymmetric bedload formulation, which is currently limited to non-breaking waves.  Table 1 specifies the filenames for accessing the data used in the above analysis. The datasets have been published on https://stellwagen.er.usgs.gov/, which is a Trusted Digital Repository for distributing time series data. These files can be downloaded from the following weblinks: 1. Fire Island [40,41]: https://stellwagen.er.usgs.gov/FIREISLAND14-a.html. 2. MVCO [42]: https://stellwagen.er.usgs.gov/mvco_15-a.html. 3. Matanzas Inlet [43]: https: //stellwagen.er.usgs.gov/matanzas_18-a.html. All weblinks were last accessed on 29 January 2022. Acknowledgments: We thank Kai Parker at the U.S. Geological Survey (USGS), St. Petersburg, Florida for providing his feedback to improve the quality of the paper. We thank Christie Hegermiller, Tom Hsu and Yashar Rafati for valuable discussions throughout the length of the project. We would like to acknowledge our colleague Jin-Si Over for her contribution in creating Figure 1 of the manuscript. We extend our gratitude to the team of scientists and engineers at USGS, Woods Hole Coastal, and Marine Science Center for gathering, processing, and publishing the observational datasets used in this work. The work by Tarandeep Kalra and Gibson Leavitt was done under contract to the USGS.

Conflicts of Interest:
Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Appendix A. Determining Skewed Waveform from Bulk Wave Statistics
The inputs required to determine the skewed waveform parameters from the wave spectrum include H s (significant wave height), T (representative wave period), u w (amplitude of wave orbital velocity), and h (water depth). Using these inputs, one can calculate Ursell number (U r ) [61] as: where a w = 0.5H s and k is the wave number computed from the dispersion relationship [39] that depends on T and h.
The third moments of skewness S u and asymmetry A u from 34,000 time series of natural near-bed orbital motion were computed according to [34]. They combined S u and A u into a non-dimensional skewed parameter B and found the phase angle ψ between skewness parameters depended on U r . B and ψ as a function of U r can be expressed as: where p 1 = 0.0, p 2 = 0.85, p 3 = −0.471, p 4 = 0.297, p 5 = 0.815 and p 6 = 0.672. S u and A u can be computed from B and ψ as follows: The analytical expressions [3] for time series of velocity and acceleration by curve fitting using parameters r and ∅.
where u w is the amplitude of wave orbital velocity and ω is the angular frequency expressed as: r is an index of nonlinearity, ∅ is a waveform parameter whose values determined the profile of the waveform (sawtooth shaped, 1st-order cnoidal wave etc.) and t ranges from 0 ≤ t ≤ 2π. To obtain r and ∅, we utilized the skewed parameters, B and ψ obtained above from [34], written as: where ∅ is defined as: where the parameter b Ruessink is computed from B as: A detailed discussion [3] on the implications of the choice of parameters r and ∅ to derive various nonlinear wave profiles. While the velocity in crest and trough along with skewness parameters can be obtained from analytical expressions [3], a more consistent way from a model implementation perspective was to use the formulations mentioned in [46]. From here onwards, the next set of equations used in step 1 are described in [46] that replace r by b (different from b Ruessink ) and utilize ∅ with its sign reversed to get the time series for velocity and acceleration and are written as: where b can be written as: where parameter P depends on r as follows: By knowing the maximum and minimum values of velocity, an analytic expression for velocity skewness (R) can be obtained from b and ∅ and is expressed as: Similarly, acceleration skewness (β) is obtained from r and ∅ by using the following equations: where β r_0 is written as: and where F 0 is written as: F 0 = 0.59 + 0.14(2r) (−6.2) , i f r ≤ 0.5 (A21) After finding velocity skewness (R) and acceleration asymmetry (β) parameters, the time of zero up-crossing and down-crossing of the waveform can be calculated from parameter c as: where c = bsin(∅). The dimensional wave period in crest (T c ) can be found from t zd and t zu (dimensionalized by dividing with angular frequency ω).
The dimensional wave period in acceleration during crest and trough cycle is calculated based on t zu and t zu (Figure 2b) as: where t m corresponds to the maximum or minimum phase in the expression defining time series of acceleration and can be written as: Next, the non-dimensional peak crest and trough velocity are computed as: The crest and trough velocity are dimensionalized using the amplitude of wave orbital velocity (û) and described as: By the end of the first step, one obtains a skewed waveform with the following parameters: velocity skewness (R) and acceleration asymmetry (β), accelerative/total wave period (T cu , T c , T tu and T t ), and peak velocity (u crest and u trough ) in each half cycle.

Appendix B. Calculation of Wave-Driven Bedload
Calculation of wave-driven bedload due to asymmetric waveforms was done using the Santos formulations [38]. After obtaining the skewed waveform, calculations are performed to get non-dimensionalized shear stress (Shields parameter) for the full wave cycle (timeaveraged absolute Shields parameter) and then the half wave cycle. For the full wave cycle, the ripple dimensions (height and length) [62] are determined from the maximum mobility number (φ) that is defined as:φ where u w is the amplitude of orbital velocity and the denominator term s = ρ s /ρ, where ρ s and ρ are the densities of sediment and water, respectively, g is the gravitational constant, and d 50 is median grain size. The denominator term (s − 1)gd 50 appears repeatedly in the following text and is only defined here for the sake of brevity. Above aφ of 240, the ripple regime transitions into a sheet flow regime. Once the bedform is determined and ripple dimensions are obtained, the time-averaged Shields parameter along with the current and wave friction factor for the full wave cycle are computed iteratively. The iterations need to be performed because the time-averaged Shields parameter depends on the friction factor calculation which in turn depends on the time-averaged Shields parameter. Note that for the current friction factor and near-bottom current velocity (u d ) calculation, the edge of the wave boundary layer needs to be determined [63]. Since only the wave-driven bedload is considered in this work, u d and corresponding current friction factors are equal to zero. Next, the effect of surface waves that modify the Lagrangian motion at the sediment bed by increasing the crest period and shortening the trough wave period is accounted for. The change in wave period is dependent on the wave propagation velocity (c w ) and the horizontal grain displacement during each half cycle ûT π and can be parameterized as: where ζ = 0.55 [64], ∆T c and ∆T t are the change in crest and trough period that is added and subtracted in crest and trough wave periods to obtain the new half-cycle periods (T c_new , T t_new ). Based on the new half-cycle periods, the accelerative components are also modified proportionally to get T cu_new and T tu_new and can be written as: After the wave period modification from the above two steps, the magnitude of crest and trough Shields parameter is computed. Each half cycle Shields parameter is then modified by the effects of progressive surface waves that account for the vertical gradient of horizontal momentum caused by the vertical orbital motion of surface waves. This phenomenon is parameterized through the contribution of a wave-averaged Reynolds stress term (τ wRe ). The effect of τ wRe is to increase the crest cycle stress (direction of wave propagation) and decrease the trough cycle stress. τ wRe is expressed as: where α w = 0.424, f wd is the full-cycle wave-current friction factor and c w is the wave speed. Following the calculation of Shields parameter in each half cycle, the next step is to obtain a phase lag parameter that determines the distribution of sediment load between each half cycle. The phase lag parameter is comprised of two terms that account for the ratio of sediment stirring height l bed f orm 2(T i −T iu )w si and effects of progressive surface waves that induce horizontal sediment flux gradients 1−εû i c w [58] in each half cycle. The phase lag parameter in each half cycle can be written as: where α and ε are set to 8.2 and 1.7,û i is the peak half cycle velocity, c w is the wave, l bed f orm corresponds to ripple height or sheet flow thickness, T i is the duration in a particular half cycle while T iu is the duration of accelerating flow in that cycle and w si is the sediment settling velocity. w si is calculated using the settling velocity for still water [39] and modified due to the vertical orbital velocity in each half-cycle caused by progressive surface waves [38,58]. If the phase lag parameter in each half cycle is less than one, then all the sediment entrained in each half cycle is transported within that half cycle. In that case, the bedload transport mechanism is only controlled by the Shields parameter calculated in that half cycle. If the phase lag parameter is greater than one, then the sediment entrained in one half cycle can be transported to the other half cycle. The resulting volumetric bedload transport per unit width ( → q b ) is given by: where s = ρ s /ρ, ρ s and ρ are the densities of sediment and water respectively, d 50 is median grain size, → θ is the non-dimensional bed stress, (Shields parameter), → θ c and → θ t are the Shields parameter in the crest and trough, respectively and T c and T t indicate crest and trough cycle wave periods, respectively. T cu and T tu refer to the wave period of accelerating flow within crest and trough cycles. Ω cc is the sediment bedload that is entrained during the wave crest and transported during the crest period, Ω tc is the sediment bedload that is entrained during the wave crest and transported through the trough period, Ω tt is the sediment bedload that is entrained during the wave trough and transported during the trough period, Ω ct is the sediment bedload that is entrained during the wave trough and transported during the crest period.