NUMERICAL INVESTIGATION OF FINE SEDIMENT TRANSPORT IN AN OSCILLATORY CHANNEL

Recent findings on a diverse range of muddy seabed states revealed by 3D, turbulence-resolving simulations are first reviewed. These transitions have critical implications to offshore delivery of fine sediment in the ocean and wave dissipation. Assuming a small particle Stokes number, the Equilibrium approximation to the Eulerian two-phase flow equations is applied. The resulting simplified equations are solved with a high-accuracy pseudo-spectral scheme in an idealized oscillatory bottom boundary layer (OBBL). For a typical energetic muddy shelf, the Stokes Reynolds number Re∆ is no more than 1000 and all of the scales of flow turbulence and their interaction with sediments are resolved. With increasing sediment availability or settling velocity, the seabed state evolves from well-mixed sediment distribution, to the formation of lutocline and a complete laminarization of the OBBL. More recently, we further include rheological stress in the simulations in order to study the interplay between turbulence and rheology in determining the flow regimes and hydrodynamic dissipation. To include rheological stress, we extend the numerical model with a hybrid spectral and compact finite difference scheme. A sixth-order compact finite difference is implemented in vertical direction to keep the spectral-like accuracy. The model is validated with analytical solutions using simple Newtonian rheology in laminar condition. Preliminary results at Reδ=600 reveal that when rheology is incorporated, high viscosity can trigger earlier laminarization of OBBL. When OBBL is laminarized, sediments settle and higher concentration is accumulated near the bed that further enhances viscosity and hydrodynamic dissipation. Our preliminary finding that rheology encourages laminarization may explain why large attenuation of surface waves over muddy seabed is ubiquitous and the highest dissipation rate is often observed during the waning stage of a storm.


INTRODUCTION
The state of muddy seabed in the coastal environments is the critical information to further evaluate the hydrodynamic dissipation, and their impact to coastal habitat due to pollutant, gas, and nutrient exchange.There are two challenges in modeling wave-induced fine sediment transport.Firstly, when sediment concentration is large, sediments can damp the flow turbulence and the resulting sediment concentration profile is significantly different from the commonly known Rouse profile.Secondly, for fluid mud driven primarily by waves on the continental shelves, the wave boundary layer itself is not fully turbulent.For example, the near bottom velocity measured by Traykovski et al. (2000) during significant fluid mud events is of free-stream oscillatory velocity amplitude U ~no more than 0.55 m/s and a wave period around T=10 sec (wave frequency ω=2π/T=0.628).The corresponding Stokes Reynolds number, defined as ν ∆ δ U Re = where ω ν ∆ 2 = is the Stokes length scale and ν is the kinematic viscosity of water, is no more than 1000 (or wave Reynolds number Re w <5×10 5 ).According to Hino et al. (1976), the wave boundary layer at such low Reynolds number is only intermittently turbulent.Hence, the transitional nature of the turbulent boundary layer along with the attenuation of turbulence by the presence of mud makes modeling of wave-driven fluid mud transport difficult when using conventional Reynolds-averaged turbulence closure.
Utilizing the fact that fine sediments in water are of very small particle response time, the equilibrium Eulerian approximation (Ferry and Balachandar 2001) has been applied to the two-phase formulation of turbidity currents for dilute flow (Cantero et al. 2007(Cantero et al. , 2009)).Ozdemir et al. (2010) recently adopted this formulation for fine sediment transport in the wave boundary layer and retained only the leading order term in the equilibrium Eulerian approximation.The resulting equations incorporate the feedback effect of sediment on carrier flow turbulence only due to sediment-induced density stratification.Using a highly accurate pseudo-spectral scheme to resolve all the scales of the carrier fluid turbulence and the turbulence-sediment interaction at Re δ =1000, their simulation results reveal the existence of four flow regimes (or possible seabed states) of wave-induced fine sediment transport due to the variation of sediment availability and settling velocity (Ozdemir et al. 2011): 1) a well-mixed sediment concentration with no modulation in turbulence in very dilute flow (~1 g/L or smaller), 2) formation of sharp sediment concentration gradient in the water column, i.e., lutocline, at near bed sediment concentration of O(1~10) g/l, but flow below the lutocline remains turbulent, 3) nearly laminar concentration profile of both the sediment and fluid phases but followed by burst events due to shear instability during flow reversal at near bed concentration of several tens per liter, 4) at O(100) g/l or greater, a complete laminarization throughout the wave cycle.The existence of these flow modes has critical implications to our capability in assessing the state of the muddy seabed and to further understand various applications related to fluid mud transport.For instance, regime I represents a significantly confined, but highly mobile transport close to the seabed.Understanding the transition between regime I and regime II and to further quantify the evolution of lutocline in regime II allow us to estimate the amount of offshore fine sediment transport via wave-supported gravity flow (e.g., Traykovski et al. 2000).Moreover, regimes III and IV represent significantly calm, gooey and less mobile fluid mud.Field observation suggests when regimes III and IV occur, significant damping of surface wave energy occurs.
To further investigate the interplay between the collapses of turbulence, the rapid settling of sediment associated with laminarization, and more importantly, the initiation of rheological stress due to the rapid increase of sediment concentration near the bed, we extend the numerical model of Ozdemir et al. (2010Ozdemir et al. ( , 2011) ) with the capability of modeling rheological stress.The governing equations and the modification of the numerical scheme necessary for including rheological stress are briefly reported in Section 2. Model results revealing the effect of rheological stress based on a Newtonian viscosity (viscosity parameterized as fast increasing function of sediment concentration) are reported in Section 3. Conclusions and future works are discussed in Section 4.

