A Numerical Study of Sheet Flow Driven by Skewed-Asymmetric Shoaling Waves Using SedWaveFoam

: SedWaveFoam, an OpenFOAM-based two-phase model that concurrently resolves the free surface wave ﬁeld, and the bottom boundary layer is used to investigate sediment transport throughout the entire water column. The numerical model was validated with large-scale wave ﬂume data for sheet ﬂow driven by shoaling skewed-asymmetric waves with two different grain sizes. Newly obtained model results were combined with previous nonbreaking and near-breaking wave cases to develop parameterization methods for time-dependent bed shear stress and sediment transport rate under various sediment sizes and wave conditions. Gonzalez-Rodriguez and Madsen (GRM07) and quasi-steady approaches were compared for intra-wave bed shear stress. The results show that in strongly asymmetric ﬂows, considering the separated boundary layer development processes at each half wave-cycle (i.e., GRM07) is essential to accurately estimating bed shear stress and highlights the impact of phase-lag effects on sediment transport rates. The quasi-steady approach underpredicts ( ∼ 60%) sediment transport rates, especially for ﬁne grains under large velocity asymmetry. A modiﬁed phase-lag parameter, incorporating velocity asymmetry, sediment stirring, and settling processes is proposed to extend the Meyer-Peter and Mueller type power law formula. The extended formula accurately estimated the enhanced net onshore sediment transport rate observed under skewed-asymmetric wave conditions. progressive wave streaming, momentary bed failure, and touch-down of wave breaking-induced turbulence on the bed fundamentally impact sediment transport processes; and additional research on a range of grain sizes and wave conditions is recommended. An additional analysis using a variety of new SedWaveFoam cases in conjunction with existing laboratory data from oscillatory water tunnel and wave ﬂume experiments should be performed for a comprehensive parameterization, such as offshore sediment transport occurring under purely velocity-skewed waves.


