Parameterization of intrawave ripple-averaged sediment pickup above steep ripples

Near-bed sediment pickup is critical for predictions of intrawave suspension and in turn net sediment transport in coastal models. In the present study, numerical results from a two-dimensional Reyn-olds-averaged Navier-Stokes model are used to assess the functional relationship of intrawave ripple-averaged sediment pickup above steep ripples. The numerical model provides intrawave time histories of ripple-averaged near-bed velocities and turbulence, which are qualitatively interrogated to determine pickup functional relationships. Several speciﬁc sediment pickup formulations are implemented within the numerical model: expressions relating pickup to near-bed velocity or near-bed turbulent kinetic energy via the bed shear stress; and expressions relating pickup to near-bed shear production of turbulent kinetic energy. These are then tested via model-data comparisons of near-bed suspended sediment concentration. The results show that the traditional functions relating sediment pickup to near-bed velocity cannot lead to reasonable intrawave suspension predictions above vortex ripples under a ripple-averaged framework. Instead, relating sediment pickup to near-bed turbulence quantities, such as turbulent kinetic energy or shear production of turbulent kinetic energy, signiﬁcantly improves the numerical predictions for these conditions.


Introduction
Numerical coastal models cannot resolve all spatiotemporal scales. This well-known issue leads to the necessary parameterization of processes that occur at the unresolved small scales. For sediment transport, this is especially pertinent for near-bed processes, and how the unresolved small scales are represented is fundamental for the predictive ability of numerical models. In particular, accurate prescription of intrawave sediment pickup is critical for estimating sediment transport in nearshore environments. Hsu and Liu [2004] showed that prescribing concentration near-bed boundary condition based on actual measured values does lead to accurate predictions of intrawave suspended sediment concentrations, and they therefore pinpointed the parameterization of the near-bed sediment boundary condition as the key component responsible for poor model skill. Good predictive ability for intrawave variations is in turn essential for calculating net suspended sediment transport. Indeed, the intrawave correlation between velocity and concentration can generate or significantly modify net sediment transport. This has been observed for fine sand for which so-called phase-lag effects can be significant [e.g., Dohmen-Janssen et al., 2002;O'Donoghue and Wright, 2004].
Even though Hsu and Liu [2004] focused on flat beds, sediment pickup is a greater concern above rippled beds, which are ubiquitous in coastal environments. In particular, models that are able to resolve intraripple dynamical processes are typically applied (i) for reduced near-bed domains only covering a few ripples' wavelengths [Malarkey and Davies, 2002;Chang and Scotti, 2004;Zedler and Street, 2006;van der Werf et al., 2008] or (ii) for scaled down systems, in which surface wave propagation is also resolved [e.g., Dimas and Kolokythas, 2011;Huang et al., 2011]. Conversely, models applied at larger coastal scales (e.g., larger than the surface wave wavelength) usually cannot resolve intraripple dynamics [e.g., Amoudry et al., 2013] and therefore have to rely on ripple-averaged approaches. Nevertheless, such models can be useful in investigating coastal sediment processes [e.g., Hsu et al., 2006;Ma et al., 2014] and require appropriate parameterizations of near-bed physical processes above rippled beds, including intrawave near-bed sediment boundary condition.
Such boundary condition may prescribe either a concentration or a sediment flux value and further details on both approaches can be found in Amoudry and Souza [2011]. The so-called sediment pickup functions seek to provide information on the erosive flux of material from the bed to the water column. Classical sediment pickup functions depend on the bed shear stress [e.g., Amoudry and Souza, 2011]. However, the bed shear stress is typically obtained by adopting an adequate model. A common example assumes that the near-bed horizontal flow velocity follows a logarithmic profile and expresses the bed shear stress following the rough wall log law. Even though such model formulations for combined pickup and bed shear stress can lead to satisfactory results for the intrawave suspended concentration above plane beds [e.g., Amoudry et al., 2005], they typically cannot reproduce suspension peaks that occur at flow reversal [Nielsen et al., 2002], which have traditionally been related to convective processes. While they are only secondary over plane beds [e.g., Nielsen et al., 2002], these peaks can become the dominant entrainment mechanism over rippled beds when ripples are steep [e.g., Thorne et al., 2003] via the so-called vortex entrainment. In a recent study, Zhang et al. [2014] concluded that conventional one-dimensional vertical models, which are implicitly ripple-averaged above rippled beds, can perform well above flat beds but cannot reproduce behavior above rippled beds.
Although high-resolution modeling of sediment transport in seabed boundary layers has resulted in better understanding of the fundamental physical processes, upscaling it into larger-scale models is still somewhat scarce. For example, Vittori [2003] investigated the dependence of sediment suspension on flow quantities from results of direct numerical simulation of flow and particle dynamics. The best correlation coefficient for parameterization of pickup was found for a dependence on the shear production of turbulent kinetic energy, while a large hysteresis was observed for the traditional dependence on the bed shear stress. These findings have, however, not been implemented and tested in practical models that operate at a coarser resolution and do not resolve small-scale convective processes and vortex entrainment.
In the present study, we solely focus on an intrawave ripple-averaged framework, and all quantities are therefore implicitly ripple-averaged. We hypothesize that (ripple-averaged) pickup formulae that depend on flow variables other than the near-bed velocity (via the bed shear stress) perform significantly better at reproducing intrawave (ripple-averaged) near-bed suspension patterns above steep ripples. In particular, we will focus on the different terms appearing in the balance of turbulent kinetic energy. In section 2, we briefly present the physical processes responsible for entrainment of sediment in suspension, which will be parameterized. We introduce our modeling approach in section 3, including a description of the numerical model. We then use numerical results to assess, a priori, functional relationships for sediment pickup in section 4 before implementing and fully testing several formulations in section 5.

Vortex Entrainment Process and Parameterization
Seabed ripples can modify the characteristics of momentum transfer and sediment transport in the nearbed boundary layer. Such an effect on momentum transfer was alluded to by the measurements of Lofquist [1986] (also reported in Nielsen [1992, Figure 1.2.3]). Ripples are typically generated by the action of waves and/or tidal currents and their size is typically related to the flow energy. Following the classification of Davies and Thorne [2008] and Thorne et al. [2009], two-dimensional steep ripples occur for low energy flows. For more energetic flows, the ripple steepness decreases leading to two-dimensional and threedimensional transitional ripples. Finally for high energy flows, dynamically plane or plane beds are found. These different types of sediment beds result in different physical mechanisms dominating momentum transfer and sediment transport in the near-bed boundary layer.
Above dynamically plane and plane beds, near-bed momentum transfer and sediment transport are mainly caused by turbulent diffusion. Above steep ripples, coherent periodic vortex structures become important and result in vortex shedding and vortex entrainment. Numerous observational and numerical studies [see e.g., Davies and Thorne, 2008;van Rijn et al., 2013, and references therein] have contributed to a reasonably good understanding of the processes involved. During the onshore wave half-cycle, a vortex is generated at the lee side of the steep ripple. As the onshore flow reverses offshore, this vortex is transported upward and offshore past the ripple crest. A similar process happens during the other wave half-cycle, albeit typically with a weaker vortex due to velocity asymmetry. These vortices locally entrain and trap sediment, resulting in intrawave suspension occurring twice per wave cycle around flow reversal as the vortices are ejected, Journal of Geophysical Research: Oceans 10.1002/2015JC011185 which may impact significantly the net transport of suspended sediment [Bijker et al., 1976;Davies and Thorne, 2008].
Even though a significant number of intrawave and intraripple models have been successfully used to reproduce hydrodynamics and sediment transport above ripples, descriptions and models of vortex entrainment within an intrawave ripple-averaged framework remain relatively scarce in comparison. Nielsen [1988] followed a qualitative argumentation to propose sediment pickup as an explicit function of time. Davies and Thorne [2005] introduced a sediment pickup, again as an explicit function of time, following Stokes expansion of the overall problem. These approaches lead to satisfactory results because the sediment pickup is artificially forced to occur at the correct intrawave phase. However, it still remains unclear how such methods, prescribing pickup as an explicit function of time, may fit in with the increasingly popular numerical models which rely on time-integration of primitive equations for hydrodynamics, on Reynoldsaveraged Navier-Stokes turbulence modeling, and on suspended sediment mass balance for sediment transport. Instead, composite functions of time, following which the time-dependent pickup is a function of a time-dependent physical variable, should typically be preferred in such models. Without the introduction of an arbitrary phase, the pickup would then follow the following form: where / represents a time-dependent physical variable. An example of such a formulation is the aforementioned time-dependent pickup depending on the time-dependent bed shear stress (i.e., / in equation (1) represents the bed shear stress). It is important to note that using sediment pickup formulations either following Nielsen [1988] and Davies and Thorne [2005] or following equation (1) is not related to differences in the fundamental physical processes responsible for pickup of sediment from the bed. Instead, it corresponds to different parameterizations at different levels of complexity of the same physical processes.

Numerical Modeling Methods
The assessment of intrawave ripple-averaged sediment pickup formulations will be achieved via near-bed model-data comparisons. The data consist of measurements taken during experiments in the Deltaflume of Delft Hydraulics (now Deltares), Netherlands, for which a sediment bed 0.5 m thick and 30 m long was installed in the flume. Details of the experimental setup are provided in Williams et al. [2000] and Thorne et al. [2002]. The large scale of the flume enabled sediment process studies to be conducted at full scale under controlled conditions for weakly asymmetrical waves.
Electromagnetic current meters (ECM) fitted to the side of the flume measured the wave flow velocities at five locations above the sediment bed and sidewall-mounted pump samples provided wave-averaged suspended sediment measurements. Intrawave suspended sediment was measured by an acoustic backscatter system (ABS) deployed on an instrumented platform. Details of the analysis procedure are given in Thorne et al. [2002], Davies and Thorne [2005], and Thorne et al. [2009]. Ripple dimensions (wavelength and height) were measured using an Acoustic Ripple Profiler [Bell et al., 1998], which provided subcentimetric measurements of bed location, and are used here as input to the model. Over the course of the observational period, ripple dimensions remained approximately constant even though ripples tended to migrate up to a full wavelength [Thorne et al., 2002;Davies and Thorne, 2005]. For the few cases specifically investigated here, the ripples migrated by approximately one full wavelength, which allowed the suspended concentration field to be studied over a full ripple wavelength. In turn, the ABS data were processed to provide the ripple-averaged, intrawave, height-varying concentration field (for details, see Davies and Thorne [2005]) at high spatiotemporal resolution (1 cm vertical bins and 4 Hz frequency) The ECM data are used to validate the hydrodynamic model, thus confirming that it does reproduce adequately the wave propagation and the resulting mean flow patterns. We then use comparisons between the ripple-averaged ABS data and model results to assess sediment pickup formulations. Since direct measurements of sediment pickup are difficult, uncommon, and not available for the experiments considered, a key issue is how sediment pickup can be isolated based on the existing measurements and model results.
The suspension data cover a region up to approximately 1 m above the rippled bed. In general, intrawave suspension patterns will depend on both intrawave sediment pickup and intrawave mixing. However, the relative importance of pickup and mixing is not constant throughout the water column. In particular, the Journal of Geophysical Research: Oceans 10.1002/2015JC011185 contribution due to turbulent mixing ought to be minimum at the bed. This supports our practical choice to assess different sediment pickup formulations by comparing time histories of sediment concentration calculated in the first numerical grid above the bed with ABS-derived experimental values at the same elevation (taking the model bed location to correspond to the ripple crests).
The numerical model that we use is based on the COrnell BReaking wAves and Structure (COBRAS) model [e.g., Liu, 1998a, 1998b], which solves the two-dimensional Reynolds-averaged Navier-Stokes equations using a k2e turbulence scheme. A suspended sediment transport model has been added to it to reproduce experiments conducted in the Deltaflume [Amoudry et al., 2013]. Overall, this model provides temporal and spatial variations of velocity, sediment concentration, and turbulence quantities.

Double-Averaging Approach
Our modeling approach is implicitly based on ripple-averaged quantities. This means that a doubleaveraging technique is employed for the derivation of the governing equations, which results in additional covariance terms in the governing equations. Taking hÁi as the (turbulence) ensemble average with Á 0 as the related fluctuation, and Á as the ripple average withÁ as the related fluctuation, the additional covariance terms are due to the following: Replacing u j by c leads to a similar relationship for sediment fluxes. The first term on the right-hand side is the (double) mean flow contribution, the second term relates to the mean stress resulting from the intraripple motions, the third and fourth terms may be seen as a large-scale and small-scale components to the traditional Reynolds-averaged turbulent stress. Only the last three right-hand side terms require further closure. More detailed background on theoretical and mathematical concepts of double averaging may be found in Gimenez-Curto and Lera [1996], Nikora et al. [2007], and Rodriguez-Abudo and Foster [2014]. The large-scale turbulent stress (third term on the right-hand side in equation (2)) is typically modeled using a classical turbulence closure scheme (k2e here). The other two terms requiring closure are due to unresolved intraripple motions, and are briefly discussed now with respect to our model.
One point of view is to neglect these two terms both for momentum and sediment, which is a common simplification partly due to parameterization not being fully developed yet [Nikora et al., 2007]. However, recent experimental results indicate that these terms may be of similar magnitude to the typical turbulent stress in the near-bed region [Rodriguez-Abudo and Foster, 2014]. An alternative point of view would be that these additional stresses/fluxes can be modeled in a similar way to the traditional turbulent eddy viscosity/ diffusivity hypothesis used in RANS models. In theory, the turbulent viscosity/diffusivity hypothesis is used to model covariance terms arising following an ensemble average, and it may seem to be a significant leap of faith to use it to model covariances arising from spatial averaging, which would also require certain assumptions on the relationship between the different averaging techniques. However, these additional stresses/fluxes are due to the unresolved intraripple motions, and an interesting parallel is the representation of the effect of unresolved mesoscale eddies in ocean circulation models [e.g., Burchard et al., 2008], which is routinely done relying on eddy coefficients [e.g., Smagorinsky, 1963]. Since the additional stresses/ fluxes arising from the ripple average are also due to spatially unresolved motions, their contribution could then be introduced via an additional eddy coefficient.
Our current modeling approach follows the model of Amoudry et al. [2013], in which modified eddy viscosity and diffusivity were introduced to reproduce wave-averaged suspension above vortex ripples. We consider here that these near-bed modifications provide a crude but simple parameterization of the additional terms arising from the double averaging: the anisotropic part of hũ i ihũ j i 1hu i 0 u j 0 i1hũ i 0ũ j 0 i proportional to the double-averaged mean rate of strain tensor via the modified eddy viscosity (or diffusivity for sediment fluxes).

Wave Hydrodynamics and Suspended Sediment Model
Using classical tensor notation and now considering all quantities to be implicitly ripple-averaged, the governing equations for the (turbulence) ensemble-averaged flow field in the water, which is assumed to be incompressible, are Journal of Geophysical Research: Oceans where hu i i is the mean velocity in the ith direction, hpi the mean pressure, u 0 i the turbulent velocity, q is the fluid density, and g i is the ith component of the gravitational acceleration. The fluid density is taken to be a constant, which implies that the effect of sediment on buoyancy is neglected (justified by the small concentration observed). hs ij i is the viscous stress tensor.
The Reynolds stress tensor is modeled following the classical isotropic eddy viscosity hypothesis: where d ij is the Kronecker delta and hS ij i is the mean rate of strain tensor. The turbulent eddy viscosity, m t , is calculated from the turbulent kinetic energy k and the turbulence dissipation rate e and includes the nearbed modification implemented in Amoudry et al. [2013] m where C l 50:09, A 0 the orbital amplitude, x the angular frequency, k s 525g r 2 =k r is the roughness with g r being the ripple height and k r the ripple length, and m t0 50:004 following Nielsen [1992]. The multiplying coefficient used to calculate roughness from the ripple dimensions typically varies from around 10 to almost 30 [e.g., Amoudry and Souza, 2011], and we take it to be 25 here to follow previous work based on the same experimental data [Thorne et al., 2002;Davies and Thorne, 2005;Thorne et al., 2009].
The turbulent kinetic energy (TKE) k and the turbulence dissipation rate e are computed by solving their respective balance equations where m is the kinematic viscosity, P is the shear production, and in which the model constants take the well-known default values r k 51:0; r e 51:3; C 1e 51:44, and C 2e 51:92. Equation (7) is the modeled balance equation for turbulent kinetic energy, the terms of which are used in sections 4 and 5. In addition to shear production P and dissipation rate e, the second term on the left-hand side represents advective transport of TKE and the last term on the right-hand side diffusive transport of TKE. Since the effect of sediment on density is neglected due to the low concentrations considered, no buoyancy production appears in the TKE balance. Equation (8) is the balance equation for the turbulent dissipation rate e, which is taken to follow the classical k2e model.
The suspended sediment concentration, c, is computed by applying the conservation of mass of suspended sediment, which results in the following advection-diffusion-settling equation: where w s is the sediment-settling velocity, which is only accounted for in the vertical direction (x 3 or z) and which is taken to be constant both in space and time. K s 5m t =r c is the turbulent sediment diffusivity, and r c is the Schmidt number, which is typically less than one, and taken here to be r c 50:5. allowing for nonzero turbulent kinetic energy at flow reversal. The near-bed dissipation rate is then obtained following the law of the wall: where e b and k b are near-bed dissipation rate and TKE computed at the elevation z b above the bed. The bottom boundary condition for the sediment concentration is implemented via the sum of an erosive flux of sediment, E, from the sediment bed into the water column and of a settling flux for deposition on the bed. Deposition is due to gravitational settling; it is thus occurring for all values of bed shear stress and is computed from the sediment settling velocity and the near-bed sediment concentration. Specific practical formulations for E following equation (1) will be explicitly introduced in section 5.1.

Model Implementation
The two-dimensional model is implemented to reproduce the wave hydrodynamics observed in the Deltaflume. The model domain represents the flume in two dimensions (along-flume x and vertical z) and covers 250 m in the wave propagation direction and 6 m in the vertical direction. Waves are generated within the numerical domain using the internal wavemaker of Lin and Liu [1999] and sponge layers are employed at both lateral boundaries to minimize wave reflection. Since we focus on ripple-averaged approaches, the bottom boundary is taken to be a rough flat bed, the roughness of which (k s ) is determined from the measured ripple dimensions. A uniform grid size of Dx525 cm and Dz55 cm, which is sufficient to resolve the propagation of the wave, is used in the present implementation.
We will focus on numerical results for the following experimental conditions: the still water depth is 4.5 m, the wave height is h w 51:1 m, the wave period is 5 s, and the wave length is therefore approximately 30 m. The sediment bed consisted of medium sand of median diameter 329 lm, and the settling velocity of sand in suspension is taken to be w s 50:017 m/s following Thorne et al. [2002]. The observed ripples were longcrested and approximately two-dimensional with the following dimensions: g r 50:05 m and k r 50:42 m. The calculated bed roughness is k s 50:15 m.

Model Validation for Wave Flow Velocities
Prior to assessing the near-bed modeled suspension patterns, it remains important to confirm that the model can reproduce the wave hydrodynamics appropriately, so that the near-bed treatment (whether for near-bed hydrodynamics or near-bed concentration) can reasonably be isolated as the key issue toward good predictions. Figure 1 presents model-data comparisons for the horizontal flow velocity at five elevations above the bed. Values for the root mean square error of the velocity are reported in the respective panels. The peak wave velocities are almost always (under-)predicted well within 10% of the measured values. This comparison confirms the good representation of wave hydrodynamics by the model.

Intrawave Measured Near-Bed Concentration
The near-bed ABS measured sediment concentration is presented in Figure 2 and displays several important characteristics. As mentioned previously, all numerical and experimental near-bed quantities are taken to correspond to the first numerical grid, so are computed 2.5 cm above the model bed location. The modeled free-stream velocity refers to the velocity 30 cm above the bed. The peaks in concentration are clearly not linked to maximum wave flow (see Figure 2, bottom), even taking into account the phase lead of near-bed velocity and bed shear stress with respect to the free-stream velocity. Instead, the maxima in near-bed concentration occur around the flow reversals. It is important to note that the free-stream velocity is plotted in Figure 2: as such, the near-bed flow reversals will occur earlier in the wave cycle (see e.g., Figure 4). The second important observation is the significant asymmetry between the two half-cycle concentration peaks. Both of these features are conceptually consistent with the vortex entrainment mechanism and may be critical for net (wave-averaged) sediment transport.
Journal of Geophysical Research: Oceans 10.1002/2015JC011185 Nielsen [1988] stated that ''the pickup function must reflect the two puffs of sand into this area, brought by the travelling vortices'' in order to derive an explicit function of time. Following a similar argument, some important characteristics of the intrawave sediment pickup can be inferred from the measured near-bed concentration. In our ripple-averaged framework, a concentration peak, as can be observed in Figure 2, would be the consequence of a slightly earlier peak in sediment pickup. This implies that intrawave pickup variations would then exhibit (i) two peaks slightly before near-bed flow reversals, and (ii) asymmetry in the peaks' magnitude. In turn, various pickup formulations of the form of equation (1) can be assessed first  based on the desired intrawave temporal behavior, then on their ability to reproduce the correct near-bed concentration features, in particular timings of primary and secondary concentration peaks and level of asymmetry between the two peaks. Timing and asymmetry of peaks are defined following the notation in Figure 3.

Intrawave Temporal Behavior of Near-Bed Hydrodynamics and Near-Bed Turbulence
As mentioned previously, the observations provide a series of requirements for the intrawave temporal behavior of pickup, and of any physical variable on which it depends (i.e., /ðtÞ in equation (1)). These can be assessed via interrogation of the hydrodynamic numerical results. It is important to note that the bed shear stress is typically obtained from a model and, as such, a dependence of pickup on stress will actually result in functions of other hydrodynamic variables. The intrawave behavior of bed shear stress [e.g., Lofquist, 1986] is therefore not investigated per se in this section, instead we will focus on primitive hydrodynamic and turbulence variables. Key characteristics of the intrawave behavior of these variables are summarized in Table 1.   Figure 4 presents the intrawave behavior of the near-bed velocity components. The near-bed horizontal velocity shows the phase lead with respect to the free-stream velocity that is expected above a rough bed. Near-bed flow reversals (u b 5 0) occur approximately at t=T50:19 and t=T50:69, which corresponds to a phase lead of about 208. Following equation (1), the typical dependence on u b via the bed shear stress obviously cannot produce intrawave pickup with peaks around flow reversals. While other mathematical expressions for Eðu b Þ could results in flow reversal peaks (e.g., Eðu b Þ / 1= ðu b 1dÞ with d a small number to avoid singularity), little asymmetry would be produced. The near-bed vertical velocity (w b ) does present extrema around flow reversals, but again little asymmetry is observed between the two half-cycles ( Figure 4 and Table 1). Altogether, it seems improbable that a sediment pickup formulation dependent on either nearbed velocity component would result in appropriate intrawave suspension behavior.
The intrawave behavior of near-bed turbulence is presented in Figure 5 for turbulent kinetic energy, dissipation rate, and eddy viscosity. All three variables display similar behavior: within a wave cycle, near-bed turbulence increases as the onshore flow decelerates, reaches a first maximum around the near-bed flow reversal, and then decreases during the accelerating offshore flow. A similar behavior occurs during the second half-wave cycle but with less intensity. Overall, this results in asymmetric peaks around flow reversals for turbulent kinetic energy, dissipation rate, and eddy viscosity. Such behavior for TKE is qualitatively near-bed vertical velocity; k b : near-bed TKE; k b;t : time derivative of k b ; e b : nearbed turbulent dissipation; P b : near-bed shear production; D b : near-bed diffusive transport; A b : near-bed advective transport; m t;b : near-bed eddy viscosity. For velocities, the asymmetry is taken on the absolute values. consistent with recent observations [Hurther and Thorne, 2011], even though a full quantitative assessment is not possible for the conditions prescribed here as near-bed turbulence was not measured. The turbulence asymmetry between the two half-cycles reproduces known behavior for asymmetric waves that induce steady streaming [e.g., Trowbridge and Madsen, 1984;Scandura, 2007;Kranenburg et al., 2012], and is linked to differences in instantaneous flow profiles and near-bed gradients between the two half-cycles due to wave asymmetry [van der A et al., 2011]. Due to the nonlinear nature of wave boundary layers, even weakly asymmetric waves can induce such asymmetry.
Following Vittori [2003], we also present in Figure 6 the intrawave time histories for the different near-bed terms appearing in the balance equation of TKE (equation (7)): time derivative @ t k b , shear production P b , dissipation rate e b , advective transport A b , and diffusive transport D b . Table 1 summarizes the timings of both half-cycle peaks and their asymmetry for all near-bed variables presented in Figures 4-6. As expected from the behavior of k b ðtÞ (Figure 5b), the time derivative is alternatively positive and negative within a wave cycle. Similarly, the advective transport term also changes sign during a wave cycle. During the decelerating onshore stage for near-bed flow (20:07 < t=T < 0:19), production, dissipation, and transport all increase in magnitude at first; production peaks before near-bed flow reversal while e b ; D b , and A b peak at flow reversal. Overall, this results in @ t k b reaching its maximum about halfway through this stage. During the accelerating offshore stage for near-bed flow (0:19 < t=T < 0:43), production, dissipation, and transport terms all decrease, and result in a primarily negative TKE time rate of change (@ t k b ). Again, a similar and weaker behavior is observed during the second half-wave cycle.
Altogether, these numerical results enable an a priori assessment of potential functional relationships of the type of equation (1). Qualitatively, intrawave sediment pickup has to display two peaks around flow reversals with significant asymmetry. It is important to note that sediment pickup as defined previously should not change sign during a wave cycle, this makes a dependence on either @ t k b or A b unlikely at least from a practical point of view. More quantitative conditions can also be deduced from the timing and asymmetry of the ABS-derived C vol and the desired variable / should satisfy the following conditions: t 1 < 0:15; t 2 < 0:80, A > 1. Out of the remaining variables considered in this section, only a dependence on shear production would satisfy all constraints. All other variables exhibit peaks that are too late or without asymmetry. This a priori assessment is consistent with the findings of Vittori [2003], which concluded that the best option was for a pickup as a function of shear production of TKE.

Implementation of Sediment Pickup
Based on the previous numerical results showcasing the behavior of near-bed hydrodynamics and turbulence, we implement three specific pickup formulations. Two of these follow the conventional approach of specifying the pickup as a function of the bed shear stress, which is calculated in two different manners (as a function of u b or as a function of k b ). Even though our a priori analysis seem to contradict such a choice, it remains interesting as a baseline test. The third formulation implements pickup as a function of the turbulent kinetic energy shear production P b , as suggested by Vittori [2003] and the results of the previous section.

Bed Shear Stress Related Pickup
A common approach to determine intrawave sediment pickup has been to modify formulations obtained for steady flows. This approach results in expressions for which the pickup E depends on the timedependent bed shear stress: EðtÞ5E½hðtÞ, where h is the time-dependent Shields number (non-dimensional bed shear stress). In turn, the bed shear stress is also modeled. It is again common to assume near-bed logarithmic velocity profiles and to calculate the bed shear stress as follows: where u ? 5 ffiffi ffi s p b =q is the friction velocity related to the bed shear stress s b and j50:41 is the von Karman constant. u b is the horizontal velocity computed at elevation z b above the rough bed (in the first grid cell). k s is related to the ripple's dimensions (see section 3.2). The time-dependent pickup will thus be related to the time-dependent horizontal flow velocity at the bottom computational grid E u 5E½u b ðtÞ 2 . Other approaches have been proposed to calculate the bed shear stress in oscillatory flows from the near-bed turbulent kinetic energy (TKE) [e.g., Hagatun and Eidsvik, 1986;Amoudry et al., 2005], which results in E k 5 E½k b ðtÞ with k b the turbulent kinetic energy at the bottom computational grid.
The specific expressions implemented in the numerical model follow Garcia and Parker [1991] and relate pickup to a reference concentration, itself calculated following van Rijn [2007]. Overall, this results in the following set of equations for calculating the velocity-dependent pickup E u and the TKE-dependent pickup E k : where and where g r and k r are the ripple height and length. In both cases, the bed shear stress is enhanced because of variations between ripple crests and troughs (second term on the right-hand side). C 0 is obtained as follows: in which z ref 5 0.01 m in our case. In this set of equations, w s is the sediment settling velocity, d the median sediment particle size, s5q s =q f the specific gravity with q s the particle density, and q f the fluid density,

10.1002/2015JC011185
g the gravitational acceleration, m the fluid kinematic viscosity, h c the critical Shields parameter for incipient motion, which is determined from the bed sediment median diameter following van Rijn [1993], C l 50:09 one of the turbulence model constant.

Pickup Function of Shear Production of Turbulent Kinetic Energy
We introduce here an alternative functional relationship in the form of equation (1). Based on the results of Vittori [2003] and of section 4, the TKE shear production seems to be the best candidate. We therefore propose the following pickup function: where P b is the near-bed TKE shear production and n is a power. E P0 is a constant, the dimension and magnitude of which depend on n, and which is practically determined to result in appropriate wave-averaged suspension profiles. The simplest relation is linear for n 5 1, but the numerical results of Vittori [2003] showed a nonlinear dependence of pickup on P b , which is also tested via a quadratic dependence (i.e., n 5 2).

Numerical Results
Four pickup functions (E u , E k , E P with n 5 1, and E P with n 5 2) are tested and assessed based on comparison of near-bed suspended concentration (Figure 7). Data from the ABS and from two independent pump sampling systems are included for the waveaveraged profiles. The intrawave nearbed concentration comparison is presented with the mean value removed, i.e., forc b 5c b 2c b , where the overline represents averaging over the wave period. Table 2 presents the following metrics for the model-observation comparison of the near-bed intrawave concentration: mean absolute error forc b ; Figure 7. Model data comparison for near-bed concentration using different pickup functions for a wave height h w 51:1 m, and ripple dimensions g r 50:05 m, k r 50:42 m. Numerical results use either E u (black), E k (green), E P n 5 1 (blue), or E P n 5 2 (red). (a) Wave-averaged concentration profiles from numerical models (lines), ABS data (circles), and two pump sampling systems (squares and diamonds).
(b) Free-stream (solid) and near-bed (dashed) horizontal velocities. (c) Near-bed (wave-average removed) concentrationc b 5c b 2c b from models (lines) and ABS data (circles, with the shaded area representing the 95% uncertainty interval [Bendat and Piersol, 1986]). (square of the) correlation coefficient between modeled and measured near-bed concentrations; index of agreement following Willmott [1981].
All four pickup functions tested contain a pickup constant: w s C 0 for E u and E k ; E P0 for E P . In practice, the value of this constant determines wave-averaged suspension and thus influences the wave-averaged concentration profile and the value of c b . In contrast, the intrawave time-dependent behavior, in particular timing of peaks and peak asymmetry, results from the time dependency due to velocity, TKE, or shear production and it is relatively insensitive to the pickup constant. This constant therefore is a tunable parameter toward achieving reasonable predictions of wave-averaged suspension profile. Small differences between the four pickup formulations can be observed for the wave-averaged concentration profiles (Figure 7a). Even though these differences could be reduced by further fine-tuning the pickup constant, perfect agreement for c b has not been pursued here as it is not the focus of the present study and the intrawave behavior (i.e.,c b ) would not be significantly affected. Only E P0 is actively tuned to match the wave-averaged concentration profile, E u and E k remain untuned here.
The intrawave behavior for the near-bed concentration is a confirmation of the a priori assessment presented in section 4. Using E u as the pickup function leads to intrawave near-bed concentration which is severely incorrect both in terms of timing of peaks and peak asymmetry: peaks occur at maximum free-stream velocity instead of near-bed flow reversal, and little asymmetry is produced. This is in turn reflected in the poor statistical metric values in Table 2. An important reason for the failure is that the assumption of a logarithmic velocity profile (equation (11)) is clearly incorrect above steep ripples, where vortices are generated and shed. All three other formulations result in a shift of pickup and near-bed concentration peaks toward near-bed flow reversals. However, E k still results in concentration peaks, especially the first one, which occur too late within the wave cycle. E P formulations lead to further improvements in terms of timing of peaks, as shown by the significantly increased values for r 2 , and the best results overall. Even though all turbulence-related pickups produce asymmetry in the near-bed concentration peaks, the amplitude of the intrawave variations is still underpredicted in all cases.

Discussion on Generalizing the Approach
We test the same four intrawave pickup approaches (E u , E k , E P with n 5 1, and E P with n 5 2) further by implementing the model for a different set of wave and ripple conditions. The experimental conditions for   Figure 8 the model-data comparison for the near-bed concentration using the different pickup formulations E u , E k , and E P (n 5 1 and n 5 2) as introduced in the previous section. Once again, data from the ABS and two independent pump sampling systems are presented for the wave-averaged profiles. For the E p formulations, only the pickup constants (E P0 ) have again been tuned to match the wave-averaged concentration profile measured by the ABS, and the values are different for the two cases presented here.
There are important differences between the two cases presented in Figures 7 and 8. The second case, for larger waves, exhibits larger discrepancies between the different measured wave-averaged concentration profiles (Figure 8a). The measured time-varying near-bed concentration (Figure 8c) is also clearly more complex. There still is a strong concentration peak (t=T % 0:25) occurring around flow reversal and thus evocative of the vortex entrainment process. There seems, however, to be some level of discrepancy in terms of the timing of this peak between the two wave conditions (i.e., peak at t=T50:15 versus at t=T % 0:25). Two other peaks are also present in Figure 8c (t=T % 0:4 and t=T % 0:6)) and they do not seem to be related to vortex entrainment. Their provenance is less evident but probably results from a combination of pickup at maximum flow and advection [e.g., Nakato et al., 1977]. Even though the peak at t=T % 0:4 does occur around maximum offshore near-bed flow, it does not correlate well with the numerical results using E u .
The wave-averaged concentration profile obtained using E u is significantly lower than the other three pickup formulations, which we attribute to a lower prediction of stress following equation (13) in this specific case. Nevertheless, the concentration profile still matches well with some of the data from pump samples, we therefore consider it is acceptable and it does not significantly affectc b . The intrawave behavior of the numerical results is remarkably similar in both cases tested, in particular with respect to the timing of concentration peaks. A quantitative assessment of the model-data comparison is provided in Table 3. This additional case confirms again the inadequacy of the E u pickup formulation and using alternative pickup formulations (E k and E P ) leads to significant improvements. However, definite conclusion on the relative performance of the turbulence-related pickups (E k and E P ) remains difficult primarily because of the complexity of this case, for which vortex entrainment does not seem to be the only process responsible for concentration peaks. This is further compounded by the discrepancy and scatter in the observed timing of the main concentration peaks between different cases, as previously observed by Sleath and Wallbridge [2002]. Other than experimental uncertainty, this apparent lack of independence of the peak timing on wave conditions may be the result of a chaotic and/or nonlinear nature of the vortex entrainment process, which may complicate its parameterization.
Within the context of generalizing the intrawave pickup approaches discussed so far, it is finally important to note that vortex entrainment is only one of the processes responsible, as may be inferred from the measured near-bed concentration in Figure 8. Another example is above flat beds, for which suspension peaks related to flow reversals are only secondary [e.g., Nielsen et al., 2002], and for which the main suspension events can be reasonably well reproduced by traditional pickup functions. More complete parameterization of sediment pickup would thus need to account for additional processes, and could combine dependencies on TKE production and near-bed velocity as a function of the bed characteristics.

