LARGE-EDDY SIMULATION OF OSCILLATORY FLOW AND MORPHODYNAMICS OVER RIPPLES

The objective of the present study was to study the morphodynamical development of ripples in a movable bed. The methodology is based on the coupling of fluid flow, sediment transport and morphodynamics. A well-resolved largeeddy simulation (LES) is employed for the simulation of the three-dimensional turbulent oscillatory flow and the corresponding bed and suspended sediment transport over a rippled bed. The evolution of the bed form is obtained by the numerical solution of the Exner equation based on the spanwise-mean flow and sediment transport conditions. The Immersed Boundary method is implemented for the imposition of fluid and sediment boundary conditions on the moving bed surface. Results are presented for ripple creation and propagation from a quasi-flat bed, as well as results of initially sinusoidal ripples adapting to water conditions, based on the mobility number, ψ. The numerical model demonstrates phenomena of ripple creation, propagation and migration, resulting in ripple lengths in agreement with those predicted by empirical equations. It was shown that under the same hydrodynamic forcing, the bed tends to reach the same equilibrium state, regardless of the initial bed form.


INTRODUCTION
Interaction of wave-induced oscillatory flows with the sandy bottom leads to the modification of the bed surface and the generation of coherent small-scale bed forms, generally known as ripples.The geometry of these structures has a strong impact on the wave-induced bottom boundary layer processes, which control sediment transport in coastal areas.Thus, the accurate prediction of the shape of ripples under certain flow conditions is an important issue in the coastal environment.
In oscillatory flow over ripples, coherent vortices are generated at the lee side of each ripple crest during each half-cycle and ejected into the outer flow at flow reversal.According to this vortex formation -ejection mechanism, Bagnold (1946) characterized these bed forms as "vortex ripples".Clifton and Dingler (1984) separated vortex ripples in three classes.Orbital ripples whose length, L r depends linearly on the wave orbital amplitude, a o , suborbital ripples whose length depends both on a o and the grain diameter, D g , and anorbital ripples whose length depends only on D g .Vortex ripples have steepness h r /L r > 0.1 where h r is the ripple height.
As far as ripple creation and evolution is concerned, significant progress has been made in understanding bed form dynamics.Blondeaux (1990) formulated a theory for the formation of sand ripples under sea waves, based on linear stability analysis, and successfully validated the predictions in comparison to experimental data.Andersen et al. (2001) introduced a one-dimensional model in order to describe the dynamics of ripple patterns under oscillatory flow.They focused on realistic conditions for the mass exchange between adjacent ripples, and they predicted the existence of a finite band of stable ripple wave numbers.Giri and Shimizu (2006) presented a depth-resolving, two-dimensional, morhodynamic model for free surface flows over bed forms.They used a nonlinear k-ε turbulence closure for reproducing flow and vorticity field.They noted that the use of an exponentially stretched grid in vertical direction improves the simulation results and they successfully simulated some important physical features of bed form evolution, such as merging and asymmetric geometry.Marieu et al. (2008) developed a morphodynamic model, based on a modified non-oscillatory central scheme for ripple growth and slope avalanching, in order to perform simulations of sandy bed evolution.They examined ripple growth from a quasi-flat bed, as well as from random initial bed forms, showing that the bed tends to reach the same equilibrium shape regardless of the initial bed shape.
In the present work, a morphodynamic model is being presented to study ripple creation and growth, as well as ripples adapting to changing water conditions.An LES method is used for the fluid flow and sediment transport simulation.The Immersed Boundary (IB) method was implemented for the imposition of fluid and sediment boundary conditions on the bed surface.The bed morphology evolution, driven by both suspended and bed load transport rate, was simulated using the spanwisemean flow conditions.A brief description of the methodology is presented in the following paragraphs.