Introduction
Sandbars are nearshore bathymetric features where sand forms submerged or partly exposed ridges. Sandbar evolution is known to be affected by the point of wave breaking [1], infragravity waves [2], and longshore currents [3]. Generally, energetic storms drive rapid offshore sandbar migration [4], whereas onshore-directed sandbar recovery is expected during quiescent periods [5]. Thornton et al. [6] showed that energetics-based sediment transport models, e.g., [7], were unsatisfactory for sandbar migration predictions, especially under mild wave conditions. Nearshore morphological evolution is closely connected to cross-shore sandbar migration, e.g., [8,9]. Critically, large-scale hydromorphological models (i.e., XBeach, [10]) require extensive calibration to correct for missing or inadequately resolved sediment transport mechanisms, e.g., [11,12]. Therefore, improved understanding of onshore and offshore sediment fluxes driven by waves and currents needs to be achieved to advance parameterizations for physical processes by including additional physics [12].
Nearshore sediment transport driven by progressive surface waves is influenced by numerous factors. Positive velocity skewness (i.e., larger velocity magnitude and shorter This study investigated the phase-lag effects on time-dependent sediment transport rates under skewed-asymmetric surface waves using the new numerical modeling strategy. SedWaveFoam was adopted to simulate sheet flow driven by shoaling waves on the sandbar crest with two different grain sizes and compared with the measurement data from the BARSED experiment [23,41,42]. SedWaveFoam explicitly resolves the complex and interconnected sediment transport mechanisms under realistic surface waves, allowing for the relative importance and temporal behavior of each process to be investigated. Brief model formulations are introduced in Section 2. Section 3 presents the experimental and model setups. Section 4 illustrates the general two-phase flow feature and comprehensive model validations. Section 5 discusses parameterizing the time-dependent bed shear stress and sediment transport rate, through combining previous SedWaveFoam results [43,46] to cover varying degrees of velocity skewness, asymmetry, and grain size. Lastly, the main conclusions of this study are summarized in Section 6.

Governing Equations
The capabilities of SedWaveFoam include resolving sediment transport processes based on a two-phase flow formulation, i.e., sedFoam [24,47]; free surface tracking while minimizing the numerical diffusion, i.e., interFoam [44,48]; and wave generation and absorption modeling at the inlet and outlet boundaries for the free surface waves, i.e., waves2Foam [45].
The mass conservation equations are based on the volume of fluid method [49]. Flow quantities are Reynolds-averaged to avoid resolving three-dimensional (3D) turbulence. In a multi-phase flow system, mass conservation equations for three phases considered here (air, water, and sediment) can be written as [50] ∂φ k ∂t where i = 1, 2 for two-dimensional (2D) flow, i.e., streamwise (i = 1) and vertical (i = 2) components. The variables u k and φ k denote the velocity and volumetric concentration of each phase where the superscript k denotes air (a), water (w), and sediment (s) phases. The summation of φ k imposes the global mass conservation equation as φ a + φ w + φ s = 1. In SedWaveFoam, air and water are considered to be two immiscible fluids following [44,48]. Thus, they are combined as an air-water mixture phase, i.e., φ f = φ a + φ w , where the superscript f stands for the air-water mixture fluid phase. The air-water interface, i.e., free surface, is numerically tracked by the interface tracking scheme [43,46] ∂φ w ∂t where u f = (φ a u a + φ w u w )/φ f and u r is the relative velocity between the air and water phases, obtained using the interface compression method [44,48]. The Reynolds-averaged momentum equations for the air-water mixture and miscible sediment phases are ∂ρ s φ s u s i ∂t + ∂ρ s φ s u s i u s where air-water mixture density, ρ f , is calculated as ρ f = (ρ a φ a + ρ w φ w )/φ f . Here, we specify ρ a = 1 kg/m 3 , ρ w = 1000 kg/m 3 , and ρ s = 2650 kg/m 3 for each phase, respectively.
The variable p f is fluid pressure, g = −9.8 m/s 2 is the gravitational acceleration, and δ i2 is the Dirac delta function. The surface tension effect on the free surface is calculated based on the surface tension coefficient (σ t = 0.074 kg/s 2 ), the surface curvature (γ), and the volume fraction gradient of the air phase. The fluid stress, τ f ij , is obtained using a two-equation k-ε model [43,46]. The interface momentum transfer consists of drag force and turbulent suspension [24,47] and satisfies Newton's third law, i.e., M f s The particle pressure, p s , and particle shear stress, τ s ij , are based on particle collision at low to moderate sediment concentration [51,52] and enduring contact at high sediment concentration [53,54]. More detailed formulations are provided in [24,43,46,47].

Experimental Setup
The BARSED experiment was conducted in the large wave flume (overall length of 104 m, width of 3.7 m, and depth of 4.6 m) at the O.H. Hinsdale Wave Research Laboratory at Oregon State University in 2015 [55]. The flume is equipped with a piston-type wave maker with active wave absorption. A fixed, barred beach profile was installed at roughly 1:3 scale, based on beach profiles observed during the Duck94 experiments in Duck, NC, USA [4,56,57]. A recessed, sand-filled sediment pit with dimensions 3.66 × 3.66 × 0.17 m (length × width × depth) was installed on the sandbar crest (see Figure 1). The still water depth was 2.45 m at the wave maker (h wm ) and 1.00 m at the sandbar crest (h pit ). Two different sand grain sizes were used during BARSED trials (specific gravity, s = 2.65): wave cases beginning with S1 used sand having d 50 = 0.17 mm (d 16   The free surface elevation, η, was measured along the profile with 11 twin-wire capacitance and 6 downward-facing ultrasonic wave gauges at 17 streamwise locations (100 Hz). The horizontal position and η of the wave maker face were also measured with a twin-wire capacitance gauge (100 Hz). Velocity and sediment concentration profiles were measured with an array of sensors installed at the center of the sediment pit. Near-bed velocity and sediment concentration profiles were measured at 0.001 m vertical resolution over a range of 0.03 m using acoustic Doppler profiling velocimeters, ADPVs (100 Hz), and conductivity concentration profilers, CCPs (8 Hz; [58]), respectively. The measurement region of the ADPVs and CCPs intersected the sediment bed such that ∼0.02 m of the 0.03 m vertical range was above the initial bed level. The velocity profile spanning the entire water column was measured with a stack of six vertically-separated (0.20 m vertical spacing) acoustic Doppler velocimeters, ADVs (100 Hz). The lowest ADV was positioned roughly 0.1 m above the initial bed level. Suspended sediment concentration profiles were measured from the sediment bed up to 0.47 m with a dual-probe array of 20 miniature fiber optic backscatter sensors, FOBS (8 Hz; [59]). Streamwise and vertical pore pressure gradients were measured with a T-shaped pore pressure transducer array, PPTA (100 Hz), buried 5∼10 mm beneath the initial bed level [23]. The PPTA consisted of four (three) horizontally (vertically) separated pore pressure sensors each separated by 0.02 m.
A broad range of wave period, wave height, and sediment size configurations were observed during BARSED, with the intention of measuring varying wave breaker types at the sandbar crest. Naming conventions in previous publications have identified each configuration by SaTbHc, where a is sediment size 1 or 2, b is the wave period in seconds, and c is the offshore wave height in centimeters. Individual trials were composed of 10∼40 monochromatic waves per trial and configurations with more than one trial were ensemble-averaged to quantify variability [55]. This study used S1T7H40 and S2T7H40 to quantify varying sediment responses to the same wave forcing ( Table 1). The waves generated in these trials shoaled over the constructed beach profile to skewed-asymmetric shapes with wave heights of H pit = 0.59 m at the sandbar crest. The transformed waves propagated further onshore and broke about 10 m shoreward of the center of sandbar crest ( Figure 1b). Thus, wave-breaking-induced turbulence was unlikely to affect the sediment transport processes on the sandbar crest under the prescribed wave forcing of the T7H40 wave cases.
The measured velocity and sediment concentration data were phase-averaged over individual waves combining 4 (1) trials for the S1T7H40 (S2T7H40) case. In S1T7H40, the individual waves with spurious noise from the sensors were removed, resulting in phase-averaged data over 67 individual waves. For S2T7H40, 10 individual waves were used for phase-averaging. Full experimental details and data analysis methods are described in [23,41,42].  [46]. ‡ Model results from [43].

Model Setup
A 2DV numerical flume was constructed as a geometric equivalent to the physical experiment [23,41,42]. The streamwise and vertical directions are denoted as x (positive shoreward) and z (positive upward), where the middle of sediment pit and the top of the initial sediment bed are defined as x = 0 and z = 0, respectively ( Figure 1). Flow in the spanwise direction (y) is assumed to be statistically homogeneous. Initially, the numerical flume was discretized with rectangular uniform grids of 16 mm width (dx) and 1 mm height (dz). The grids were deformed following the beach profile. In order to resolve sediment transport processes properly, the grids were scaled down to dx = 2 mm and dz = 1 mm near the sediment pit using three layers of unstructured triangular meshes (i.e., divided by 2 3 ). The length of the sediment pit was shortened by 1.1 m compared to the physical experiment to maintain rectangular refined grids for the region with high φ s . In total, 4.6 million grid points were used (see Figure 1a in [46]). The immobile beach profile and bottom boundaries were specified as a no-slip for u k components parallel to the wall, and a no-flux boundary for wall-normal u k components and scalar quantities (e.g., φ k and p f ). The top boundary was specified as an atmospheric open boundary calculating the inflow flux based on p f . Two spanwise lateral boundaries are treated as an empty condition in OpenFOAM.
Monochromatic waves were generated in the numerical model by imposing φ f and u f at the inlet boundary. Consistently with the physical experiment, a 50th order stream function was applied to obtain φ f and u f at the inlet using T = 7.0 s, h wm = 1.001 m, and H wm = 0.4 m where H wm is the wave height at the wavemaker. A relaxation zone slightly longer than one wave length was added on the offshore side of the numerical flume to assure numerical stability and to absorb reflected waves in the relaxation zone [45], so as to reproduce the active wave absorption in the physical experiment ( Figure 1). General wave propagation features are indistinguishable between the two modeled cases; hence, only the S1 case is shown in Figure 1 for brevity.

Model Validation
The spatial evolution of the sediment phase under the wave crest of the seventh wave is shown in the zoomed in x-z plane snapshots of volumetric sediment concentration, φ s , at the sediment pit for each sediment size (Figure 2a,b). Due to the net onshore sediment transport driven by the skewed-asymmetric waves, scour (accumulation) of φ s was observed on the offshore (onshore) side of the sediment pit in both S1T7H40 and S2T7H40 cases. However, the center of the sediment pit was not affected by the scour and accumulation of sediment observed at the edges. Horizontal wave-averaged sediment flux profiles, q s (x, z) = φ s (x, z, t)u s (x, z, t) , were computed in the span of −1 m < x < 1 m for each successive wave to iteratively identify a quasi-equilibrium in space and time. Notably, represents the wave-average operator. For both S1T7H40 and S2T7H40, I A > 0.99 was achieved for q s (x, z) under the sixth and seventh waves with a vertical coordinate shift; hence, by the seventh wave, a quasi-steady state assumption is valid. Metric, I A quantifies a similarity of trend and absolute accuracy between the two datasets, where I A = 0 and I A = 1.0 denote complete disagreement and perfect agreement, respectively. The near-bed temporal evolution of φ s up to the seventh wave at x = 0 is shown in Figure 2c,d.

Wave Height
Modeled free surface elevation, η(x, t), was compared with measurements at four selected cross-shore locations indicated with dotted lines in Figure 3a. It should be noted that the model results were obtained based on Reynolds-averaged governing equations (see Section 2). Thus, the modeled seventh wave was selected for model validations comparing the ensemble (or phase) averaged measured data, and for the further analysis, as was also done in [46]. The instantaneous water depth (h(x, t)) in the numerical experiments was calculated from the vertically integrated volumetric concentration of the water phase, h(x, t) = ∑ φ w (x, z, t)dz, and the free surface elevation was then calculated as η( is the still water depth at each location. Normalized root-mean-square error (NRMSE) was the metric used to evaluate the numerical model accuracy against BARSED measurements. NRMSE was computed via dividing the RMSE between the measured data and model results by the range of measured data (i.e., maximum value -minimum value). For instance, each η time series was normalized by the maximum measured H at each location. For the free surface elevation comparisons, the ramp up portion of the signal (t ≤ 30 s) was ignored to exclude errors that arose from different wave generation methods between the physical experiment and numerical model. The NRMSEs of η(x, t) validation for both cases are provided in Table 2. Overall, η(x, t) was well-captured by the model, with the NRMSE being less than 10% at the four sites discussed here.  Table 2. A summary of normalized root-mean-square-errors for selected variables in model validation.

Fluid Velocity
Velocity measurements at the center of the sediment pit from each trial in BARSED were separated into individual waves based on the timing of zero up-crossings in measured pressure time series near the sediment bed, and were subsequently aggregated into cases based on sediment size and prescribed wave forcing (e.g., S1T7H40). Two main components of fluid velocity (streamwise: u f ; vertical: w f ) are compared at three selected ADV locations ( Figure 4 for u f ; Figure 5 for w f ).The phase-averaged measured data are compared with the NRMSEs of Reynolds-averaged u f and w f from the model results in Table 2. In general, good agreement was obtained for the u f time series comparison with NRMSEs less than 4% for both cases at three selected vertical locations ( Figure 4). However, a slight phase difference existed between the measured and modeled data for w f , resulting in NRMSEs that were increased to 22.1% and 15.3% near the sediment bed for S1T7H40 and S2T7H40, respectively (Figure 5c,f). At higher vertical locations, NRMSEs were below 5% (Figure 5a,b,d,e).  The sediment bed elevation was actively changing for each wave during the BARSED experiment, due to non-zero net sediment transport. Thus, a local vertical coordinate system (z * ) was introduced at x = 0, referencing the sediment bed level at t/T = 0 [41,42]. Nearbed flow quantities, such as u f (z * , t/T), k f (z * , t/T), and φ s (z * , t/T), were based on the z * coordinate system for the phase-averaging. Following the physical experiment, the inflection point of modeled φ s at t/T = 0 based on a curve fitting method suggested by [37] was also defined as z * = 0. The free stream velocity, u f ∞ (t/T), was defined at z * top,max -that is, the maximum value of z * top (t/T) where z * top (t/T) is the vertical elevation of maximum absolute u f (z * , t/T) obtained at a given time. For each case, z * top,max were 28.6 mm (S1T7H40) and 36.5 mm (S2T7H40), respectively. Due to the larger grain roughness, the thickness of wave bottom boundary layer (WBBL) was larger in the S2T7H40 case, compared to that in the S1T7H40 case.
The model-data comparison of u f ∞ yielded NRMSEs of 2.5% ( Figure 6a) and 2.2% (Figure 6b) for S1T7H40 and S2T7H40, respectively ( Table 2). The degrees of velocity skewness (Sk u ) and velocity asymmetry (As u ) were quantified using u f ∞ . The definitions follow [41,42]: where H is the Hilbert transform. Modeled (measured) Sk u and As u were 0.56 (0.58) and 1.03 (1.15) for S1T7H40 and 0.56 (0.60) and 1.04 (1.16) for S2T7H40. In general, higher-order statistical numbers closely matched between the measured data and model results. Computed values for Sk u and As u suggest sediment transport was driven by velocity-skewed asymmetric waves in both S1T7H40 and S2T7H40. For reference, Sk u and As u for prior nonbreaking [43] and near-breaking [46] SedWaveFoam wave flume case studies are also provided in Table 1.
The vertical profiles of u f (z * , t/T) at the center of the sediment pit under the wave crest and trough were compared (Figure 6c-f). Velocities in the sheet flow layer were not measured with the ADPV due to high acoustic signal attenuation in the dense sediment layer. In the S2T7H40 case, ADPVs were positioned slightly higher; hence, near-bed velocities up to 9.3 mm could not be obtained. With an available dataset, the numerical model was able to predict u f (z * , t/T) well with NRMSEs (reference to u f ∞,max ) of 5.9% (Figure 6c, wave crest) and 6.7% (Figure 6d, wave trough) for S1T7H40, and 9.1% (Figure 6e, wave crest) and 7.5% (Figure 6f, wave trough) for S2T7H40 (Table 2).

Turbulent Kinetic Energy
The model-data comparison of near-bed turbulent kinetic energy (TKE) profiles under the wave crest and trough is shown in Figure 7. For each individual wave in the physical experiment, the high-pass Butterworth filter was applied to extract turbulent velocity fluctuations with a cutoff frequency of 1 Hz. It should be noted that the sensitivity of the chosen cutoff frequency for the TKE estimation was reported in previous studies [46,57]. To exclude the effects of evolving currents in the wave flume (e.g., undertow), the demeaned velocity,ũ f i , was used for the filtering, whereũ The measured TKE was calculated using the high-pass filtered demeaned velocity,ũ f : where i = 1, 2, 3 and overbar denotes the ensemble-average operator. Based on the two-equation k-ε turbulence model, k f was directed computed in SedWaveFoam [43,46].

Volumetric Sediment Concentration
Normalized volumetric sediment concentration profiles, φ s /φ s max (z * , t/T), are compared at the wave crest and trough (Figure 8). Scattering of measured φ s /φ s max is illustrated with the gray envelopes, showing ±1 standard deviation. The measured (phase-averaged) and modeled φ s max were 0.624 and 0.613, respectively. The numerical model predicts φ s /φ s max (z * , t/T) with NRMSEs of 5.1% (Figure 8c, wave crest) and 7.8% (Figure 8d, wave trough) for S1T7H40 and 4.1% (Figure 8e, wave crest) and 4.3% (Figure 8f, wave trough) for S2T7H40 (Table 2). Two potential uncertainties were noted in the observation data. First, the CCP sensors covering φ s /φ s max 0.03 tend to smooth out the profiles near the shoulder where φ s /φ s max 0.75 [61]. Secondly, background sediment concentration exists from the beginning of wavemaker motion in very dilute region (φ s /φ s max ∼ 0.01), or FOBS calibration uncertainties. Overall, the quality of comparison is similar to that reported in previous SedWaveFoam studies [43,46].

Horizontal Pressure Gradient
Pore pressure gradient measurements are presented here as 10-wave ensembles from individual trials (trial 74 for S1T7H40; see [23]) to ensure consistent PPTA burial depth. A third-order finite difference formula [62,63] was employed to compute the horizontal pressure gradient from the measured pore pressure at four cross-shore locations (x = −20, 0, 20, and 40 mm). The pore pressure transducers were buried 5∼10 mm below the initial sediment bed level. For the model results, a central difference method was applied to p f at z * = −6.8 mm (z * = −6.4 mm) for the S1T7H40 (S2T7H40) case, where the distance between each grid point was dx = 2 mm. Erosion processes for the S1T7H40 and S2T7H40 cases were slightly different due to different grain sizes; thus, the z * coordinate systems in the two cases were not identical. The horizontal pressure gradient is commonly nondimensionalized by the immersed weight of the particle; i.e., the Sleath parameter, S [64] is given as The agreement of S between the measured data and model results was reasonable (NRMSE = 10.5%) for the S1T7H40 case ( Figure 9b). However, the peak of S was 14.0% underpredicted around the flow reversal, presumably due to the slightly different erosion depth. In general, modeled S for S2T7H40 (red dash-dot curve in Figure 9b) was smoother compared to that for S1T7H40 (blue dashed curve in Figure 2b). Furthermore, S from S2T7H40 at t/T = 0.157 was 56.1% less than S from S1T7H40. Although z * were not identical for S1T7H40 and S2T7H40, this result is counter to the common assumptions used to derive S based on the boundary layer approximation when pore pressures are not directly measured (i.e., ∂p f /∂x ≈ −ρ f ∂u f /∂t). The almost identical u f ∞ from both cases (Figure 9a) should theoretically lead to equivalent S if the linear approximation is appropriate. The dissimilarity in the modeled S from two cases may suggest that a larger Shields parameter [65], θ (see Section 5.2), decreases the threshold S value to trigger the near-bed instability [46]. The maximum Shields parameters, θ max , for the S1T7H40 and S2T7H40 cases, were 1.68 and 1.34, respectively (Table 3). Figure 9. Time series of (a) modeled free stream velocity (blue dashed curve, S1T7H40; red dash-dot curve, S2T7H40) and (b) Sleath parameter, S, from measured data for S1T7H40 (gray solid curve, phase-averaged; gray envelope, ±1 standard deviation) and model results (blue dashed curve, S1T7H40 at z * = −6.8 mm; red dash-dot curve, S2T7H40 z * = −6.4 mm). and u f ∞ . The subscripts Q and G denote the quasi-steady approach and GRM07 method, respectively. By definition, t c = T 1 , t r = T 1 + T 2 , and t t = T 1 + T 2 + T 3 . † Model results from [46]. ‡ Model results from [43].

Parameterization
The following sections use model results to investigate key intra-wave features in the bed shear stress (Section 5.1) and phase-lag effect contributing to sediment transport rates (Section 5.2). Novel bed shear stress and sediment transport rate formulations are proposed using relatively accessible parameters (i.e., u f ∞ ) and provide improved accuracy when compared to the traditional quasi-steady approach. The goal of this study is to illustrate the unsteadiness in bed shear stress and sediment transport rate parameterizations, and thus alternate methods incorporating additional parameters, e.g., local acceleration [17,22,66], are not compared for brevity. The previously reported sheet flow results for nonbreaking [43] and near-breaking [46] waves are also presented in this study to provide a comprehensive parameterization method to ultimately predict the sediment transport rate applicable for various wave conditions. The velocity skewness (Sk u ) was similar for all cases (see Table 1). However, velocity asymmetry (A s ) for the S1T7H40 and S2T7H40 cases lay in between those from nonbreaking [43] and near-breaking [46] waves. The effect of grain size on phase-lag was also considered by comparing S1T7H40 and S2T7H40 cases that had different grain sizes under same wave forcing.

Bed Shear Stress
Bed shear stress, τ b , is a fundamental component of all sediment transport formulae and is critical to estimating sediment transport rates. It is essential to accurately quantify time-dependent τ b (t/T) to improve intra-wave sediment transport predictions. Timedependent τ b (t/T) is conventionally computed based on a quasi-steady approach [28,29,67]: where f w,Q is the wave friction factor and l τ,Q is the time-lag between τ b and u f ∞ , in which the subscript Q denotes the quasi-steady approach. The two parameters, f w,Q and l τ,Q , are empirical additions to the τ b,Q (t/T) formula with the potential to upscale to relevant error quantities in sediment transport rate calculations. In the quasi-steady approach, f w,Q and l τ,Q are calculated as constants throughout the wave cycle [68] where X Q = |u  Figure 10), ω = 2π/T is the wave angular frequency, and k s,Q is the equivalent Nikuradse sand roughness. In this study, k s,Q was adjusted to match the peak of parameterized τ b,Q at the wave crest with the directly modeled τ b,max for each wave condition (see Table 3). Alternatively, Gonzalez-Rodriguez and Madsen [33] suggested that the development of WBBL under asymmetric waves is different for the wave crest compared to the wave trough. Based on the Gonzalez-Rodriguez and Madsen (GRM07) method, f w,G and l τ,G are different for each phase of the wave, where the subscript G represents the GRM07 method. The duration between the zero up-crossing and u f ∞,c (zero down-crossing and u f ∞,t , where u f ∞,t represents the minimum free stream velocity), is defined as a quarter sinusoid for the wave crest (trough), denoted as T 1 (T 3 ) ( Figure 10). Likewise, the duration between u f ∞,c and zero down-crossing (u f ∞,t and zero up-crossing) is defined as T 2 (T 4 ). The summation of T 1 and T 2 (T 3 and T 4 ) is defined as T + = T 1 + T 2 (T − = T 3 + T 4 ), representing the duration of positive (negative) free stream velocity under the wave crest (trough). At the timings of maximum and minimum free stream velocity, i.e., at t = t c and t = t t , two wave friction factors ( f w,G (t c /T) and f w,G (t t /T)) and two time-lag parameters (l τ,G (t c /T) and l τ,G (t t /T)) were separately calculated using the equations similar to Equations (10) and (11) [33]: 60 ω 11 − 2.0 log 10 X G for 0.2 < X G < 10 3 0 for 10 3 < X G , (13) with It should be noted that k s,G is different from k s,Q due to different X G ; thus, k s,G was separately calculated to match the directly modeled τ b and parameterized τ b,G at t = t c ( Table 3). The time-dependent friction factor, f w,G (t/T), and the time-lag, l τ,G (t/T), were then obtained by a linear interpolation between t c and t t . Using the GRM07 method, time-dependent bed shear stress, τ b,G (t/T), was calculated as [33]: The inter-comparison of resulting f w and l τ between the methods based on the quasisteady approach and GRM07 is illustrated in Figure 11. The previously reported nonbreaking wave case (case ID: D-J&H, [43]) that has a similar magnitude of u f ∞,max (see Figure 11a-c) and T (see Table 3) is contrasted with the present cases. The major difference between D-J&H and other cases is in the velocity asymmetry. Of note, S1T7H60 [46] is excluded here for brevity. For the skewed-asymmetric waves, f w,G (t/T) between the wave crest and trough varied 23.1% for S1T7H40 (Figure 11d) and 26.1% for S2T7H40 (Figure 11e), respectively, whereas it only varied 4.6% for the nearly-symmetric nonbreaking waves (Figure 11f; also see Table 3). The variability of l τ,G (t/T) was more significant under the skewed-asymmetric waves. For instance, the quasi-steady approach yielded values for l τ,Q that are about three times larger than l τ,G (t c /T) (Figure 11g,h). In addition, l τ,G (t t /T) was estimated as zero due to longer T 3 via the GRM07 method (Table 3). On the contrary, l τ,G (t/T) fluctuates around l τ,Q = 0.26 for the D-J&H case; thus, it was not too deviated from the quasi-steady approach (Figure 11i).
Parameterized, time-dependent τ b computed with the different formulations for f w and l τ are compared with modeled bed shear stress (Figure 11j-l). The directly modeled τ b was obtained from the closures for particle collision [51,52] and enduring contact [53,54]. The parameterized τ b,Q from the quasi-steady approach was overestimated l τ under the skewed-asymmetric waves (Figure 11g,h); hence, an unrealistic phase lead was predicted in the parameterized τ b,Q , compared to the directly modeled τ b (Figure 11j,k). In addition, τ b,Q at t = t t were 29.2% (S1T7H40) and 20.1% (S2T7H40) overestimated using the quasi-steady approach due to the constant f w,Q (Figure 11d,e). On the other hand, GRM07 estimated τ b,G reasonably well. NRMSEs were calculated as 8.15% using the quasi-steady approach and 4.69% with the GRM07 method for the S1T7H40 case (Figure 11j). A similar trend was achieved in the S2T7H40 case where NRMSEs were calculated as 6.02% for the quasi-steady approach and 2.59% for the GRM07 method (Figure 11k). On the contrary, the improvement was less apparent for the symmetric nonbreaking waves (Figure 11l). The NRMSE between the directly modeled and parameterized shear stress was 4.79% using the conventional quasisteady approach, and 4.06% with the GRM07 method for D-J&H. The lower NRMSEs for GRM07 suggest the time-dependent f w,G (t/T) and l τ,G (t/T) improve the overall estimate for the time-dependent bed shear stress under asymmetric waves. Figure 11. Time series of (a-c) free stream velocity, (d-f) wave friction factor, (g-i) time-lag, and (j-l) bed shear stress (black solid curves, directly modeled; blue dashed curves, quasi-steady approach; red dash-dot curves, GRM07) for the present (a,d,g,j) S1T7H40 and (b,e,h,k) S2T7H40 cases, and a (c,f,i,l) previously reported D-J&H case [43]. Vertical dotted lines represent t = t c (circles, f w,G (t c /T) and l τ,G (t c /T)) and t = t t (squares, f w,G (t t /T) and l τ,G (t t /T)).
While the GRM07 method led to more accurate τ b estimates for asymmetric waves, relying on u f ∞ to compute τ b still has limitations. For example, a relatively sharp reduction of τ b after the flow peak (t/T ∼ 0.2) was due to the density stratification induced by the intense sediment suspension (see Figure 2) and cannot be captured using the any of the parameterization methods. The decrease of modeled τ b slightly after the wave crest was more apparent in the S1T7H40 case, compared with S2T7H40, because sediment suspension was more intense due to the smaller grain size (S1: d 50 = 0.17 mm; S2: d 50 = 0.27 mm).
It has commonly been believed that a significant increase of an onshore-directed sediment transport rate under the free surface waves is associated with the progressive wave streaming-induced net bed shear stress, τ pws ; e.g., see [17,66]. Contrarily, by utilizing the inter-comparison between SedWaveFoam and 1DV sedFoam (i.e., without the free surface), Kim et al. [43] showed that τ pws is relatively small, and the modulated wave-averaged velocity (i.e., mean current) was the major mechanism causing the increase of onshoredirected sediment transport rate. In this study, τ pws was calculated as 0.07, 0.23, −0.24, and 0.18 Pa for S1T7H40, S2T7H40, D-J&H, and S1T7H60, respectively, where the reduced τ pws in the S1T7H40 and S1T7H60 cases was influenced by momentary bed failure. Hence, the contribution of τ pws was excluded in this study for simplicity.

Sediment Transport
The validity of quasi-steadiness for sediment transport under skewed-asymmetric waves is verified in this section. In the sheet flow regime (θ 1), coarse sediments (e.g., 0.5 mm < d 50 < 2 mm [69]) are assumed to instantaneously respond to changes in u f ∞ [7,29]. On the other hand, for medium to fine sands (0.063 mm < d 50 < 0.5 mm, [69]), pick-up and settling processes may not be bounded in each half wave-cycle due to slower settling velocities [36][37][38]. In other words, sediment transport processes under the positive (T + ) and negative (T − ) phases (see Figure 10) may not be equivalent with large velocity asymmetry, particularly for fine sands.
The time-dependent total (i.e., the summation of bedload and suspended load) sediment transport rate, Q s (t/T), directly from the model, is computed as Ribberink et al. [70] suggested that the non-dimensionalized total sediment transport rate Φ s (t/T) can be parameterized as a function of the Shields parameter (θ) using a power law: where s = ρ s /ρ f is the specific gravity of the sediment (s = 2.65), t * = t − l τ (t/T) is the shifted time frame considering the time-lag between u f ∞ and τ b (see Section 5.1), and sgn(x) is the signum function (1 if x > 0, 0 if x = 0, −1 if x < 0). The Shields parameter, θ, represents non-dimensionalized bed shear stress (τ b ), and is defined as [65] The directly modeled τ b is employed for the analysis in this section. While the critical Shields parameter, θ cr , for incipient sediment motion, is defined as [71] θ cr = 0.14 d 50 sg where ν = 10 −6 m 2 /s is the kinematic viscosity of water. The empirical coefficients, m and n, in Equation (17), were m = 11 and n = 1.65 for the bedload-dominant sheet flow condition under oscillatory forcing in a flow tunnel [29,70]. The merit of using SedWaveFoam is that we can investigate the sediment transport driven by realistic free surface waves while resolving inter-connected mechanisms such as progressive wave streaming, wave shape streaming, and momentary bed failure. Kim et al. [43] demonstrated that Equation (17) can be extended for free surface waves. The enhanced onshore sediment transport induced by progressive wave streaming is incorporated using an alternate m value for each half wave-cycle (i.e., T + and T − ), and the non-dimensionalized total sediment transport rate is proposed as [43] Φ s where t r is the moment of zero down-crossing flow reversal (see Figure 10). Here, m + (m − ) represents the enhanced onshore (reduced offshore) sediment transport due to the effect of progressive wave streaming under the crest (trough). The subscript 1/2 denotes the half wave-cycle approach applied for the non-dimensionalized total sediment transport rate parameterization with m + (m − ) corresponding to the peak of positive T + (negative T − ) half wave-cycle [43]. The empirical parameters, m + and m − , are determined as where t * c = t c − l τ (t c /T) and t * t = t t − l τ (t t /T). Figure 12 shows the relationship between modeled Φ s and θ − θ cr at the wave crest (t = t c ) and trough (t = t c ) of all four cases. The estimated m + and m − values are provided in Table 4. The mean of m + and m − values out of all four cases adopted here was calculated asm = 21.53. The estimated m + and m − values do not significantly deviate fromm ( Figure 12); hence,m ∼ 21.53 is suggested when Φ s information is unavailable. In this study, m + and m − were assumed to be less relevant to the different velocity asymmetry and grain size. On the other hand, the relationship between Φ s (t/T) and θ(t * /T) during the T 1 (0 < t < t c ) and T 3 (t r < t < t t ) phases is shown in Figure 13, where the slopes represent the best-fit n values at each phase obtained using a linear regression method. The best-fit n i for each phase (corresponds to T i where i = 1, 2, 3, 4) is given in Table 4. In general, the best-fit n i values were significantly lower than the typical value of 1.65, particularly under the skewed-asymmetric waves. This implies that more intense sediment transport can occur with weak bed shear stress (i.e., small n i with θ < 1 predicts higher Φ s ) under skewed-asymmetric waves. Values in the parentheses represent R 2 coefficients of determination. † Model results from [46]. ‡ Model results from [43]. Figure 13. The relationship between non-dimensionalized total sediment transport rate, Φ s (t/T), and the phase-shifted Shields parameter, θ(t * /T), during the T 1 (blue circles) and T 3 (green squares) phases for (a) S1T7H40, (b) S2T7H40, (c) S1T7H60 [46], and (d) D-J&H [43]. Blue dashed and green dash-dot lines represent the best-fit n 1 and n 3 , respectively. Dotted lines denote n = 1.65.
Dohmen-Janssen et al. [36] introduced a phase-lag parameter, Pl, to incorporate unsteadiness effects. A semi-unsteady model was developed to include the phase-lag effect on the net, i.e., wave-averaged, sediment transport rate [36]. van der A et al. [38] expanded the phase-lag parameter based on the half wave-cycle concept [39] to account for suspended sediment exchange between each half wave-cycles. We further modified the phase-lag parameter to be, Pl i , separately defined for each quarter wave-cycle: in which i = 1, 2, 3, 4, where δ s,c (δ s,t ) is the modeled sheet flow layer thickness at the wave crest (trough), and w s is the particle settling velocity (following Stokes' law, w s = (s − 1)gd 2 50 /18ν). The modified Pl i represents whether the suspended sediment at the characteristic elevation (i.e., δ s,c and δ s,t ) from the previous half wave-cycle is able to settle prior to the beginning of the next quarter wave cycle.
The relationship between the best-fit n i for all four cases and the newly proposed phaselag parameter (Pl i ) defined in Equation (22) is shown in Figure 14. In general, n i decreases with increasing Pl i (Figure 14). Using a nonlinear least squares method, empirical formulae with three different types were tested. The goodness of fit was computed as R 2 = 0.79 for the rational type (gray solid curve in Figure 14), R 2 = 0.74 for the exponential type (gray dashed curve in Figure 14), and R 2 = 0.72 for the linear type (gray dash-dot curve in Figure 14). The rational type empirical formula yielding the highest R 2 was where i = 1, 2, 3, 4. The dimensionless, time-dependent sediment transport rate can be re-written in terms of m for the positive (m + ) and negative (m − ) phases of the wave cycle, and n i for each quarter wave-cycle, such that phase-lag effects are taken into account: where the subscript 1/4 denotes the quarter wave-cycle approach. The performance of the newly proposed quarter wave-cycle formula (Equation (24)) is compared with the half wavecycle approach Equation (20) in Figure 15. The primary difference between these methods is that Equation (20) utilizes a constant n = 1.65 without considering the phase-lag effect reducing n. The temporal evolution of Φ s 1/4 (t/T), calculated by Equation (24) incorporating four n i , agreed well with Φ s (t/T) directly obtained from SedWaveFoam results for all four cases (NRMSEs = 5.5%, 3.5%, 4.1% and 3.4% for the S1T7H40, S2T7H40, S1T7H60, and D-J&H cases, respectively). The resulting net sediment transport rates, Φ s 1/4 , were also close to Φ s directly obtained from the model results (see Table 4). On the other hand, the agreement of parameterized Φ 1/2 (t/T) with Equation (20) is poor (NRMSEs = 8.9%, 4.5%, 10.9% and 3.8% for the S1T7H40, S2T7H40, S1T7H60, and D-J&H cases, respectively), particularly around the flow peak when the phase-lag effect is amplified (Figure 15a,c). Furthermore, the parameterized Φ 1/2 was significantly underpredicted, compared to the model results (Table 4). Figure 15. Time series of directly modeled Φ s (t/T) (gray solid curves) compared to those parameterized by Equation (20) (blue solid curves) and Equation (24) (red dashed curves) for (a) S1T7H40, (b) S2T7H40, (c) S1T7H60 [46], and (d) D-J&H [43] cases.
In summary, the applicability of quasi-steady approach (i.e., n = 1.65) is limited for predicting the sediment transport driven by asymmetric waves. A significant underprediction of net onshore sediment transport rate roughly up to 60% was observed for the finer grain size cases under asymmetric waves (S1T7H40 and S1T7H60). The newly proposed method specifically accounts for phase-lag effects using varying n at each quarter wave-cycle, and substantially improves the prediction of both time-dependent and net sediment transport rates. Compared to the quasi-steady approach, the error in net sediment transport rate was maximal at 10.7% (minimum 2.9%) using the method proposed in this study. The analysis presented here represents a first step toward integrating enhanced onshore sediment transport mechanisms due to the phase-lag effect by modifying the quasisteady approach based on two-phase flow numerical model results and experimental data of [41,42,72]. For a complete parameterization and fully calibrated empirical coefficients (e.g., m + , m − , n i , Pl i ), subsequent studies may be required applying with additional wave conditions and grain sizes.

Conclusions
A numerical investigation of sheet flow under skewed-asymmetric shoaling waves was conducted using a free surface-resolving two-phase model, SedWaveFoam. SedWave-Foam simultaneously resolves sediment transport processes and flow fields in the entire water column beneath surface waves propagating over complex bathymetry. The integrated modeling approach can improve the understanding of wave-sediments interaction and the quality of predictions such as bed shear stress and sediment transport rate. In the present study, SedWaveFoam was further validated with two different grain sizes (S1T7H40 and S2T7H40) in a large-scale laboratory wave flume dataset [23,41,42]. SedWaveFoam accurately reproduced flow characteristics (i.e., free surface elevation, and streamwise and vertical velocities). Good agreement was obtained with WBBL processes, including the vertical profiles of streamwise velocity, TKE, and volumetric sediment concentration. The horizontal pore pressure gradient time series agreed well with the measured data for the S1T7H40 case.
Practical parameterizations for time-dependent bed shear stress and sediment transport rate were proposed using the model results. The bed shear stress parameterization method was compared to previously suggested methods based on the quasi-steady approach and GRM07 [33]. Generally, the quasi-steady approach overestimated bed shear stress in the wave trough but underestimated sediment transport rates, particularly for skewed-asymmetric waves. The GRM07 method more accurately characterizes the separated WBBL development and improves the bed shear stress prediction; however, bed shear stress reduction induced by density stratification cannot be captured by simply using the free stream velocity.
The phase-lag effect beneath each wave type was investigated to identify its contribution to sediment transport processes. For the symmetric nonbreaking waves, the Meyer-Peter and Mueller type power law sediment transport rate formula (e.g., [43]) yielded a good agreement in the prediction of time-dependent sediment transport rate and wave-averaged net sediment transport rate. However, the empirical coefficient n in the power law formula was shown to vary with the phase-lag effect under the skewedasymmetric waves. An empirical formula was proposed to calculate n i for each quarter of wave-cycle, as a function of extended phase-lag parameter (Pl i ) associated with the sheet flow layer thickness, settling velocity, and time taken to settle. A modified intra-wave approach to incorporate the phase-lag effect onto the time-dependent sediment transport rate prediction appears promising, particularly for fine grain sediments or flows exhibiting large velocity asymmetry. The newly proposed parameterization was able to capture enhanced net onshore sediment transport rate consistent with the directly modeled results, about 60% larger than that from the quasi-steady approach.
The analyses presented here considered three wave conditions (shoaling, near-breaking, nonbreaking) and three typical grain sizes (d 50 = 0.17, 0.24, 0.27 mm). We expect that progressive wave streaming, momentary bed failure, and touch-down of wave breakinginduced turbulence on the bed fundamentally impact sediment transport processes; and additional research on a range of grain sizes and wave conditions is recommended. An additional analysis using a variety of new SedWaveFoam cases in conjunction with existing laboratory data from oscillatory water tunnel and wave flume experiments should be performed for a comprehensive parameterization, such as offshore sediment transport occurring under purely velocity-skewed waves.