Conclusion
Both the a priori assessment in section 4 and the two model-observation comparisons in sections 5 and 6 indicate that relating sediment pickup to the near-bed horizontal velocity (via the bed shear stress for example) does not lead to reasonable results above vortex ripples in the ripple-averaged framework. While we focused on a two-dimensional model, similar conclusions have been found using one-dimensional vertical models [e.g., Zhang et al., 2014]. It is clear that using pickup functions of the type of E u would always lead to maximum pickup at maximum near-bed flow, and near-bed sediment concentration incorrectly peaking near maximum flow with little to no asymmetry between the two half-cycle peaks.
Relating pickup to near-bed turbulent kinetic energy improves the numerical results in terms of peak timing and asymmetry. This is due to the modeled near-bed TKE peaking near flow reversals, which is qualitatively consistent with observations above rippled beds [e.g., Hurther and Thorne, 2011]. We also considered each term in the turbulent kinetic energy balance and our a priori results suggest that the best approach is to parameterize sediment pickup as a function of TKE production confirming the recommendation of Vittori [2003]. Overall, the qualitative behavior of the modeled production-dependent pickup is what can be expected based on small-scale observations of the physical processes: entrainment occurs close to flow reversal and there is a strong asymmetry between onshore-to-offshore and offshore-to-onshore reversals. Further work is, however, still necessary in the light of our attempt at applying the new pickup formulation E P to different conditions. From an experimental and observational point of view, more data would be required to help determine the timing of peaks, and to provide information on the near-bed turbulence to inform the parameterization and validation of double-averaged turbulent stresses/fluxes. From a numerical modeling point of view, further work toward determining the form of E P is still necessary, especially concerning the values for both E P0 and n. We believe that this would probably be best achieved from analysis of small-scale numerical results (i.e., from Large-Eddy or Direct Numerical Simulations), instead of trying to directly match ripple-averaged RANS results to experimental data.