Large Eddy Simulation of Sediment Transport Over Rippled Beds

. Wave-induced boundary layer (BL) ﬂows over sandy rippled bottoms are studied using a numerical model that applies a one-way coupling of a “far-ﬁeld” inviscid ﬂow model to a “near-ﬁeld” large eddy simulation (LES) Navier– Stokes (NS) model. The incident inviscid velocity and pressure ﬁelds force the LES, in which near-ﬁeld, wave-induced, turbulent bottom BL ﬂows are simulated. A sediment suspension and transport model is embedded within the coupled ﬂow model. The numerical implementation of the various models has been reported elsewhere, where we showed that the LES was able to accurately simulate both mean ﬂow and turbulent statistics for oscillatory BL ﬂows over a ﬂat, rough bed. Here we show that the model accurately predicts the mean velocity ﬁelds and suspended sediment concentration for oscillatory ﬂows over full-scale vortex ripples. Tests show that surface roughness has a signiﬁcant effect on the results. Beyond increasing our insight into wave-induced oscillatory bottom BL physics, sophisticated coupled models of sediment transport such as that presented have the potential to make quantitative predictions of sediment transport and erosion/accretion around partly buried objects in the bottom, which


Introduction
Rippled seabeds frequently occur in coastal waters with sandy bottom, and the geometry of such ripples strongly affects wave-induced bottom boundary layer (BL) processes. For this reason, many studies have attempted to model flow and sediment transport over ripples, using methods ranging from discrete particle models, in which individual particles are represented (Calantoni et al., 2005), to simply adjusting the effective bottom roughness (Nielsen, 1992).
Here we simulate wave-induced flows over vortex ripples using a previously developed and validated hybrid hydrodynamic model in which a "near-field" Large Eddy Simulation (LES) model, solving Navier-Stokes (NS) equations, is forced by a "far-field" model, solving inviscid Euler equations (Grilli et al., 2009;Harris and Grilli, 2012). Additionally, in the present work, wave-induced sediment suspension and transport, forced by the computed flow, are simulated with a model combining a semi-empirical reference concentration along the seabed and a standard equation simulating sediment transport and accretion. While much simpler than a discrete particle model, the LES of three-dimensional (3-D) flows over a complex boundary still requires significant computational time. The hybrid modeling approach makes it possible to limit the 3-D-LES computational domain to that necessary and sufficient for simulating the salient physics in a given problem.
The LES near-field model used in the present work is an extension of that reported by Harris and Grilli (2012). While still based on a modification of the LES model of Cui and Street (2001), this model improves upon earlier work (e.g., Gilbert et al., 2007) by considering the turbulent bottom BL flow as the (potentially large) perturbation of an inviscid flow over the same domain. The perturbation scheme consists of first dividing the total pressure and velocity fields into inviscid and viscous parts and then in rewriting the governing NS equations for the perturbation fields only, assuming the inviscid flow is known from computations in the far-field model. This yields new forcing terms in the perturbation flow equations, which are function of inviscid flow fields representing the incident wave forcing (similar to,for example, Kim et al., 2005;Alessandrini, 2007). This (one-way) coupling approach makes it possible using a variety of fully realistic nonlinear and irregular wave forcings of the BL flow, besides the commonly used simple oscillatory or linear wave flows (see, for example, Dean and Dalrymple, 1991). Harris and Grilli (2012), for instance, simulated the nearshore transformation of far-field waves over a (possibly) complex bottom with a fully nonlinear potential flow (FNPF) formalism in the physical space. The latter is often referred to as a numerical wave tank (NWT), for which efficient and accurate boundary element models have been developed for two-(2-D) and three-dimensional (3-D) problems (Grilli and Subramanya, 1996;Grilli and Horrillo, 1997;Grilli et al., 2001Grilli et al., , 2003. Harris and Grilli (2012) validated the hybrid coupled model analytically for laminar wave-induced BL flows and experimentally for turbulent oscillatory bottom BL flows. However, due to the lack of accurate reference data, the model has not yet been applied to and validated for both more complex wave forcing and/or bottom geometry; similarly, the current model has not yet been used and validated for modeling sediment transport. These extensions and validations are the object of the present paper.
As the height and length of ripples that form in coastal waters on a sandy bed are dependent on the local wave environment, physically reproducing these vortex ripples with similar dimensions in a laboratory, generated by progressive water waves, would require a very large experimental setup. Since, to a first order, vortex ripples are forced by horizontal water oscillations over the seabed, most vortex ripple experiments have been performed in oscillatory water tunnels (e.g., Ribberink and Al-Salem, 1995), whose flow is simply forced by a piston motion at one end. This is a vastly simpler laboratory setup than using a large wave tank, which still captures much of the dynamics and salient physics of the wave-induced BL flow, including the shedding of vortices from ripple crests. Hence, despite their idealization, such results are used in this paper to validate our model of flow and sediment transport over vortex ripples (van der Werf et al., 2007). One important limitation in this type of setup is that it does not allow for creating and thus measuring of either the BL steady streaming (Longuet-Higgins, 1953) or Stokes drift due to nonlinear wave flow asymmetry, which are both higher-order nonlinear effects.
As indicated above, a variety of complex forcings of the 3-D flow in the smaller near-field domain can and have been simulated in the larger far-field domain of the hybrid model, from simple spatially homogeneous oscillatory flows to spatially variable flows induced by nonlinear waves shoaling over a complex bottom topography, such as occurs in nature (e.g., Harris and Grilli, 2012). Here, in order to validate simulations against laboratory experiments of flows and sediment transport over vortex ripples, we will only use the simplest forcing of an oscillatory uniform flow. Such a forcing can be analytically defined without the need to run the FNPF model, but, once validated, the model can be used to simulate and study the effects of much more complex and realistic wave forcing. Additionally, for modeling experiments of oscillatory BLs inside laboratory water tunnels, the assumption of periodic boundary conditions may be suitable for the near-field perturbation flow, as long as turbulence is sufficiently well resolved. This simplification permits simulations over even smaller domains of simple shape, which adequately represent flow conditions in much larger experimental setups. This makes it possible to use a more refined numerical grid for the same computational effort.
Note that completely independent developments of the LES model of Cui and Street (2001), other than those that led to the work of Gilbert et al. (2007) and to the present work on suspended sediment transport, were pursued by Street (2001, 2006) and, to study bedform evolution, by Chou and Fringer (2008, 2010. The latter authors extended the model to consider an evolving bed, and, by devoting sufficient computer power and time, they were able to directly simulate the formation of vortex ripples on a sandy bed rather than assuming an initial perturbed shape as will be done here. The results from Chou and Fringer (2010), however, were only subject to limited comparisons with experimental data and, in particular, no direct comparisons of velocity fields, suspended sediment concentrations, or sediment transport rates with observations were made.
In the present work, the 3-D LES will be used to gain physical insight into oscillatory bottom BL flows, assuming a realistic wave forcing (such as afforded by the present hybrid hydrodynamic model). However, such sophisticated models have the potential to make quantitative predictions of sediment transport and erosion/accretion around partly buried objects in the bottom, an important problem in scour around and burial of pipelines (e.g., Brørs, 1999;Liang and Cheng, 2005), cobblestones (e.g., Voropayev et al., 2003), short cylinders (e.g., Voropayev et al., 2007;Testik et al., 2005Testik et al., , 2006Cataño-Lopera and García, 2006;Bower et al., 2007;Hatton et al., 2007;Trembanis et al., 2007), and bottom sea mines (e.g., Inman and Jenkins, 2002;Guyonic et al., 2007). Moreover, most bed morphology models proposed to date have been limited to 2-D problems (e.g., Jensen and Fredsøe, 2001;Soulis, 2002), and only a few models have recently been proposed which attempt to simulate 3-D scour (e.g., Smith and Foster, 2005). Once fully validated, the coupled FNPF-3-D-LES model used here could be applied on a larger scale to these problems, while accounting for fully nonlinear and shoaling effects in incident wave fields, as opposed to only considering uniform oscillatory flows or linear waves, as in work published to date.
In the following sections, we briefly present the hybrid inviscid/3-D LES and embedded sediment suspension and transport model equations and numerical methods. Then, we present detailed validation of the model against laboratory experiments in a water tunnel, for sediment transport induced by an oscillatory uniform flow over vortex ripples. This validation includes the comparison of computations and measurements of velocity and suspended sediment concentration fields.
u ≈ u I u = u I + u P Figure 1. Schematic of hybrid model of vortex ripples and boundary layer flow. In most of the domain, the wave-induced flow can be (and is) considered inviscid (left), but above rippled beds, turbulent vortices are the dominant cause of momentum transfer (right), produced by the mostly oscillatory flow. A fully viscous/turbulent flow is modeled in this region using LES.

Measurements and models of vortex ripples
Oscillatory flows over vortex ripples were initially modeled by assuming an inviscid fluid (e.g., Davies, 1979). More recently, though still idealized, one-dimensional eddy viscosity models were used as a practical method of modeling suspended sediment concentration and fluxes over ripples (e.g., Davies and Thorne, 2005). Because the dynamics of such flows are dominated by the coherent eddies formed at the ripple crests, with stochastic turbulence being a secondary process, discrete-vortex models have met with some success (e.g., Hansen et al., 1994;Malarkey and Davies, 2002). Models based on Reynolds averaged NS equations have also been commonly used (RANS; e.g., Eidsvik, 2006); Chang and Scotti (2004), for instance, compared RANS techniques with LES for modeling flows over ripples. Direct simulations of NS equations (DNS) of flows over ripples have also been performed (e.g., Scandura et al., 2000;Blondeaux et al., 2004), but there are stringent limits on the flow Reynolds number that can realistically be computationally achieved. Vortex ripples are found in a range of dimensions, but are characterized by flow separation in the lee of each ripple crest (e.g., Fig. 1). Bagnold (1946) described these shapes and the flow above them. Such flow separation spawns recirculating eddies, which are ejected away or released from the ripple crests at flow reversal. Thus, every half waveperiod, the wave-induced oscillatory flow induces sheet vortices over each ripple crest, which dominate momentum and sediment transport in the BL. Ripple formation has now been extensively studied, both for the more commonly considered long-crested ripples (e.g., Blondeaux, 1990;O'Donoghue and Clubb, 2001;Testik et al., 2005;van der Werf et al., 2007) and for 3-D ripples that form when waves approach the coast at an angle (e.g., Roos and Blondeaux, 2001). In addition, the relationship between ripple geometry and oscillatory flow parameters has been well established (see, for example, Wiberg and Harris, 1994).
Modeling ripple formation with a 3-D-NS solver (such as an LES) would require substantial computational efforts (e.g., Chou and Fringer, 2010). However, since ripple geometry under periodic flows rather quickly becomes quasisteady, in order to study fundamental physical processes and validate numerical models for those, our focus can be limited to studying the velocity field, suspended sediment concentration, and sediment transport rates over a rippled bed of specified (albeit realistic) geometry.
Ripples have also been studied in a large variety of field experiments, and measured suspended sediment concentrations over rippled beds were compared with existing models of ripple characteristics (e.g., Vincent and Green, 1990;Green and Black, 1999;Grasmeijer and Kleinhans, 2004). Several experiments have looked at the evolution of sand ripples over time as wave conditions change (e.g., Hanes et al., 2001;Hay and Mudge, 2005), as well as considering the effective roughness or wave friction factor of the ripplecovered bed (Hay, 2008). Using a multi-instrument tripod, Traykovski et al. (1999) made detailed measurements of current and vertical profiles of suspended sediment concentration; using a sidescan sonar, they simultaneously measured the bedform geometry evolution over 6 weeks of observations, which included the passage of several tropical storms. Even more detailed field measurements and analyses of mean flow and turbulent statistics were conducted by Williams et al. (2003), including hydrodynamic conditions, bedforms, and suspended sediment concentration. Detailed particleimage velocimetry (PIV) measurements of flow fields in the coastal bottom BL have also been made (e.g., Nimmo Smith et al., 2002), but these are not currently as well suited for comparison with somewhat idealized numerical simulations as laboratory observations.
As indicated, laboratory experiments can provide more controlled conditions for studying and measuring flows over ripples, but few have measured both the entire flow and suspended sediment concentration fields, while also reproducing the same types of flow conditions as seen in the field. Ribberink and Al-Salem (1995) made detailed time-dependent measurements of flow velocity and suspended sediment concentration, but in sheet flow conditions over a flat bed. Faraci and Foti (2001) studied the evolution and migration of rolling grain ripples over a seabed, which are on a smaller scale than vortex ripples and are generated not as a result of the lee vortex that appears each half-cycle in vortex ripples but due to the motion of sediment along the seabed. Thorne et al. (2002) measured ripples in a large wave flume, including bedform morphology and suspended sediment concentration. These experiments were limited, though, as they only included a few flow measurements, and these were only obtained from electromagnetic current meters, which do not resolve the vortices that dominate the momentum transfer in the BL. Marin (2004) measured both the flow field and Eulerian drift over ripples under progressive waves, but at low Reynolds number and with a fixed bed. Furthermore, Rousseaux (2008) made PIV measurements to make a detailed study of the vortex dynamics above ripples, but did not focus on the suspended sediment concentration. It is less common to measure the flow field and the suspended sediment concentration, over a mobile bed, with full-scale ripples, such as van der Werf et al. (2007). This latter data set is used in our work for validating our LES computations.
Several modeling approaches for similar problems have been recently proposed, although not in a coupled/hybrid environment allowing for more complex wave forcing and seabed geometry to be studied. Van der Werf et al. (2008) modeled the same laboratory experiments as considered here (van der Werf et al., 2007) with both k-ω and discretevortex particle-tracking models. They showed reasonable agreement with measurements both in terms of the velocity field and sediment transport. Models of suspended sediment transport over ripples using a LES method similar to ours (although independently developed) were proposed by Zedler and Street (2006) and Chou and Fringer (2010). Using an LES, Chou and Fringer (2010) also simulated the longerterm evolution of ripples on the seabed, but with a less detailed comparison with experimental results than presented here.
In our approach, the total velocity and pressure fields are expressed as the sum of irrotational (thus kinematically inviscid) and near-field viscous perturbation flow components, above a rigid seabed of arbitrary geometry. The NS equations are formulated and solved for the perturbation fields only, which are forced by additional terms, as a function of the incident forcing flow fields. As compared to other stateof-the-art LES models, used in similar applications, the oneway coupled hybrid/perturbation approach used here is both more efficient and brings the ability of representing more realistic nonlinear incident wave fields. While only simple applications of the model will be presented here as part of its experimental validation, the present work serves as a test case for a method which could easily be adapted to address much more complicated scenarios than can be addressed with other models, such as sediment transport caused by irregular nonlinear waves around partially buried objects. This will be the subject of future work.

Large eddy simulation
The Navier-Stokes (NS) equations, assuming a Boussinesq approximation, for an incompressible, isothermal, Newtonian fluid, with a non-cohesive suspended sediment concentration (i.e., SSC), C, read where u i and p are the flow velocity and dynamic pressure, respectively, in a fluid and sediment mixture of density ρ, with ρ 0 the fluid (i.e., water) density and ν its kinematic viscosity, and ρ s the dry sediment density, with ρ = (1−C)ρ 0 +Cρ s and s = ρ s /ρ 0 the relative sediment density. We adopt the indicial tensor notation convention, with x 3 denoting a vertical distance measured from some reference point (usually the free surface) and δ ij the Kronecker delta: Similar to Zedler and Street (2001) and Gilbert et al. (2007), the SSC is governed by an advection-diffusion transport equation with a constant settling velocity w s : where κ = ν/σ denotes the sediment diffusivity, with σ the Schmidt number. Note that this formulation of the SSC equation assumes that the concentration is low enough to avoid particle-fluid and particle-particle interactions beyond a constant settling velocity, which implies that the sediment dynamics do not affect the fluid flow much. The validity of this assumption is discussed by Villaret and Davies (1995) and Elghobashi (1994). Elghobashi (1994) states that a sediment suspension can be considered as dilute if the volume fraction of sediment is C < 10 −3 , and that the physical coupling between the fluid and particles can be considered to be truly one way for C < 10 −6 . Using these criteria in the experiments considered here, only small regions directly next to the sand ripples would be considered to be a dense suspension, which may nevertheless cause some effects on the turbulence that are not included in our model. We note, however, that earlier simulations using an approach to sediment transport similar to ours have been successful in predicting the sediment transport in the same experiments (van der Werf et al., 2008). Following Harris and Grilli (2012), let us denote by (u I i , p I ) the velocity and pressure fields of the ocean wave forcing flow, which is considered to be inviscid outside of a thin BL near the seabed. Such a flow is well described by Euler equations: Let us then introduce a decomposition of the total viscous flow into the sum of the latter inviscid free-stream flow and a defect or perturbation flow, of velocity u P i and pressure p P : Replacing Eqs. (7) and (8) with Eqs. (1) and (2), and subtracting Eqs. (5) and (6), yields the governing equations for the perturbation fields as Here the perturbation is defined in a region encompassing the near-field bottom BL of interest, which defines the computational domain (Fig. 1).
Although formally different, for the range of problems studied here, these equations can be shown to be nearly equivalent to the forcing of the total flow with the inviscid wave dynamic pressure gradient proposed by Gilbert et al. (2007) (with the exception of the inclusion of density variations), expressed as There are two key advantages, however, to the current approach, as compared to this earlier work: (1) boundary conditions can be more clearly and accurately defined for the viscous perturbation (i.e., as vanishing or using a radiation condition away from the wall); and (2) only the inviscid velocity is needed in the NS forcing terms rather than the dynamic pressure gradient (which requires additional computations). By applying a spatial-average operator (overbar) to the governing equations, we obtain the momentum equation for the resolved perturbation as where τ ij = u i u j −u i u j is the subgrid-scale (SGS) stress tensor and χ j = u j C − u j C is the subgrid-scale suspended sediment flux. Note that SGS models typically only consider the deviatoric stress τ ij − τ kk /3, because the resolved turbulent pressure, p * , is different from the resolved hydrodynamic pressure, with p * /ρ = p/ρ + τ kk /3. These governing equations are discretized in 3-D as in Cui and Street (2001), i.e., using a finite-volume formulation with second-order accuracy in both space and time on a nonstaggered grid. Quadratic upstream interpolation (QUICK; Leonard, 1979) is used for convective terms. Second-order centered differences are used for the remaining terms. The convective terms are time integrated using the second-order Adams-Bashforth technique, and the diffusive terms with a second-order implicit Crank-Nicolson scheme. The Poisson equation for the pressure field is solved with a multigrid technique (Perng and Street, 1991). In order to use sufficiently fine discretization in the simulations, the LES was implemented for parallel computing using Fortran and the Message Passing Interface protocol, for use on large computer clusters.

Experimental post-processing
The van der Werf et al. (2007) experiment referred to as Mr5b63 is used for comparison with the LES results. This experiment was conducted in an oscillatory flow tunnel, starting with a flat bed made of sand with a median grain diameter of d 50 = 0.44 mm. The flow velocity far from the boundary, u ∞ , temporally periodic and asymmetric, is well described by where U 1 = 0.54 m s −1 and U 2 = 0.095 m s −1 , with a fundamental period of oscillation of T = 2π/ω = 5.0 s. Such a flow is aimed at simulating the near-bottom flow induced by a mildly nonlinear wave, where a negative velocity corresponds to an "offshore" flow and a positive one to an "onshore" flow. Under such forcing, the flat bed eventually evolved into a rippled bed, with a ripple wavelength of 0.41 m and height of 0.076 m, which stayed relatively steady (with a small migration rate of 18 mm min −1 ). Once the bed geometry reached a quasi-steady state, measurements were made of the velocity field (by means of PIV) and suspended sediment concentration field (with an acoustic backscatter system -ABS). Both the PIV and ABS measurements are statistical averages over several oscillations. The PIV measurements were phase-averaged over five oscillations. The ABS measurements were compiled while six ripples migrated past the instrument. The PIV measurements used the suspended sand as a seeding agent, which due to inertia and settling velocity does not exactly follow water particle trajectories. The sediment settling effect was attenuated in post-processing by forcing the velocity data to be horizontally periodic and  Table 1).
removing the horizontally averaged vertical velocity (which must be true from flow continuity). Note that, though van der Werf et al. (2007) did calibrate their data against other measurements of suspended sediment, the ABS concentration measurements are accurate only within a factor of 2, which limits the degree to which the suspended sediment transport rates can be expected to agree with LES results. Van der Werf et al. (2008) did, however, provide error estimates of the sediment transport measurements, and from that we can expect that the ABS results are not so inaccurate.
After measurements were made and the water had become still, high-resolution measurements of the ripple geometry were made with a laser displacement sensor. Six parallel profiles were measured, 40 mm apart, across the oscillatory tunnel width, with each profile measured every 5 mm, at a vertical resolution of 0.05 mm. In many of the early theoretical solutions or models of flow over ripples (e.g., Benjamin, 1959;Lyne, 1971;Longuet-Higgins, 1981;Tanaka, 1986), the ripple geometry was transformed into a flat bed through a conformal mapping. In such an approach, a complex series expansion such as used by Shum (1992) can provide a reasonably accurate representation of any measured ripple geometry: where z = x 1 + ix 3 ; ζ = ξ + iχ (with ξ and χ the transformed horizontal and vertical coordinates); N is the number of terms in the series; and k is the wavenumber, equal to 2π/λ, for a given ripple wavelength λ. The series coefficients (α n = a n + ib n ) needed to reproduce the measured ripple shape (Fig. 2) were computed by van der Werf et al. (2008) and will be used in the following computations (Table 1). The PIV measurements (Fig. 5) show the velocity structure through the typical period of oscillation, and suspended sediment concentrations measured by ABS (Fig. 7) show that sediment is being suspended by this flow. Also, at times when there are clearly high velocities on the leading edge of the ripple (e.g., at ωt = 60 • ; Fig. 5), and one would assume the  sand bed stress to be very high, local sediment concentrations are not particularly high relative to the rest of the ripple (although it is possible that it is, but limited to a thin layer that is not resolved in the observations). This has implications for the forcing and boundary conditions that are applied to the LES model, as described in the next section. Note that the PIV measurements presented in this section (Fig. 5) have not been corrected for the fall velocity of the sediment and are only presented for qualitative comparison. Also, note that the ripple shape measured by the laser displacement sensor does not perfectly correspond to the shape of the ripples based on the PIV and ABS measurements (this is particularly clear in Fig. 7). Some of the gaps between the ABS measurements and the measured ripple profile could be due to the high concentrations of suspended sediment in the bottom BL or acoustic reflections from the boundary that prevents measurements from being recorded.

Boundary conditions and forcing
The hybrid 3-D-LES model is used to compute the perturbations fields (u P i , p P ) over the near-field computational domain, corresponding to one ripple profile defined by Eq. (18) on the bottom (Fig. 2), of wavelength λ, which repeats itself when specifying lateral (streamwise) periodic boundary conditions. Based on the flow parameters, the orbital diameter of the motion far from the bed is quite close to one wavelength. Certainly, if the orbital motion were any greater, the domain size would be a major limitation to the present setup. Initial tests with larger domains of two or three ripple wavelengths, however, did not show significantly different behavior.
According to our hybrid modeling approach, the flow in this model is forced by specifying the inviscid velocity field u I i , such as defined in Eq. (7), which here should represent the inviscid part of the free stream flow used in the laboratory experiments (Eq. 16). As indicated before, because of the simple geometry and uniform flow considered here, rather than computing the velocity field u I i (x j , t) using an inviscid numerical model (such as a FNPF-NWT; see Harris and Grilli, 2012), one can analytically calculate it based on a conformal mapping that transforms the near-field domain in coordinates (x 1 , x 3 ) into a rectangle in coordinates (ξ, χ ). The latter coordinate system also defines the transformed LES model grid. This approach is very similar to that of Longuet-Higgins (1981), who modeled flows over ripples by a combination of inviscid flow (found by conformal mapping) and discrete vortices.
As in Harris and Grilli (2012), in order to increase the resolution of the numerical solution near the bottom, the LES grid is vertically stretched with a stretching ratio α = 1.1. Hence, the transformed grid, of dimensions (L 1 , L 3 ) and size (N 1 , N 3 ) in the vertical plane, is defined as for n 1 = 1, . . ., N 1 and n 3 = 1, . . ., N 3 (corresponding to the number of computational cells in the streamwise and vertical directions). We then use the conformal mapping to find the analytic expression of the inviscid velocity, based on the definition of the ripple shape in Eq. (18), as with u α (t) a slightly modified inviscid free stream velocity, related to the far-field (free-stream) velocity u ∞ (t) given by Eq. (16). Indeed, because both the latter velocity and the ripple shape are asymmetric in the laboratory water tunnel, if the ripples were in the open ocean, a non-zero Eulerian drift would be induced at the edge of the BL. In a closed water tunnel, however, a pressure gradient will form as a result to prevent any net water flux. Hence, this makes u α (t) slightly different from u ∞ (t). In order to include this effect without having to model the entire water tunnel in their simulations, van der Werf et al. (2008) forced the velocity at a certain height to match Eq. (16). A similar technique was used by Holmedal and Myrhaug (2006). Here we instead forced the average horizontal velocity to match u ∞ , similar to how the physical water tunnel is forced. This yields At the upper boundary of the LES computational domain, in contrast to the zero-gradient boundary conditions used by Harris and Grilli (2012), a free-slip boundary condition is specified, for which the normal (i.e., vertical) gradient of the horizontal velocity and the vertical velocity are both set to zero. As indicated before, in the free-stream x 1 direction, periodic boundary conditions are used for all the relevant fields in order to approximate an infinitely long oscillatory water tunnel. On the other lateral boundaries (the spanwise direction x 2 ), a no-slip condition is applied, similar to that induced by the side walls of the water tunnel. This has the effect of stimulating turbulence production in the initial flow oscillations computed by the model.
A log layer is specified along the bottom boundary, for which the local friction velocity u * is defined as where κ is the von Karman constant, taken to be 0.41, and u s is the locally resolved velocity in the direction tangent to the boundary (i.e., the resolved velocity vector at the grid point next to the boundary, with the normal velocity vector subtracted), and z 1 is the distance from the boundary to the center of the nearest grid cell. As discussed below, we test several different values for z 0 . Sediment motion at the seabed is governed by bedload transport, the settling of suspended sediment, and sediment pickup. These processes can be described by nondimensional parameters, including the density ratio, s, and the Shields parameter (i.e., a dimensionless bottom shear stress), θ , defined as with θ cr the critical Shields parameter and d 50,s the median suspended sediment grain diameter. Because of the grain size distribution, the median suspended grain size is smaller than the median grain size of all sediment within the water tunnel. Van der Werf et al. (2008) estimates d 50,s as 0.244 mm, which we use here. The onset of sediment motion on the seabed is defined by comparing the Shields parameter to its critical value. The latter is obtained from van Rijn (1993): where D * = d 50,s [(s − 1)g/ν 2 ] 1/3 , which gives a critical Shields parameter of 0.0314 for the present flow calculations. Note that this formulation neglects any effect of bed slope. Eventually, the simulations could be improved using (n) tracers corresponding to a variety of sediment size classes, and solving (n) advection-diffusion equations, instead of one class for the median suspended sediment diameter.
The bottom boundary condition for suspended sediment concentration is similar to that of the k-ω model of van der Werf et al. (2008); when the local instantaneous Shields parameter is below the critical value, zero sediment flux perpendicular to the ripple surface is assumed; at higher values the reference concentration relationship proposed by van Rijn (1984) for non-cohesive sediment with grain sizes between 0.2 and 2 mm is used: where T = (θ − θ cr )/θ cr is the transport stage parameter for a Shields parameter greater than the critical value, as suggested by Nielsen (1992). Note that this approach is different to using a sediment pickup function, as also suggested by Nielsen (1992). In preliminary tests, a sediment pickup function approach induced unrealistically large suspended sediment concentration values near the bed. Both the eddy viscosity and diffusivity are set to zero at the bed, and the surface stress is applied as similar to Harris and Grilli (2012) (with µ = ρν, the dynamic viscosity).

Subgrid-scale model
The governing equations for the LES contain subgrid-scale terms τ ij and χ j , which are modeled with the dynamic mixed model of Zang et al. (1994), based on the stress decomposition proposed by Germano (1986). Note that, as in Cui and Street (2001), the spatial gradient of the eddy viscosity is neglected in the discretized governing equations. A complete description of the SGS model, as applied to this numerical technique, can be found in Harris and Grilli (2012). As indicated before, owing to the assumed low SSC values, the effects of suspended sediment on turbulent fields are neglected in the LES model, and hence SSC is not explicitly included in the SGS closure scheme. In the present applications, we find that the SSC is indeed not often high enough to affect turbulence (i.e., above the 10 −6 limit given by Elghobashi, 1994). Hence, we are dealing with a dilute suspension, except when extremely close to the ripple surface. Additionally, for dilute suspensions with an SSC below 10 −3 , particle-particle interactions are negligible, so we consider the dynamic mixed model suitable. Finally, note that Chou and Fringer (2010) have argued that the effects of SSC on subgrid-scale physics are implicitly modeled in the LES model to some degree, through their effects on the resolved fields (via density fields, ρ(x i , t), in the Navier-Stokes equations).
As in Harris and Grilli (2012) and following Chow and Street (2004) and Chow et al. (2005), the eddy viscosity near the bottom boundary (wall) in the SGS model is increased in order to augment the near-wall shear stresses. By refining the resolution near the bottom boundary, we obtain computational cells with large aspect ratios, with a fine vertical resolution, but without resolving turbulence on these small scales, so the SGS model improperly predicts a very low eddy viscosity. Thus, under the assumption that near the bottom the flow can be approximated by a log-layer and that the eddy viscosity determined by the SGS model is negligible, the eddy viscosity is augmented as (ν T ) total = (ν T ) SGS + κu * zcos 2 π z 4 √ J /2z 1 (28) for z < 2 √ J /2z, with z being the distance from a point to the seabed, as before, z 1 being the distance of the center of the first grid cell to the boundary, and J being the Jacobian of the transformation used in deriving the discretized governing equations. This is an extension of the technique used by Harris and Grilli (2012) to curvilinear boundaries, and for a Cartesian grid, the near-wall thickness of 2 √ J /2z would reduce to 2 x 1 .

LES setup
The LES model described above is used to simulate the laboratory experiments of van der Werf et al. (2007). To do so, a modest grid size is used that has N 1 × N 2 × N 3 = 32×32×32 points, spanning a length of λ = 41 cm (one ripple wavelength), with on average a 50 cm height and 30 cm width. Similar to Harris and Grilli (2012), the simulation is run for 10 periods of flow oscillation T (i.e., 50 s), using a time step of 1.0 ms (i.e., 50 000 time steps).
Preliminary results showed significant differences for simulations with different surface roughnesses. For beds with fixed sediment, a surface roughness of z 0 = d 50 /12 is often assumed, but, as in Zedler and Street (2006) and Chou and Fringer (2010), larger roughnesses around z 0 = d 50 are expected for mobile beds because of grain saltation. Although there are some empirical relations relating the Shields parameter to surface roughness (see, for example, Camenen et al., 2006), for simplicity we considered fixed roughnesses. To show the sensitivity of the results on z 0 , we considered values of d 50 /12, d 50 /4, and d 50 . As instantaneous results will be found to be quite similar for various z 0 , unless mentioned otherwise results will be shown for z 0 = d 50 /4.
For processing results, we are interested in four types of averages: the phase-averaged results (i.e., the results averaged for a set of ωt values, separated by 2π ); the periodaveraged results (i.e., time-averaged over T ); the periodand ripple-averaged results (i.e., averaged results over ripple length λ at a given vertical height); and finally the cumulative average (i.e., period-, ripple-and vertically averaged; such as the total suspended sediment flux). For each of these types of results, we are interested in the velocity field, the SSC, and the sediment fluxes. For simplicity, we will denote the above-defined averages of, for example, q as q (ωt, x 1 , x 3 ), q (x 1 , x 3 ), q (x 3 ), and q , respectively. We then compute the horizontal averages by reinterpolating the results onto a uniform grid. For comparison with experimental data, we

Wall stress
The LES simulations predict that full 3-D turbulence quickly develops above the vortex ripple (Fig. 3). Accordingly, in order to compare with the essentially 2-D laboratory observations, spanwise averaging is applied to all of the results.
It is observed that the LES results quickly achieve a quasisteady periodic solution. This can most easily be seen in the spatially averaged wall stress (Fig. 4). Convergence is further demonstrated below, in terms of vertical profiles of horizontal velocity, as well as in the overall suspended sediment flux.
Note that when simulating 50 s of physical time and running the model on eight processors, the simulation takes approximately 16 h of clock time, or 128 CPU hours. This compares to the 45 120 CPU hours required for the bedform evolution simulations of Chou and Fringer (2010), although a direct comparison of computational efficiency is not possible, since their computations are for a more complicated physical scenario.

Velocity field
To compare LES and experimental results, the computed velocity field is plotted in Fig. 5 in a manner similar to the velocity vectors measured with PIV, i.e., for six phases separated by 60 • . Spanwise averaging was applied to the LES results, and the figures shows the last (10th) period of oscillation of the simulation. Comparing both, we see that computational results agree well with measurements. At a 0 • phase, when there is no flow in the far field, a strong offshore vortex occurs, although it is not as well formed in the LES results as in experiments. At 60 • , the flow is in the onshore direction, with lower velocities near the bed. At 120 • , there is a large lee (onshore) vortex. At 180 • , the flow in the far field is weak, but near the bed there is a moderate offshore flow. By 300 • , a clear lee (offshore) vortex has formed.
It is also useful to analyze results for the period-and spanwise-averaged velocity, which drives much of the sediment transport in the bottom BL. This is shown in Fig. 6, where we see that the LES results appear similar to the PIV measurements. The largest difference (right panel) is that the onshore (right-side) vortex is slightly different in the LES results. But, overall, the present LES model achieves a rather remarkable agreement with experiments, for the average velocity field, and one which is quite a bit better than that reported by van der Werf et al. (2008) in their Fig. 6. In their results, computed with a RANS (k-ω) model, they did not obtain a period-averaged vortex on the offshore side and they predicted a more symmetric period-averaged flow in their discrete-vortex model.

Suspended sediment
The phase-averaged suspended sediment field, plotted in Fig. 7 for six different phases during the 10th period of oscillation, shows good qualitative agreement with the ABS measurements. Primarily, there is a layer of very high SSC close to the ripple surface, which moves with the flow occurring above it. The major difference with the ABS measurements is an overprediction of SSC above the ripple crest. These results are also an improvement compared to the simulations of van der Werf et al. (2008). For example, their k-ω model significantly underpredicted how much sediment will be suspended above the ripple crests, and their discrete-vortex particle-tracking model shows suspended sediment clouds in different locations than in the observations near the times of maximum velocity (i.e., near 60 and 240 • ; see their  This may be because their discrete-vortex particle-tracking model only has sediment released from the bed at the crest.

Sediment flux
One of the main goals in simulating flow and sediment dynamics over vortex ripples is naturally to obtain accurate suspended sediment transport rates. If we integrate the total suspended sediment transport, q s = uC , we can then compare the LES results to the experimental data in terms of sediment fluxes. We ignore here the minor subgrid-scale effects in the post-processing, (i.e., the difference between uC and uC ).
To best understand total sediment fluxes, we can first compare the instantaneous observed (Fig. 8) and predicted ( Fig. 9) sediment fluxes. The largest difference appears to be driven by a high predicted suspended sediment flux onshore at ωt = 60 • . We can see from the previous results that this flux is mostly driven by overpredicting the SSC over the ripple crest, rather than differences in the velocity field.

Vertical profiles
In addition to considering instantaneous velocity and suspended sediment concentration, we can also consider the ripple-averaged vertical profiles of horizontal velocity,  suspended sediment, and suspended sediment flux (Fig. 10). For comparison, we show the results for all three surface roughnesses considered (i.e., z 0 = d 50 /12, d 50 /4, and d 50 ).
For z 0 = d 50 /12 we see that the velocity profile does not agree with the observations at all, instead showing large offshore velocities significantly above the ripple crest. At higher roughnesses, though, the velocity profiles are reasonable, and certainly within the range of results reported by van der Werf et al. (2008) in their Fig. 7, using different models. Note, for example, that their k-ω model appears to predict a timeaveraged water flow through the tank, which would not exist in an oscillatory water tunnel.
If we consider the vertical profile of suspended sediment concentration plotted in Fig. 10 (middle panel), the SSC profile is in good agreement for all cases (note, however, the logarithmic scale for this figure). As expected, larger surface roughnesses result in higher values of suspended sediment, but we see that z 0 = d 50 /4 shows the best agreement with observations. Notably the rate of decay with height of the SSC above the ripple crest is reasonable, in comparison to van der Werf et al. (2008), who underpredict the SSC above the ripple (e.g., for x 3 /λ > 0.4).
The computed suspended sediment flux profile shows the largest discrepancy with experimental results. We do see a maximum suspended sediment flux offshore just above the ripple crest, but there is a substantial onshore sediment flux that does not match observations, particularly for the z 0 = d 50 case. We can further compare the total suspended sediment transport, which was observed to be −10.6 ± 1.7 mm 2 s −1 . By averaging the results over the 6th-10th oscillations (from 25-50 s), we find that the z 0 = d 50 /4 predicts a suspended sediment transport of −2.80 mm 2 s −1 , which, among results for various roughnesses, yields the closest agreement with observation (Table 2). Table 2. Period-and ripple-averaged suspended sediment flux (as compared to the experimental result of −10.6 ± 1.7 mm 2 s −1 ) for varying surface roughnesses for each period of oscillation (i.e., the mean over 5 s) of the simulations.

Near-wall modeling issues
In view of various LES results obtained, we see that neither the fixed bed assumption of z 0 = d 50 /12 nor the z 0 = d 50 assumption for flow over ripples provides an accurate representation of the flow. Rather, selecting z 0 = d 50 /4 yields the best agreement of model results with experiments. As surface roughness varies with the Shields parameter, more complicated parameterizations may be necessary such as proposed by Camenen et al. (2006). This highlights the need for validating such sophisticated sediment transport models against a variety of experimental conditions. Note that the boundary conditions applied on the seabed, which has the underlying assumption that the surface stress can be predicted from a logarithmic velocity profile, are based on the premise that the flow is steady, when it is clearly not. More importantly, this condition does not take into account the effects of pressure gradients, which are extremely important for separated flows, such as seen here. Doing so would clearly require a more sophisticated wall model than a simple log-layer assumption. For instance, the modified log-layer assumption derived by Fourrière et al. (2007) could be applied, where both the local pressure gradient and the surface roughness are considered in deriving the mean velocity profile. A similar equation has been found by Loureiro et al. (2008) and Loureiro and Freire (2009) to be experimentally correct.
Additionally, the near-wall modeling is influenced not just by the actual boundary condition but also by the RANS-like near-wall eddy viscosity expressed by Eq. (28), which has been used previously by Zedler and Street (2006) and Harris and Grilli (2012). This transition between a smooth RANS solution to a well-resolved turbulent velocity field for an LES is actually a significant problem with hybrid RANS/LES schemes. This can be improved by using techniques such as controlled forcing or applying synthetic turbulence (see, for example, Keating et al., 2006). Actually, in the results presented here, the turbulent fluctuations above the ripple crest are mostly due to the lateral no-slip boundary conditions. This was verified in preliminary testing by using spanwise periodic boundary conditions and observing that no turbulent eddies occurred. While others have used initial turbulent conditions to trigger turbulence in similar simulations (e.g., Zedler and Street, 2006), an improved near-wall turbulence approach would provide a more general solution.

Conclusions
A new hybrid LES approach for modeling the Navier-Stokes equations was applied to the simulation of wave-induced sediment transport over sand ripples. This hybrid technique is likely to be particularly useful for modeling coastal flow processes occurring near the seafloor, under complex nonlinear incident wave forcing. In that case, one may only need to solve the full Navier-Stokes equations in a relatively small region above the seabed. Harris and Grilli (2012) have already shown this approach to be accurate for modeling turbulent oscillatory boundary layers over flat beds, and practical for coupling the LES model to numerical wave tanks. In this paper, we compared our simulation results for flow and sediment transport over vortex ripples with the experimental data of van der Werf et al. (2007).
We obtained good agreement for the velocity field, including the instantaneous velocity and the period-averaged velocity, as well as a reasonable agreement for the vertical profiles of period-and ripple-averaged horizontal velocity. We obtained a reasonable agreement of suspended sediment concentration, although the SSC above the ripple crest is higher than in observations; as a result, the overall suspended sediment flux is quite different from what is observed, although within the same order of magnitude. For a surface roughness of z 0 = d 50 /4 we predicted a suspended sediment transport rate of −2.80 mm 2 s −1 , as opposed to the observed −10.6 ± 1.7 mm 2 s −1 . This could possibly be improved with some minor changes to the model setup, particularly in the near-wall subgrid-scale model and surface roughness. As changing the surface roughness in the model substantially changed the total sediment flux, but did not appear to affect the velocity or suspended sediment concentration profiles as much, sediment flux may be a good indicator for validating LES of vortex ripples in the future, as opposed to just qualitative behavior or averaged velocity or suspended sediment concentration separately.
A similar modeling effort was reported by van der Werf et al. (2008), who, with their k-ω model, were able to obtain a suspended sediment transport rate only 26 % lower than that observed. This does not necessarily indicate that their model is quantitatively better, as their k-ω model appears to predict a time-averaged mass flow through the water tunnel, which is not realistic. While many of the results of such models, as well as our own hybrid LES, qualitatively agree with observations, there are substantial variations between models, as well as the changes in vertical profiles of velocity and suspended sediment flux due to different surface roughnesses. Considering all three models (the two models of van der Werf et al. (2008) and our hybrid LES) show significant variation even with a priori knowledge of the laboratory conditions that may not be known in the open ocean (e.g., the median suspended grain size, settling velocity), this highlights the need for more advanced models, such as the bedform evolution model of Chou and Fringer (2010), to be well validated before being used in general applications.
Future work may extend upon the present results, in particular, by improving the turbulence model used to produce better predictions of suspended sediment transport, and eventually include a moving seabed, allowing the shape of the ripples to evolve over the course of the simulation.