MODEL FORMULATION
We idealize the case of wave-induced turbulence at the bottom boundary layer in a typical coastal setting.For typical fine sediment transport in the continental shelf, the settling velocity is around 0.1~1 mm/s (Hill et al. 2000).Hence, we consider silt with grain size of d=24 µm and a specific gravity of s=2.65.This results in a still fluid Stokes settling velocity of about W s =0.5 mm/s.The particle response time, defined as ν τ 18 2 sd p = can be computed to be 8.5×10 -5 s.For the wave condition considered in the study with Re δ =400~1000, the timescale of the Stokes layer, , is more than an order of magnitude larger than the particle response time.Thus, it can be established that the particles are sufficiently small that their Stokes number, defined as St=τ p /τ f , is much less than 1.As a result, here we ignore the inertial effect of particles (Ferry and Balachandar, 2001) and simply define the nondimensional particle velocity as the sum of fluid velocity, f i u and the particle settling velocity: (1) where the superscripts f and s stand for the fluid and particle phases, respectively and V s is the nondimensionalized particle settling velocity.In this study, the free-stream velocity amplitude U ~and Stokes length scale are chosen as the characteristic velocity scale and length scale to nondimensionalize the flow quantities.Moreover, we define the volume-averaged concentration over the entire volume ∀ of the computational domain, to normalize sediment volumetric concentration ϕ * .Hence, for the remaining of this paper, the normalized sediment volumetric concentration is used.By substituting the algebraic relationship shown in equation ( 1) into the Eulerian-Eulerian two-phase equations and applying the Boussinesq approximation, which is valid for dilute sediment transport in water, the resulting governing equations can be greatly simplified (e.g., Cantero et al. 2007).The continuity and momentum equations of the fluid phase can be written as: where i=1, 2, 3 representing streamwise, spanwise and vertical direction, respectively and p is the nondimensionalized dynamic pressure, A represents a mean (normalized) pressure gradient that drives the flow and f(ϕ) is a Newtonian rheology function, which will be discussed later.In the present problem of an oscillatory channel flow, the time variation of the dimensional free-stream velocity is given as sinusoidal flow with amplitude U ~ and angular frequency ω.In the free stream, away from the bottom boundary layer, the momentum equation can be simplified to inviscid flow to obtain the flow forcing (5) The third term on the right-hand-side of equation ( 4) represents the feedback effect of the sediment phase on the carrier fluid via sediment-induced density stratification.Here the bulk Richardson number Ri is defined as: In this study, we consider the simple Newtonian rheology proposed by Einstein (1906) (7) This linear relationship is valid for mono-sized spherical particles in dilute condition.For this study, the Φ is chosen to be 1% and dilute assumption is satisfied.
Finally, the normalized sediment volumetric concentration is calculated by the conservation of sediment mass: (8) The term on the right-hand-side represents an effective diffusive flux.Here Sc is the Schmidt number of the sediment phase.The governing equation of particle motion, when applied in the Lagrangian framework, contains no diffusive effect, i.e.Sc=∞.However, in the present Eulerian simulations, a diffusive term is required not only for numerical stability but also to account for subscale particle motion through sub-grid interactions.
The velocity boundary conditions at both the top and the bottom boundaries are no-slip wall boundary conditions.Periodic boundary conditions are imposed at the lateral boundaries of the computational domain.For particle concentration, the following relation is imposed at both top and bottom boundaries of the computational domain: (9) This boundary condition imposes no net transport of sediment across the top and bottom boundaries.Periodic boundary conditions for sediment concentration are applied along the streamwise and spanwise directions.These boundary conditions in turn imply that the total volume of particles remains constant in the domain throughout the computation (Φ is constant throughout the computation).This integral property has been confirmed as a verification test in the present computations.
The size of the computational domain is selected to be 60Δ×30Δ×60Δ along the streamwise, spanwise and vertical directions, respectively.Thus, the distance from the bottom boundary to the center plane of the channel is chosen to be 30 times the Stokes layer thickness.The selected domain size is exactly the same as that of Spalart & Baldwin (1989) and their domain is the largest and the most conservative among the different numerical studies.By examining the two point correlation function, Ozdemir et al. (2010) showed that this domain size is sufficiently large to capture the main turbulent characteristics at Re Δ =1000 in the case of clear fluid.For Re ∆ =600, the computational domain is discretized into 128×128×193 elements.Hence, the total amount of grid points is around 3.2 millions.
To ensure high numerical accuracy, the numerical scheme adopted by Ozdemir et al. (2010Ozdemir et al. ( , 2011) is based on a pseudo-spectral scheme.However, it is difficult to further incorporate particle stresses (rheology) in a pseudo-spectral scheme.Recently, this code is revised with a sixth-order compact-finite difference scheme in the vertical direction in order to investigate the interplay between turbulence and rheology in determining the resulting flow regimes and flow dissipation.A general compact finite difference scheme for one-dimensional problem has the form (Lele 1992) (10) where u j are the function values given at the grid points x j on set I n and I m , and u j (p) are the values of the p-th derivative of the function at the corresponding grid points.To resolve the turbulent boundary layer, Chebyshev collocation points are used, where the grid points cluster at both the top and bottom boundaries. (11) For our study, a tri-diagonal centered sixth order compact finite difference scheme is applied on the interior grid points, where the points sets I n = {i-1, i+1}, and I m ={i-2, i, i+2,}.For the grid points near the top and bottom boundaries, a one-sided scheme with the same order of accuracy is implemented.The detailed derivation of the compact finite difference scheme on non-uniform grid can be found in Shukla & Zong (2005).The first and second order derivatives, therefore, can be written as (12) where A, A 2 are tri-diagonal matrices, and B, B 2 are banded matrices.
The flow field is advanced using Crank-Nicolson scheme for the diffusion terms.The non-linear advection terms are calculated by the Arakawa method with the 3/2 de-aliasing law, and integrated by a third-order low-storage Runge-Kutta scheme.More details on the implementation of this numerical scheme can be found in Cortese and Balachandar (1995).
Although the rheology closure adopted in the present study is quite idealized comparing to realistic mud, it allows as to begin to evaluate the effect of Newtonian rheology on the flow regimes via an enhanced viscosity due to sediment concentration.More important, the simple linear rheology closure used in equation ( 7) allows us to derive analytical solution for laminar condition in order to validate the numerical scheme.Excellent agreement between the numerical results and analytical solutions is obtained.

RESULTS
Simulation is carried out at Re ∆ =600, which represents typical wave condition of typical low energy muddy continental shelves with U ~=0.4 m/s and wave period T=7.0 sec.Sediment availability is specified as Ri=1.0×10 - with a nondimensional settling velocity V s =4.5×10 -4 .According to prior studies, such as Spalart & Baldwin (1989) and our own direct numerical simulations (for clear fluid), flow at Re ∆ =600 is not fully turbulent.However, flow instabilities occur during flow reversal that further leads to chaotic motion before the following peak flow.Hence, ensemble-averaged sediment concentration profiles are similar to the regime II condition discussed in Ozdemir et al. (2010Ozdemir et al. ( , 2011) ) where the lutocline feature is observed that further separates the lower turbulent region from the upper quasi-laminar region (see Figure 1 and Figure 4(b1) and ( b2)).
It should be first mentioned here that because the flow is not fully turbulent for Re ∆ =600, adding sediments further damps flow turbulence via sediment-induced density stratification (Ozdemir et al. 2011) and we found that when computing the simulations for sufficiently amount of wave periods, flow eventually laminarizes.This feature is quite different from that observed for Re ∆ =1000 in more turbulent condition.In the present simulation of Re ∆ =600, flow eventually laminarize after 25 wave period.Hence, the key issue to be investigated is whether adding rheological effect may delay or encourage laminarization.
To illustrate the turbulent state of the calculated flow field, the turbulent coherent structure is identified by Q-method (Hunt et al. 1988), where the positive second invariant Q of velocity gradient is used: (13) A small value (very close to 0) is typically used to identify the coherent structure (Q=0.2 is used here).To further illustrate the mixing and suspension of sediment near the bed, the iso-surface of normalized sediment concentration of ϕ=2.5 is shown.These flow visualizations during peak flow and flow reversal are presented in Figure 2 and Figure 3, respectively.Comparing to the model results without rheology, we find that the present Newtonian rheological stress with an enhance viscosity tends to attenuate more turbulence and eventually trigger earlier laminarization.Preliminary results reveal that when rheology is incorporated (Figure 2(b1), (b2)), wave boundary layer during flow peak is less turbulent comparing to that without rheology (Figure 2(a1), (a2)).Similar features are also observed during flow reversal (see Figure 3).Reduced turbulence at Re ∆ =600 due to rheology, eventually leads to earlier laminarization.For the present simulation with Einstein rheology, laminarization occurs quite early.It only takes four wave periods for laminarization (comparing to twenty wave periods when rheology is excluded).Laminarization further causes more significant settling, higher sediment concentration accumulated near the bed, and finally results in enhanced hydrodynamic dissipation.Our preliminary finding that rheology encourages laminarization may explain why large attenuation of surface waves over muddy seabed is ubiquitous and the highest dissipation rate is often observed during the waning stage of a storm (A.Sheremet (U.Florida) and P. Traykovski (WHOI), personal communications).More quantitative comparison can be made by examining the ensemble-averaged flow velocity, sediment concentration and turbulence kinetic energy (TKE) profiles during flow peak and flow reversal (see Figure 4).In the present simulation, ensemble average is carried out via spatial averaging over the x-y plane.Sediment concentration profiles are more or less similar throughout the wave cycle, which is representative of fine sediment with small settling velocity.A sharp negative concentration gradient, namely the lutocline, is observed in the upper region of the wave boundary layer.TKE is basically shot down in the lutocline region and above but significant turbulence is observed below the lutocline (see Figure 4(c2)).During flow peak, it is quite clear that TKE near the bed with rheology incorporated is 30%~50% lower than that without rheology (see Figure 4(c1)).Due to lower turbulence, sediment concentration gradient very near the bed is also larger when rheology is included.More interesting features are observed during flow reversal.When rheology is not considered, the region below the lutocline is turbulent with significant TKE.On the other hand, when rheology is included, TKE is almost zero throughout the region below the lutocline except very near the bed where strong TKE is observed.In other words, when rheology is included, the region of significant turbulent production is greatly suppressed.Since the intensity of turbulence scales with the length scale of the turbulent region, this explains why turbulence is much lower when rheology is incorporated.More detailed analyses on the TKE budget and the transitions between the more turbulent region and quasilaminar region are necessary in order to pinpoint the mechanisms causing the reduced turbulence and subsequent earlier laminarization due to rheological effects.

CONCLUSION
Turbulence resolving simulations reveal the existence of four flow regimes from well-mixed condition, to the formation of lutocline and a complete laminarization due to increasing sediment-induced density stratification.By implementing a Newtonian rheology of Einstein (1906) where the viscosity increases with sediment volumetric concentration, we demonstrate that the enhanced viscosity tends to damp more turbulence and in some cases to trigger earlier laminarization.Hence, there is a strong and complex interplay between boundary layer turbulence, turbulence modulation by sediment and rheological effect during wave-mud interaction.Future studies will focus on identifying the key mechanisms to explain the damping of turbulence due to Newtonian rheology and to quantify the energy dissipation.Non-Newtonian rheology (shear thinning or Bingham plastic) should be incorporated to investigate more realistic mud transport processes.Higher order inertia effects should be further included in order to evaluate the effect of preferential concentration and extended the applicability of the present numerical model to coarser grains.

Figure 1 .
Figure 1.A summary on the existence of four flow regimes revealed by Ozdemir et al. (2010).These regimes shown above are due to different sediment availability under the same wave intensity (Stokes Reynolds number Reδ=1000) and settling velocity (0.5 mm/s).