Hydrodynamics
In LES, the large flow scales are solved directly, while the subgrid scales are modeled using subgrid-scale (SGS) stresses.The governing non-dimensional flow equations are the continuity and the Navier-Stokes equations where x i are the three Cartesian coordinates (x, y, z), u i are the corresponding velocity components (u, v, w), t is the time, p is the dynamic pressure, τ ij are the SGS stresses, Re is the Reynolds number, and f i is the source term associated with the implementation of the IB method for the enforcement of the nonslip boundary conditions on the immersed bed surface.The SGS effect was modeled using an eddyviscosity type model, proposed by Smagorinsky (1963), defined as where v t is the SGS turbulent viscosity and S ij = (1/2) (∂u i /∂x j + (∂u i /∂x j ) is the resolved strain rate tensor.
For the oscillatory flow considered in the present study, the Re in Eq. ( 2) is based on the orbital motion amplitude α o =U o /ω, where U o is the velocity amplitude, ω=2π/T is the radial frequency and T is the period.Also, the dynamic pressure p = P + δP, is split into the dynamic pressure P of the external oscillatory flow and the dynamic pressure correction δP due to the presence of solid bodies.

Sediment Transport
In the present study, both bed load and suspended sediment load are considered.The bed load transport rate q b is expressed by the normalized variable Φ b according to the expression: where S is the sediment specific gravity and g is the gravity acceleration.Φ b is computed by the semiempirical formulation by Engelund and Fredsøe (1976), where μ d is the dynamic friction coefficient, θ is the Shields parameter and θ c is the critical Shields parameter defined according to van Rijn (1993).The Shields parameter represents the non-dimensional bed shear stress, τ b ,   where ρ is the fluid density.
The suspended sediment concentration, c, is computed using an advection-diffusion equation where w s is the settling velocity of the sediment, calculated according to Hallermeier (1981), σ t is the turbulent Schmidt number, assumed here to be equal to 1 and f c is the source term associated with the implementation of the IB method for the enforcement of the boundary conditions on the bed surface.
Here, a Dirichlet bottom boundary condition for Eq. ( 7) is used, which corresponds to the sediment concentration on the bed level available for suspension.According to the model in Fredsøe and Deigaard (1992) used here, this sediment concentration, also known as the reference concentration, c b , is given by the expression where λ is the linear concentration and c o = 0.65 is the maximum value of the volumetric concentration of spheres with equal diameter.The suspended load flux, q s , is computed by integration over the water column ( , , ) ( , , , ) ( , , , )

Morphodynamics
The bed evolution is computed by the numerical solution of the sediment conservation equation (also called Exner equation), expressed as 12 0 where h is the bed level, q x,y is the total sediment flux and n (= 0.4) is the bed sediment porosity.

NUMERICAL METHODS
The Navier-Stokes equations were discretized through a two-stage, time-splitting scheme.In the first stage, an intermediate velocity field is computed where m is the time level, Δt is the time step, H is a spatial operator, which includes the convective and the viscous terms, based on a 2 nd order in time Adams-Bashforth scheme, and P is the dynamic pressure of the external flow.In the second stage, the computation of the velocity at the next time-step is expressed as The dynamic pressure correction, δP, is computed by solving a Poisson equation, which enforces the satisfaction of the continuity equation the velocity at the next time step.
For the spatial disretization and the imposition of the boundary conditions, 2 nd order central finite differences on a structured Cartesian grid and the Immersed Boundary method (IB) were utilized.The basic idea of the IB method, which is explained in detail in Balaras (2004), is that solid boundary surfaces, like the bed surface in the present application, do not coincide with any grid lines but are immersed in the computational domain.Therefore, the computational domain includes both the fluid and the solid phases.The boundary conditions are not applied exactly on the bed surface, which does not coincide with any grid lines, but the solution is appropriately reconstructed in the vicinity of the immersed boundary.The main advantage of the IB method is that efficient Cartesian solvers can be applied for the large-eddy simulations of turbulent flows in the presence of complex boundaries.The solution is reconstructed on the so-called forcing points, which are the grid points that belong to the fluid phase and are the nearest ones to the immersed boundary.The implementation of the IB method is as follows: -The intermediate velocity is computed on every grid point by Eq. ( 12) setting the source term f = 0.
-The intermediate velocity is reconstructed at the forcing points using a trilinear interpolation to compute the appropriate value of f, which enforces the desired boundary conditions on the immersed boundary.-The dynamic pressure correction, δP, is computed by solving a Poisson equation with appropriate boundary conditions on the computational domain boundaries but not on the immersed boundary.-The final velocity is computed on every grid point by Eq. ( 13).
A similar procedure is followed for the numerical solution of Eq. ( 7) for the suspended sediment concentration and the imposition of the boundary condition in Eq. ( 8).
The morphologic evolution of the bed is solved using a 3 rd order in time Adams-Bashforth method, similar to the one used in Ferziger and Peric (1999) and Niemann et al. ( 2011) using the spanwise-mean flow conditions.After discretization, Eq. ( 10) is expressed as where Δt m is the morphological time step.The spatial gradient, ∂q x /∂x 1 , is discretized using a simple 2 nd order central difference scheme In order to avoid numerical instabilities, a filter was applied on the bed level solution, using a twostep scheme, reported by Richards and Taylor (1981) In the above function, hu is the unfiltered bed level, ĥ is the bed level after the first step of the morphological filter and h is the resulting bed level after the application of the filter.As reported in Richards and Taylor (1981), the specific filtering application results in changes of less than 2 percent, a value usually well within the accuracy of the measurements.Moreover, an avalanche function was included in the bed load computation, in order to avoid nonphysical bed slopes that exceed the angle of repose of sediment.The procedure followed here is similar to the one used in Niemann et al. (2011).

Validation
Numerical simulations of oscillatory flow and sediment transport over fixed vortex ripples were conducted and the results were compared with numerical and experimental results n the literature for validation.The hydrodynamic results were compared with the numerical results of Blondeaux and Vittori (1991) and Scandura et al. (2000).In Fig. 1, the development of spanwise vorticity at the end of an oscillatory period, starting with fluid at rest, is shown, and the results are in good agreement with the ones in Scandura et al. (2000).In Fig. 2, the time-and bed-averaged suspended sediment concentration profile resulting by the sediment transport over vortex ripples, according to the experimental conditions Mr5b63 in van der Werf et al. (2007), is shown.The present result, obtained after 100 periods of oscillatory flow simulation and time-averaging over the last 20 ones, is in very good agreement with the experimental measurements in van der Werf et al. (2007).

Ripple creation from a flat bed
Here, the morphologic evolution of a movable bed was examined under hydrodynamic/sediment forcing.Results of coupled simulations of oscillatory flow and sediment transport are presented for sediment with normalized median grain diameter D g /α o = 2•10 -3 .The conditions were chosen so that a certain number of ripples were expected to be generated along the length of the computational domain, according to the following empirical formula in Nielsen (1992) 0.34 2.2 0.345 for the correlation of ripple length to the mobility parameter for regular waves and sediments with density close to that of quartz sand (S = 2.65).Specifically, we studied the case with ψ = 20, which according to Eq. ( 17) corresponds to ripples with L r /a o = 1,2446.The streamwise and spanwise lengths of the computational domain were set equal to 2.5a o in order to capture the development of two ripples in the streamwise direction and the development of streamwise vortical flow structures.The depth of the computational domain was set equal to 4a o in order to ensure that the proximity of the upper boundary does not affect the development of the boundary layer over the ripples.The grid spacing in the streamwise and spanwise direction, Δx 1 /α o and Δx 2 /α o , is uniform and equal to 0.01, while the grid spacing in the vertical direction, Δx 3 /α o , varies from 0.002 to 0.01, as one moves away from the bed.For the solution of the Poisson equation for δP, a Dirichlet condition δP = 0 is enforced at the bottom surface of the computational domain, which is in the solid phase below the immersed bed surface, a Neumann condition d(δP)/dz = 0 (rigid-lid) at the top surface of the computational domain, and periodic conditions in the streamwise and spanwise directions.
In Fig. 3, the spatio-temporal evolution of an initially flat bed surface with a perturbation in the middle of the computational length is shown at Re = 10,000.The initial perturbation grows fast in height during the first 100 wave cycles with smaller perturbations created on each side of it.In time, more perturbations are generated and grow (ripple creation and growth), while others are getting smaller, until they finally disappear (ripple annihilation).In some cases, ripple crests are approaching each other (see at about 600T), and they finally merge (see at about 700T) into one larger ripple (ripple merging).The rippled bed converges to its equilibrium geometry, composed of 2 ripples with length of 1.25•a o , at about 700 wave cycles.After that time, the ripples migrate with a small but positive velocity of about 410 -4 U o .The resulting bed profile displays crests sharper than troughs.This is similar to observations with real ripples under oscillatory flow (Sleath, 1984;Vittori & Blondeaux, 1991).A simulation over a longer computational domain of length 5a o for the same ψ = 20 was also performed in order to diminish the influence of the streamwise periodicity on the development of ripples.This time four ripples of length of 1.25•a o , were expected to be developed in the computational domain.In Fig. 4, ripple creation and evolution is presented for the case of the longer domain.Again, the initial perturbation grows quickly in the beginning, with smaller perturbations growing on each side of it.At about 200T, the domain is covered with ripples, while at about 300T, the bed consists of 5 ripples, which is close to the predicted final geometry.After that time, the morphological process is slower, as the net sediment transport rate for ripples near equilibrium is small (Marieau et al., 2008).At roughly 400 wave cycles, the 3 rd of the 5 ripples slowly decreases in height.The rippled bed finally converges to its equilibrium geometry, consisting of 4 ripple lengths, at about 800 wave cycles, resulting in ripple lengths of 1.25•a o .These lengths are again in agreement with the wavelengths predicted by the empirical relationship Eq. ( 17).After that time, the ripples migrate with a small positive velocity of about 410 -4 U o as in the case of the shorter computational domain.

Ripple creation from a random sinusoidal bed
Simulations were also conducted starting from initially random bed forms.Fig. 5 illustrates the evolution of an initial bed, consisting of 10 sinusoidal ripples with dimensions: L r /a o = 0.5 and h r /a o = 0.07.The simulation set-up, apart from the initial bed form, is the same as in the previous section, meaning that the computational domain has length 5a o and the mobility parameter is ψ = 20.In this case, the initial bed form is more energetic than the almost flat ones presented in the previous section, so the equilibrium bed profile is expected to be reached sooner.The initial ripples are growing and merging to form larger ripples.At 350 wave cycles, the bed geometry is composed of 5 ripples, which is close to the predicted final form.As mentioned above, the morphological process is slower as the transport rate for ripples near equilibrium is small.At about 450 wave periods, the 4 th of the 5 remaining ripples slowly decreases in amplitude, until the bed form reaches an equilibrium state, consisting of 4 ripples with wave lengths of 1.25 • a o , at about 550 wave periods.This quasi-stable bed geometry is very similar to the one obtained at the end of the evolution from the quasi-flat bed (Fig. 4).This shows that the bed tends to reach the same equilibrium state, regardless of the initial bed form.
In all the above cases, the contribution of bed load to the total sediment transport rate was dominant (more than 90% of the total load) in comparison to the suspended load.

CONCLUSIONS
LES of oscillatory flow coupled with sediment transport over a movable bed were performed.Emphasis was placed on the development of a morphodynamic model to perform simulations of the evolution of a moving bed.Ripple development from quasi-flat or random beds were examined, until the bed form reaches an equilibrium shape, adapting to flow conditions based on ψ.Simulations showed that the numerical model is able to demonstrate phenomena of ripple creation, growth, annihilation, merging and migration, resulting in ripple lengths in agreement with those predicted by empirical equations.Under the same hydrodynamic forcing, the bed tends to reach the same equilibrium state, regardless of the initial bed form.The only difference is that the bed reaches the equilibrium state faster when the simulations are initiated from a random bed form than a quasi-flat one.

Figure 2 .
Figure 2. Time-and bed-averaged concentration of suspended sediment induced by oscillatory flow over vortex ripples: present numerical results (red line) and experimental measurements (black line) in van der Werf et al. (2007).

Figure 3 .
Figure 3. Ripple creation and evolution in time and space initiating from a flat bed with a perturbation in the middle of the computational domain of length 2.5ao and driven by an oscillatory flow and sediment transport with ψ = 20.

Figure 4 .
Figure 4. Ripple creation and evolution in time and space initiating from a flat bed with a perturbation in the middle of the computational domain of length 5ao and driven by an oscillatory flow and sediment transport with ψ = 20.

Figure 5 .
Figure 5. Ripple creation and evolution in time and space initiating from a sinusoidal rippled bed along the computational domain of length 5ao and driven by an oscillatory flow and sediment transport with ψ = 20.