Longitudinal dispersion of heavy particles in an oscillating tunnel and application to wave boundary layers

Abstract

The present research aims at getting an understanding of the process of dispersion of surface sediment in an oscillatory boundary layer, which may represent an idealised case of, for example, a stockpile area where excavated sediment is stockpiled temporarily (or permanently). The process is studied numerically, using a random-walk particle model with the input data for the mean and turbulence characteristics of the wave boundary layer picked up from a transitional two-equation k–$\omega$ Reynolds averaged Navier–Stokes model and plugged in the random-walk model. First, the flow model is validated against experimental data in the literature. Then, the random-walk dispersion model is run for different oscillatory flow cases and for a number of steady flow cases for comparison. The primary sediment grains of concern are fine sediments (with low fall velocity), which would stay in suspension for most of the time. Nevertheless, the dispersion of neutrally buoyant and heavier particles that spend most of their time in close vicinity to the bed are also discussed. The numerical model results are compared with the results of a series of experiments carried out in an oscillating U-tunnel facility. The results are found to be in general agreement both qualitatively and quantitatively. In the last part of the study, an example application of the present model for fine sand dispersing in a wave boundary layer under storm conditions is given.

1 Introduction

Construction operations in the marine environment may require vast amount of dredging by means of different methods. During the handling and disposal process of excavated (dredged) sediment, various amounts of this rather fine material will be spread into the marine environment. Eventually, the dredged material is either used as fill material for other parts of the construction or is disposed to a stockpile area in the sea. Since the marine environment, especially the near-bed region, is highly dynamic in terms of current and/or wave action, the disposed sediments will go through a series of transport and mixing processes (i.e. advection, diffusion, dispersion). It is of great importance to understand the mechanics of these processes for planning and managing of marine dredging operations.

Dispersion in bottom boundary layers is the most dominant of all these mixing processes qualitatively and quantitatively. In this process, the joint action of shear (mean velocity gradient) and wall-normal turbulence fluctuations in the boundary layer near the bed disperses the sediment particles horizontally. The dispersion coefficient is the measure for the rate of dispersion.

Longitudinal dispersion in steady flows has, in the past, drawn quite a lot of attention, with the pioneering work of Taylor (1953, 1954), followed by Aris (1956), Elder (1959), Fischer (1966, 1967), Sayre (1968), Chatwin (1973), Sumer (1973, 1974), Pedersen (1977), Fischer et al. (1979), Chatwin and Sullivan (1982), Allen (1982), Smith (1983) and Demuren and Rodi (1986) among others. The findings of these studies were often applied to channels, transmission pipes, atmospheric boundary layers, rivers, tidal inlets and estuaries.

Wave boundary layers, on the other hand, are formed over the seabed under waves. They are characterised by a very small vertical extent (in the order of magnitude of 20–30 cm at most), with very strong shear and turbulence, two important “ingredients” for dispersion. Unlike the steady flow situations, there are not many studies dealing with dispersion processes in unsteady flows, particularly in oscillatory and wave boundary-layer flows. Aris (1960) was the first to deal with the dispersion of a solute in periodically altering flows, albeit for the laminar case. Chatwin (1975) investigated the longitudinal dispersion of a passive neutrally buoyant dispersant in oscillatory pipe flow. Smith (1982), Yasuda (1982) and Yasuda (1984) studied the dispersion process in oscillatory two-dimensional boundary layers, which resembles tidal flow boundary layers rather than wave boundary layers. Yasuda (1989) extended his previous research to dispersion of suspended (heavy) particles. Mei and Chian (1994) conducted a theoretical study in which they paid special attention to the dispersion of small suspended particles in wave boundary layers, including the convection due to the steady streaming. These works all used a time-invariant eddy viscosity when describing the turbulent flow. Ng (2004) studied the dispersion process in oscillatory boundary layers using a time-dependent turbulent diffusivity and found that the time dependency of the eddy viscosity could not be neglected since it generated a significant difference in terms of dispersion coefficients. Recently, Mazumder and Paul (2012) carried out a numerical study in which they modelled the dispersion of suspended sediments in oscillatory boundary layers including the settlement, temporary storage and re-entrainment processes.

The present research aims at getting an understanding of the process of dispersion of the surface sediment in an oscillatory boundary layer, which may represent an idealised case of, for example, a stockpile area where excavated sediment (stockpiled temporarily or permanently) is subject to waves. Section 2 presents a single-particle analysis for the longitudinal dispersion, first for a steady turbulent flow in an open channel, followed by an extension of the latter analysis to an oscillatory turbulent flow in an oscillating tunnel. Section 3 presents the numerical model adopted in the study, comprising a k–$\omega$ model to calculate the “background” flow, and a random-walk model to calculate the longitudinal dispersion. Section 4 describes the single-particle experiments conducted in an oscillating tunnel. Section 5 presents the results including the validation of the numerical model, and comparison with the experimental data obtained in the experimental campaign. The paper ends with a final section, Sect. 6, in which the numerical model developed in the present study is applied to a real-life problem with dispersion of disposed sediment in a wave boundary layer under storm conditions. The present results and the existing knowledge (yet, highly limited) form a complementary source of information on the dispersion of sediment in wave/oscillatory boundary layers.

2 Longitudinal dispersion in an oscillating tunnel

In this section, the description of the longitudinal dispersion process will be discussed from the point of view of one-particle analysis. First, the dispersion process will be discussed for an open channel flow (Sect. 2.1) and, subsequently, the analysis will be extended to the idealised case of dispersion in an oscillating tunnel (Sect. 2.2).

2.1 Longitudinal dispersion in an open channel flow

We first consider the turbulent flow in an open channel with a smooth bottom. The flow depth is h and the x and y axes are chosen so that x is in the flow direction and y in the upward direction (Fig. 1). We will describe the process of longitudinal dispersion of sediment particles, using the so-called one-particle analysis. We release the particle from a point in the channel (Fig. 1). The particle will travel under the combined action of the mean shearing motion, turbulence and gravity. We consider that the flow conditions and the particle properties are such that the particle, although heavy, is not deposited on the bottom when it comes close to the bottom, but is reflected back into the main body of the flow.

Fig. 1

Definition sketch. Longitudinal dispersion in an open channel flow

Because of the presence of the free surface and the bottom, and more importantly the conditions ensuring that the particle remains in the flow, the horizontal velocity of the particle is necessarily a stationary random function of time as soon as the influence of the special choice of the point on the cross-section where the particle was released has been lost.

Sumer (1974) developed a statistical formulation to describe the longitudinal dispersion of heavy particles in such a flow environment. In what follows is a summary of this formulation, which will be extended to the case of oscillatory flow in an oscillating tunnel in the next section.

Let U(t) denote the fluctuation about the ensemble mean $\overline{U}$ of the component of the particle velocity (Lagrangian velocity) in the flow direction. The mean position of the particle in the horizontal direction is a distance

$$\begin{aligned} \overline{X} (t)=\overline{U} \;(t-t_0 ) \end{aligned}$$

(1)

downstream from the point of release. Here, t is the time and $t_{0}$ is the time at which the particle is released into the flow. The variance of the displacement about the mean, $\overline{X}$, is

$$\begin{aligned} \overline{X^2} (t )\cong & {} 2(t-t_0 )\;\overline{U^2} \int \limits _0^\infty S(t')\;\mathrm{d}t'\nonumber \\&-2\overline{U^2} \int \limits _0^\infty {t'S(t')\;\mathrm{d}t'} \quad \quad \text{ as }\;\;t-t_0 \rightarrow \infty \end{aligned}$$

(2)

in which $S(t')$ is the autocorrelation coefficient of the velocity U:

$$\begin{aligned} S(t')=\overline{U(t)\;U(t+t')}/\overline{U^2} \end{aligned}$$

(3)

The probability density function of the particle position in the horizontal direction tends to a Gaussian distribution with a longitudinal diffusivity $D_{1}$ (longitudinal dispersion coefficient) as $t-t_0 \rightarrow \infty$ where

$$\begin{aligned} D_1 =\frac{1}{2}\frac{\mathrm{d}\overline{X^2} }{\mathrm{d}t}=\overline{U^2} \;\int \limits _0^\infty {S(t')\;\mathrm{d}t'} \end{aligned}$$

(4)

Sumer (1974), on dimensional grounds, expressed the mean particle velocity, $\overline{U}$, and the longitudinal dispersion coefficient, $D_1$, by

$$\begin{aligned} \frac{\overline{U} }{U_\mathrm{f} }=f(\beta ),\quad \text{ and }\quad \frac{D_1 }{hU_\mathrm{f} }=g(\beta ) \end{aligned}$$

(5)

in which $U_\mathrm{f}$ is the friction velocity and $\beta$, also called the Rouse parameter, is defined by

$$\begin{aligned} \beta =\frac{w_\mathrm{s} }{\kappa U_\mathrm{f} } \end{aligned}$$

(6)

in which $\kappa$is the von Karman constant ($\kappa =0.4)$. It may be noted that Sumer (1974) calculated analytically the latter quantities, using Aris’ moment method for the range of $0\le \beta <0.7$, as will be discussed later.

2.2 Longitudinal dispersion in an oscillating tunnel

We now consider an oscillatory flow in an oscillating tunnel (Fig. 2). The oscillatory flow Reynolds number, $Re_\mathrm{w} =\frac{aU_\mathrm{m} }{\nu }$, is such that the oscillatory boundary layer is in the turbulent regime, i.e. $Re_\mathrm{w} >1.5\times 10^5$ (Sleath 1984; Jensen et al. 1989). Here, $U_\mathrm{m}$ is the amplitude of the free stream velocity $U_0 =U_\mathrm{m} \sin ( {\omega _\mathrm{w} \,t})$, $\omega _\mathrm{w}$ is angular frequency and a is the amplitude of the free stream motion, equal to $U_\mathrm{m}$/$\omega _\mathrm{w}$. Note that this is precisely the same flow environment as that used in the single-particle experiments described in Sect. 4.

Fig. 2

Definition sketch. Longitudinal dispersion in an oscillatory flow in an oscillating tunnel

Similar to the open channel case, we will describe the process of longitudinal dispersion of sediment particles in the oscillating tunnel, using the one-particle analysis. We release the particle from a point in the tunnel (Fig. 2). The particle will travel under the combined action of the phase-resolved mean shearing motion, phase-resolved turbulence and gravity. Here, too, we consider that the flow conditions and the particle properties are such that the particle is maintained in the main body of the flow.

Because of the presence of the two walls of the tunnel, and more importantly the conditions ensuring that the particle remains in the flow, the horizontal velocity of the particle will be a stationary random function of time as soon as the influence of the special choice of the point on the cross-section where the particle was released has been lost, similar to the open-channel case. Therefore, the statistical formulation developed for the open-channel flow case should be equally valid for the oscillating tunnel case as well (Eqs. (1)–(4)), with $\overline{U}$ and $\overline{X} (t)$ in Eq. (1) identically equal to zero due to the symmetric periodic motion in the tunnel.

The non-dimensional expression for the longitudinal dispersion coefficient, $D_1$ in Eq. (5), may be adopted for the oscillating tunnel case, replacing (1) the flow depth h (which is actually the thickness of the boundary layer in the open channel) with the thickness of the oscillatory boundary layer in the oscillating tunnel, $\delta$, and (2) the friction velocity $U_\mathrm{f}$ with the maximum value of the friction velocity associated with the oscillatory boundary layer in the oscillating tunnel, $U_\mathrm{fm}$. Here, $\delta$ is defined as the location of the point of maximum velocity from the bed measured at the phase $\omega _\mathrm{w}t=90^\mathrm{o}$ (i.e. when $U_{0}=U_\mathrm{m})$ (Jensen et al. 1989).

Regarding the function $g(\beta )$ on the right hand side of Eq. (5), there will be one additional parameter in the present case, namely the ratio of the flow amplitude and the tunnel half width, a / h:

$$\begin{aligned} \frac{D_1 }{\delta U_\mathrm{fm} }=g\left( \beta , \frac{a}{h}\right) \end{aligned}$$

(7)

Contrary to wave boundary layers where the boundary-layer flow is described by only one single parameter, namely the wave Reynolds number, the oscillatory boundary layer flow in an oscillating tunnel is governed by two parameters, namely the wave Reynolds number and a non-dimensional parameter involving the height of the oscillating tunnel (see Lodahl et al. 1998). In the above formulation (Eq. (7)), the wave Reynolds number is involved indirectly through the boundary layer thickness $\delta$, and the height of the tunnel is involved through the non-dimensional parameter $\frac{a}{h}$.

3 Numerical model

The numerical model has two components: (1) the flow model that will yield the phase-dependent velocity and turbulent kinetic energy profiles, and (2) the random-walk dispersion model that will allow the simulation of the dispersion process, based on the tracking of a single particle for multiple times with the input from the flow model. These components are explained below.

3.1 The flow model

Two-equation turbulence closure models for Reynolds averaged Navier–Stokes (RANS) equation have been widely used in fluid mechanic and coastal engineering for solving various problems related to turbulent boundary layers. These models either use the k–$\omega$ or k–$\varepsilon$ turbulence closure. In these models, it is possible to solve the Reynolds averaged velocities (i.e. mean velocities) along with wall shear stresses and turbulence kinetic energy on a time-averaged or phase-resolved basis.

In the present study, a modified version of the k–$\omega$ model developed by Fuhrman et al. (2013), MatRANS, has been used for the mean flow and turbulence input necessary for the longitudinal dispersion simulations. The modification of the model is a transitional flow option, which can capture the transition from laminar to turbulent flow conditions. The model solves simplified versions of the horizontal component of the incompressible RANS equations, combined with the two-equation k–$\omega$ turbulence closure model of Wilcox (2006). The considered RANS equation reads:

$$\begin{aligned} \frac{\partial \bar{u}}{\partial t}=-\frac{1}{\rho }\frac{\partial \bar{p}}{\partial x}+\nu \frac{\partial ^2\bar{u}}{\partial y^2}+\frac{\partial }{\partial y}\left( {\nu _\mathrm{T} \frac{\partial \bar{u}}{\partial y}}\right) \end{aligned}$$

(8)

where u, v and w are the streamwise, wall-normal and lateral velocity components, respectively. p is pressure, $\rho$ is specific mass of fluid, $\nu$ and $\nu _\mathrm{T}$ are laminar and turbulent kinematic viscosities, respectively. An overbar denotes the Reynolds averaging. The turbulence model consists of two respective transport equations for (1) the turbulent kinetic energy (per unit mass) $k=\frac{1}{2}\left( {\overline{{u}'^2} +\overline{{v}'^2} +\overline{{w}'^2} }\right)$, where the prime denotes fluctuating components, and (2) the specific dissipation rate $\omega$ (Fuhrman et al. 2013):

$$\begin{aligned} \frac{\partial k}{\partial t}= & {} \nu _\mathrm{T} \left( {\frac{\partial \bar{u}}{\partial y}\frac{\partial \bar{u}}{\partial y}}\right) -\beta ^*k\omega +\frac{\partial }{\partial y}\left[ {\left( {\nu +\sigma ^*\alpha ^*\frac{k}{\omega }}\right) \frac{\partial k}{\partial y}} \right] \nonumber \\\end{aligned}$$

(9)

$$\begin{aligned} \frac{\partial \omega }{\partial t}= & {} \alpha _\omega \frac{\omega }{k}\nu _\mathrm{T} \left( {\frac{\partial \bar{u}}{\partial y}\frac{\partial \bar{u}}{\partial y}}\right) -\beta _\omega \omega ^2+\frac{\sigma _\mathrm{d} }{\omega }\frac{\partial k}{\partial y}\frac{\partial \omega }{\partial y}\nonumber \\&+\frac{\partial }{\partial y}\left[ {\left( {\nu +\sigma _\omega \alpha ^*\frac{k}{\omega }}\right) \frac{\partial \omega }{\partial y}} \right] \end{aligned}$$

(10)

In Eq. (10), $\sigma _\mathrm{d}$ is defined as:

$$\begin{aligned} \sigma _\mathrm{d} =\left\{ {\begin{array}{l@{\quad }l} 0 , &{}\dfrac{\partial k}{\partial y}\dfrac{\partial \omega }{\partial y}<0 \\ \sigma _\mathrm{do} , &{}\dfrac{\partial k}{\partial y}\dfrac{\partial \omega }{\partial y}\ge 0 \\ \end{array}} \right. \end{aligned}$$

(11)

The turbulent (eddy) viscosity is defined slightly different from the previous versions of the model:

$$\begin{aligned} \nu _\mathrm{T} =\frac{\alpha ^*k}{\omega }\, \end{aligned}$$

(12)

Here, $\alpha ^{*}$ is a model coefficient adjusted by a model Reynolds number Re $_\mathrm{T}$ as follows:

$$\begin{aligned} \alpha ^*=\frac{\alpha _0^*+\frac{Re_\mathrm{T} }{R_\mathrm{k} }}{1+\frac{Re_\mathrm{T} }{R_\mathrm{k} }},\quad Re_\mathrm{T} =\frac{k}{\omega \cdot \nu } \end{aligned}$$

(13)

When $\alpha ^{*}$ is unity, the model becomes identical with that of Fuhrman et al. (2013). Similarly, $\alpha$ and $\beta ^{*}$ coefficients are also adjusted by Re $_\mathrm{T}$.

$$\begin{aligned} \alpha= & {} \frac{13}{25}\frac{\alpha _0 +\frac{Re_\mathrm{T} }{R_\omega }}{\alpha ^*\left( {1+\frac{Re_\mathrm{T} }{R_\omega }}\right) }\end{aligned}$$

(14)

$$\begin{aligned} \beta ^*= & {} \beta _0^*\frac{\frac{100}{27}\beta _\omega +\left( {\frac{Re_\mathrm{T} }{R_\beta }}\right) ^4}{1+\left( {\frac{Re_\mathrm{T} }{R_\beta }}\right) ^4} \end{aligned}$$

(15)

The values assigned to rest of the model coefficients are given as follows:

$$\begin{aligned} \begin{array}{l} \beta _\omega =0.0708 \quad \alpha _0 =\frac{1}{9} \quad \alpha _0^*=\frac{1}{3}\beta _\omega \\ \beta _0^*=\frac{9}{100} \quad R_\omega =2.61 \quad R_\mathrm{k} =3 \quad R_\beta =8 \\ \sigma _\omega =\frac{1}{2} \quad \sigma _\mathrm{do} =\frac{1}{8} \quad \sigma ^*=\frac{3}{5} \\ \end{array} \end{aligned}$$

(16)

Since the flow model does not give the turbulent fluctuations separately but only the turbulent kinetic energy (k), wall-normal turbulence fluctuating velocity, $\overline{{v}'^2}$, is estimated by the following approximation suggested by Nezu and Nakagawa (1993):

$$\begin{aligned} \overline{{v}'^2} \approx 0.34\,k \end{aligned}$$

(17)

Finally, a length scale of turbulence, $\ell$, is defined as follows:

$$\begin{aligned} \ell =\frac{\sqrt{k} }{\omega } \end{aligned}$$

(18)

The flow model is solved for the lower half of the tunnel due to symmetry around the centreline. In the model, a no-slip condition is imposed at the bottom boundary (i.e. wall), where the velocity is zero. For the turbulent kinetic energy, also a zero value is imposed at the wall (i.e. $k=0$). For specific dissipation rate, $\omega$, the following bottom boundary condition adopted from Wilcox (2006) was used:

$$\begin{aligned} \omega =\frac{U_\mathrm{f} ^2}{\nu }S_\mathrm{R},\quad y=0 \end{aligned}$$

(19)

in which $S_\mathrm{R}$ is a factor based on the roughness Reynolds number $k_\mathrm{s}^{+}=k_\mathrm{s}U_\mathrm{f}$/$\nu$ where $k_\mathrm{s}$ is the Nikuradse’s equivalent sand roughness. $S_\mathrm{R}$ is given as:

$$\begin{aligned} S_\mathrm{R} =\left\{ {\begin{array}{ll} \left( {\frac{200}{k_\mathrm{s} ^+}}\right) ^2 , &{} k_\mathrm{s} ^+\le 5 \\ \frac{K_\mathrm{r} }{k_\mathrm{s} ^+}+\left[ {\left( {\frac{200}{k_\mathrm{s} ^+}}\right) ^2-\frac{K_\mathrm{r} }{k_\mathrm{s} ^+}} \right] e^{5-k_\mathrm{s} ^+} , &{} k_\mathrm{s} ^+>0 \\ \end{array}} \right. \end{aligned}$$

(20)

Here, $K_\mathrm{r}$ is a model coefficient with the value of 50.

On the upper boundary (namely, at the centreline of the tunnel), a frictionless lid is considered, whereby the shear stress and vertical derivatives of all variables are set to zero.

The flow is driven by prescribing the pressure gradient term appearing in Eq. (8). For this study, the following values are used for the pressure gradient term:

$$\begin{aligned} \frac{1}{\rho }\frac{\partial p}{\partial x}=-\frac{\partial U_0 }{\partial t}+P_x \end{aligned}$$

(21)

$U_{0}$ is the desired free stream velocity ($U_0 =U_\mathrm{m} \sin ( {\omega _\mathrm{w} \,t}))$ and $P_{x}$ is the constant pressure gradient associated with the hydraulic slope which becomes zero in the case of pure oscillatory flow and $P_x =-U_\mathrm{f} ^2/h$ in the case of steady flow.

When the flow model was run for steady flow conditions, the model time was run for long enough such that all the flow parameters would remain unchanged with time. When oscillatory flow conditions are simulated by flow model, it was seen that after 4$\mathrm{th}$ or 5$\mathrm{th}$ wave period the model output started to repeat itself periodically. As such, the flow model was run for a warm up time of 10 wave periods and the values in the final wave period of the model were extracted for repeated use in the dispersion model.

In this version of the model, secondary (convective) terms were not included in the model equations. Therefore, when simulating the oscillatory flow, the steady streaming associated with the wave boundary layers is not considered. For other details of the MatRANS model, Fuhrman et al. (2013) can be consulted.

3.2 The dispersion model

The dispersion process is studied, using a random-walk model. The random-walk model is a particle-based model where a number of particles (N) are released from a point source (or a line or plane source) one after another, each being tracked from time $t=0$ to a time of $t=t$. A particle placed at a vertical position is convected in the vertical direction by (1) the vertical fluctuating velocity, and (2) the fall velocity of the particle itself, while it is convected in the horizontal direction by the horizontal mean velocity. The model assumes that the vertical fluctuating velocity exists for a small time interval. The magnitude and direction of the vertical fluctuating velocity are the random elements of the model. The path of an initially marked particle position is calculated as the sum of a series of such small time intervals as the particle migrates through the statistical field variables. Many such paths are used to compile the statistical properties by means of an ensemble average. In the simulations, the particles are maintained in the main body of the flow all the time. To this end, the particles coming very close to the bottom are re-entrained into the flow, as will be detailed later.

We note that no interaction between the particles (grain-to-grain collision) is taken into account in the simulations, implying that simulations are not valid for extremely high concentrations. Also, the bouncing of particles at the bed is considered to be due to near-bed turbulence alone. This is a valid approximation for sediment covering silt, fine sand and even partly medium-size sediment.

To the authors’ knowledge, Sullivan (1972) and Sumer (1973) were the first to use a random-walk model for the purpose of simulating the longitudinal dispersion processes in a turbulent open-channel flow, Sullivan (1972) for neutrally buoyant particles, and Sumer for heavy particles. Bayazit (1972) also used a similar random-walk model to calculate the settling distances of heavy particles in a turbulent open-channel flow.

In the present study, the simulations were carried out for three kinds of flow environments: (1) for a steady flow in a tunnel; (2) for an oscillatory flow in an oscillating tunnel, and (3) for an oscillatory flow in a wave boundary layer. The steady flow case is included as a reference case, which also enables the model validation as well. The model is also validated in the case of the oscillating tunnel against the laboratory experiments carried out in an oscillating tunnel, described in the next section. The third case, the wave boundary layer, is included to demonstrate the implementation of the model in a real-life situation.

The flow domain is defined as a 2D tunnel (Fig. 2). The hydrodynamic model is run for the half-space 0 first, and then the flow domain for dispersion model is defined by vertically mirroring the hydrodynamic model output about the center line.

In all the three cases, the bottom wall of the flow environment is considered to be hydraulically smooth. When particles get very close to the wall, they may get embedded in the viscous sub-layer since the turbulence fluctuations are very small there. To prevent this, the particles are not allowed to travel closer to the wall than a certain distance, $\delta _\mathrm{b}$, taken to be equal to the thickness of viscous sub-layer, similar to Sullivan (1972) and Sumer (1973):

$$\begin{aligned} \delta _\mathrm{b} ^+=5y^+ \quad \text{ or }\quad \delta _\mathrm{b} =5\frac{\nu }{U_\mathrm{f} } \end{aligned}$$

(22)

where $U_\mathrm{f}$ is the friction velocity in the steady flow case, or the amplitude of the friction velocity, $U_\mathrm{fm}$, in the oscillatory flow cases (the oscillating tunnel and the wave boundary layer cases).

The particles are initially released from $x = 0$ in the form of a line source (aligned with equal vertical intervals) at the instant $t = 0$. The initial y position of the $n\mathrm{th}$ particle, $y_{n0}$, will, for example for the case of the oscillating tunnel, be

$$\begin{aligned} y_{n0} =\delta _\mathrm{b} +\frac{2( {h-\delta _\mathrm{b} })}{N+1}n \end{aligned}$$

(23)

in which N is the total number of particles released. Once a particle is released, it will move in small discrete steps. Each step will take a time $\Delta t$, and the particle will change its horizontal (streamwise) and vertical (wall-normal) positions from x and y to $x+\Delta x$ and $y+\Delta y$, respectively. The time interval $\Delta t$ for each step is calculated by

$$\begin{aligned}&\Delta t=\frac{\ell _\mathrm{p} }{\sqrt{\overline{{v}'^2} } }\end{aligned}$$

(24)

$$\begin{aligned}&\ell _\mathrm{p} =b\cdot \ell \end{aligned}$$

(25)

Here, $\sqrt{\overline{{v}'^2}}$ is the standard deviation of the wall-normal velocity component. The statistical character of the wall-normal velocity fluctuations is represented by the Gaussian (normal) distribution whose mean value is zero ($\bar{v}=0)$ and variance is a function of the wall-normal distance for the steady flow case, $\left( \sqrt{\overline{v^{\prime 2}} } (y) \right)$, and a function of both the wall-normal distance and the phase for the oscillatory flow cases, $\left( \sqrt{\overline{v^{\prime 2}} } (\omega _\mathrm {w} \,t,y)\right)$ in Eq. (24) is a length scale of turbulence that quantifies the vertical motion of the particles and this quantity, too, is a function of the wall-normal distance for the steady flow case, $\ell _\mathrm{p} (y)$, and a function of both the wall-normal distance and the phase for the oscillatory flow cases, $\ell _\mathrm{p} (\omega _\mathrm{w} \,t,y)$. As seen in Eq. (25), the scales $\ell$ and $\ell _\mathrm{p}$ are linked by the coefficient b, which is a constant that serves as a tuning parameter for the length scale. For the present model, the tuning was performed with the values of the dispersion coefficients of steady flow case obtained by the analytical models (Elder 1959; Sumer 1974). The value of b that makes the model results match the analytical values best came out to be 7.75. The quantities $\ell$ and $\sqrt{\overline{{v}'^2} }$ are calculated from the flow model through Eqs. (17) and (18), respectively.

The vertical displacement of the particle, $\Delta y$, is calculated from

$$\begin{aligned} \Delta y=( {{v}'-w_\mathrm{s} })\cdot \Delta t \end{aligned}$$

(26)

For neutrally buoyant particles, $w_\mathrm{s}$ is obviously set equal to zero. Here, $v'$ is the vertical fluctuation velocity calculated by

$$\begin{aligned} {v}'=a_\mathrm{r} \cdot \sqrt{\overline{{v}'^2} } ( {\omega _\mathrm{w} \,t,y}) \end{aligned}$$

(27)

where the coefficient $a_\mathrm{r}$ is a random variable whose probability density function satisfies the normal distribution with a mean equal to 0 and standard deviation of 1 (i.e. a standard normal variable). Thus, $v'$ can either be positive or negative, which will also affect the direction of the particle’s vertical displacement. It is worth mentioning that, for the sake of physical consistency, the value of $\Delta y$ was not allowed to be greater than 2h.

For the particles to be maintained in the main body of the flow, the following conditions are employed:

$$\begin{aligned} y_2 =\left\{ {\begin{array}{l} y+\Delta y , \quad \text{ when }\quad \delta _\mathrm{b} <y+\Delta y<( {2h-\delta _\mathrm{b} }) \\ 2( {2h-\delta _\mathrm{b} })-( {y+\Delta y}) , \quad \text{ when }\quad y+\Delta y\ge ( {2h-\delta _\mathrm{b} }) \\ \max \left\{ {\delta _\mathrm{b},\quad \left[ {2\delta _\mathrm{b} -( {y+\Delta y})-w_\mathrm{s} \cdot \Delta t} \right] } \right\} ,\\ \qquad \text{ when }\quad y+\Delta y\le \delta _\mathrm{b} \\ \end{array}} \right. \end{aligned}$$

(28)

in which $y_{2}$ is the new position of the particle in the vertical. Similar conditions are applied for the case of the wave boundary layer, to ensure that the particle is maintained in the flow all the time.

As the particle travels up and down in the flow domain, it will be carried by the flow in the streamwise direction. The Eulerian streamwise velocity of the flow, $\bar{u}(\omega _\mathrm{w} \,t,y)$, varies as a function of wall-normal distance and also the phase. A particle released in such a flow field will experience a time dependent Lagrangian velocity, U, as it is convected downstream. The horizontal distance travelled by the particle during the time interval $\Delta t$ is

$$\begin{aligned} \Delta x=\int \limits _t^{t+\Delta t} {U\mathrm{d}t} \end{aligned}$$

(29)

For the steady flow case, Eq. (29) can easily be expressed in terms of Eulerian velocity:

$$\begin{aligned} \Delta x=\frac{1}{( {y_2 -y})}\int \limits _y^{y_2 } {\bar{u}( y)\,\mathrm{d}y} \cdot \Delta t \end{aligned}$$

(30)

For the oscillatory-flow case, however, this velocity is calculated by the following numerical integration:

$$\begin{aligned} \Delta x= & {} \int \limits _t^{t+\Delta t} {U\mathrm{d}t} \nonumber \\\approx & {} \underbrace{\frac{1}{M}\sum \limits _{m=0}^M {\bar{u}\left( {\omega _\mathrm{w} \left( t+\frac{m}{M}\Delta t\right) ,y+\frac{m}{M}\Delta y}\right) } }_{\hbox {average }\hbox {Lagrangian }\hbox {velocity }\hbox {of }\hbox {the}\hbox { particle}}\cdot \Delta t \end{aligned}$$

(31)

in which $\bar{u}( {\omega _\mathrm{w} \,t,y})$ is taken from the flow-model calculations and M is an integer of resolution that holds the following condition:

$$\begin{aligned} M\ge \max \left\{ {\frac{1000\cdot \Delta t}{T_\mathrm{w} },\, \max \left[ {1000\frac{\Delta y}{h},\, \frac{1}{5}\cdot \frac{\Delta yUf}{\nu }} \right] } \right\} \end{aligned}$$

(32)

Here, $T_\mathrm{w}$ is the wave period equal to $T_\mathrm{w} =\frac{2\pi }{\omega _\mathrm{w} }$. If the particle is bouncing from any of the flow boundaries, then Eq. (31) is evaluated separately from the particle initial position to the boundary and from the boundary to the particle final position.

In the oscillatory flow cases, the value of $\Delta t$ needs to be limited with a fraction of the wave period, $T_\mathrm{w}$, considering that the timescale of turbulent motion should be significantly smaller than the wave period. In other words, the vertical travel of the particle cannot be sustained with a constant vertical velocity for long durations comparable with the wave period, given that the statistical properties of the turbulence characteristics are a function of the phase (i.e. $\omega _\mathrm{w}t)$. For example, for each $T_\mathrm{w}$/4 time interval, the free stream velocity rises from minimum (0) to maximum ($U_\mathrm{m})$ along with the bed shear stress and associated turbulence characteristics. Thus, one would expect that $\Delta t$ cannot exceed, say, a quarter of the wave period. The comparison of the numerical model results with that of the experiments showed that limiting the travel time to be $\Delta t \le$ 0.15 $T_\mathrm{w}$ gives a better match with the experiments. We further note that the difference in the dispersion coefficient caused by the selection of limiting value as 0.15 $T_\mathrm{w}$ and 0.25 $T_\mathrm{w}$ is rather small. It should also be mentioned that when $\frac{a}{h}\ge 20$, the effect of limiting $\Delta t$ on dispersion coefficient vanishes.

The calculation procedure implemented in MATLAB environment is summarised as follows:

• 1. Run the flow model for the designated flow conditions and obtain $\bar{u}( {\omega _\mathrm{w} \,t,y})$, $k( {\omega _\mathrm{w} \,t,y})$ and $\ell ( {\omega _\mathrm{w} \,t,y})$ for the lower half of the oscillating tunnel

• 2. Calculate $\sqrt{\overline{v'^2} } ( {\omega _\mathrm{w} \,t,y})$ and $\ell _\mathrm{p} ( {\omega _\mathrm{w} \,t,y})$ matrices through Eqs. (17) and (25), respectively

• 3. Mirror the matrices in step 2, to include also the upper half of the tunnel

• 4. Line-up N particles with the designated $\beta$ at $x=0$ and $t=0$ along the vertical with intervals described by Eq. (23)

• 5. Release the N particles one by one and run the dispersion simulation for each particle until $t=t_\mathrm{end}$.

  • (a) For the given phase ($\omega _{w }t)$ and vertical location (y) of the particle, pick up the values of $\ell _\mathrm{p}$ and $\sqrt{\overline{v'^2} }$ from the relevant matrices and calculate $\Delta t$ via Eq. (24)

  • (b) Check if $\Delta t \le$ 0.15$T_\mathrm{w}$

  • (c) Calculate the random $v'$ and $\Delta y$ from Eq. (27) and (26), respectively

  • (d) Calculate $y_{2}$ from Eq. (28)

  • (e) Calculate $\Delta x$ from Eq. (31)

  • (f) Update t, x and y as $t =t +\Delta t$, $x =x +\Delta x$ and $y =y_{2}$

  • (g) If $t < t_\mathrm{end}$, continue the simulation with repeating the steps a–f

  • (h) If $t > t_\mathrm{end}$, discard the last vertical travel of the particle and arrange a final travel by designating $\Delta t=t_\mathrm{end}$ – t.

  • (i) Save the path of the particle and go to step 5 for a new particle

• 6. After the last ($N\mathrm{th}$) particle completes its path, stop the dispersion simulation and calculate the relevant particle statistics.

4 Experiments

The experiments were conducted in the oscillating U-tunnel facility located at the Hydraulics Laboratory of Technical University of Denmark. A detailed description of the tunnel can be found in, e.g. Jensen et al. (1989). The longitudinal and cross-sectional sketches of the U-tunnel are given in Fig. 3. The walls of the U-tunnel were smooth. A vertical sheet of light by means of LED spots, 3 cm in thickness at the top of the tunnel, and 4 cm at the bottom of the tunnel (Fig. 3), was established to illuminate the particles travelling in the tunnel. The length of the light sheet was 2 m. The particle trajectories were videotaped with a camera (Fig. 3) through the transparent side wall of the tunnel at a rate of 50 frames per second. The camera angle was calibrated such that the distortion introduced by the lens was corrected.

Fig. 3

Longitudinal and cross-sectional sketch of the oscillating U-tunnel and experimental setup

A spherical particle, printed in plastic using a 3D printer, was used in the experiments. The particle size, the specific gravity and the settling velocity of the particles were $d=2.9$ mm, $s=1.19$ and $w_\mathrm{s}=9.8$ cm/s, respectively.

The free stream velocity was measured by a Dantec Dynamics FiberFlow laser Doppler anemometer (LDA) at 14.5 cm away from the bed (i.e. at the centreline of the tunnel). The LDA signal also served as a reference signal for monitoring the phase of the flow during particle tracking experiments. In addition to free stream velocity, velocity profiles were determined by means of LDA measurements conducted at 30 positions across the half thickness of the oscillating tunnel. At each location, velocity measurements were taken for 50 wave periods and ensemble averaging was conducted to obtain the Reynolds averaged velocity as a function of the phase, $\bar{u}({\omega _\mathrm{w} t,y})$. These measurements are not included here for reasons of space.

For the determination of the friction velocity, $U_\mathrm{fm}$, direct measurements of bed shear stress were conducted with a Dantec Dynamics 55R46 Hot-film probe. These measurements were also conducted for 50 wave periods to obtain the ensemble-averaged values.

Three sets of experiments were conducted, with Re $_\mathrm{w}= 6$ $\times$ 10$^{5}$, 1.5 $\times$ 10$^{6}$ and 3.3 $\times$ 10$^{6}$, respectively. The period of the oscillatory flow was kept constant as $T_\mathrm{w}=9.72$ s (the natural oscillation period of the tunnel). The flow properties and the particle properties are summarised in Table 1. In this table, $\theta$ is the Shields parameter. The boundary layer thicknesses ($\delta$) given in Table 1 were calculated by means of the flow model. Note that these $\delta$ values are in good agreement when compared with the values obtained from the mean velocity profiles measured at $\omega _\mathrm{w}t=90^{\circ }$.

Table 1 The experiments conducted in the oscillating tunnel

In the experiments, the particle was released into the flow domain multiple times (Table 1, last column). Then, the particle trajectory was videotaped, and subsequently the mean particle position and the variance along with the longitudinal dispersion coefficient were calculated. Similar particle tracking experiments were previously conducted by Sumer and Oguz (1978) and Sumer and Deigaard (1981) in an open channel flow.

One of the limitations of the experimental setup was that, it was not possible to track the particles too long before they disappear out of the light sheet. Typical durations of monitoring were about one wave period. It should also be noted that the particles might be lost for some frames in the video, but then found again and subsequently tracked. The last three columns of Table 1 refer to the minimum, maximum and mean number of particles that could be tracked at each frame during one complete wave period. Further details about the experiments can be found in Steffensen (2014).

5 Results

5.1 Validation of the flow model

The flow model used in this study, MatRANS, will be validated against different steady and unsteady flow cases (Fuhrman et al. 2013). For validating the results of the present modified version of the model, it is compared with the oscillatory flow measurements of the literature.

Figures 4, 5 and 6 present the comparison of model results with the experimental data of Jensen et al. (1989) for the oscillatory flow case of $Re_\mathrm{w} =\frac{U_\mathrm{m} ^2}{\omega \nu }=6.2\times 10^6$. Note that in Figs. 5 and 6, the y-axis is non-dimensionalised by the horizontal amplitude of the free stream motion, $a={U_\mathrm{m} }/ \omega _\mathrm{w}$. In Fig. 4, the variation of the friction velocity, $U_\mathrm{f}$, during one period is given, whereas Reynolds-averaged velocity profiles for different wave phases are shown in Fig. 5. As can be seen there is a good agreement between the model results and the experimental data. The turbulence kinetic energy profiles for different wave phases are given in Fig. 6. It should be noted that the turbulence kinetic energy values from Jensen et al. (1989) were not directly measured in three dimensions, but approximated by $k\approx 0.65( {\overline{u'^2} +\overline{v'^2} })$ (Justesen 1991). Generally, the model can capture the turbulent kinetic energy profiles well. Yet, the turbulence in deceleration phases (where the pressure gradient is adverse) seems to be relatively underestimated by the model. It is worth mentioning that comparison against the steady flow data from Fuhrman et al. (2010) has likewise yielded a very good agreement both for mean flow and turbulence properties, though these results are not included herein for brevity.

Fig. 4

Comparison of model results for oscillatory flow with the measurements of Jensen et al. (1989) ($Re_\mathrm{w}=6.2\times 10^{6}$). Variation of friction velocity during one wave period ($U_\mathrm{f}/U_\mathrm{m}$ vs. $\omega _\mathrm{w}t$)

Fig. 5

Comparison of model results for oscillatory flow with the measurements of Jensen et al. (1989) ($Re_\mathrm{w}=6.2\times 10^{6}$). Left and right panel show the velocity profiles in acceleration and deceleration stages of the flow for each 15$^{\circ }$ increments of $\omega _\mathrm{w}t$. Symbols same as in Fig. 4

Fig. 6

Comparison of model results for oscillatory flow with the measurements of Jensen et al. (1989) ($Re_\mathrm{w}=6.2\times 10^{6}$). Turbulent kinetic energy profiles for each 15$^{\circ }$ increment of $\omega _\mathrm{w}t$. Symbols same as in Fig. 4

Similar to the flow model used in this study, recently Tang and Lin (2015) developed a RANS model established on the k–$\omega$ baseline formulation, which can simulate the turbulent as well as transitional boundary layer flows. They also compared their findings with the experimental data of Jensen et al. (1989), yielding a good agreement.

5.2 Longitudinal dispersion in steady flow

Throughout the present section, the results will be presented in non-dimensional forms with

$$\begin{aligned} X=\frac{x}{h},\quad Y=\frac{y}{h},\quad T=\frac{tU_\mathrm{f} }{h},\quad U=\frac{u}{U_\mathrm{f} } \end{aligned}$$

(33)

Before the results of dispersion model were acquired, a sensitivity analysis was conducted with different number of particles to optimize the number of particles, N, to be used in the simulations. This is obviously a trade-off between computational effort and higher accuracy. The results of these preliminary runs are given in Fig. 7. As can be seen, when the number of particles reaches roughly 10,000, the mean and variance of particle positions become practically stable. Therefore, $N=10{,}000$ was adopted for the dispersion simulations.

Fig. 7

Sensitivity of mean and variance of particle positions with respect to the number of particles released into the flow, N, in the random walk simulation

The flow parameters used in the simulations are as follows: The half thickness of the tunnel, $h=14.5$ cm, the mean flow velocity, V = 2 m/s, the friction velocity, $U_\mathrm{f}$ = 8 cm/s. The simulations were conducted over a time period of $0\le T\le 20$, which corresponds to $0\le t\le 11$ s with the present flow parameters. The vertical and horizontal particle positions (snapshot of particle locations) were recorded by interpolation at frequent time intervals (with 0.1 in terms of non-dimensional time). The simulations were conducted for different settling velocities, covering $0\le \beta \le 1$.

Figure 8 shows the time series of mean particle positions for neutrally buoyant particles. On this figure, the displacement obtained by integration of sectional average velocity (V) is also shown. As can be seen, the mean particle velocity (i.e. velocity of the particle cloud) is almost identical to the mean cross-sectional velocity. In other words, Eulerian velocity and Lagrangian velocity are equal, as expected in a steady ($\partial {/}{\partial t}=0)$ and uniform ($\partial {/}{\partial x}=0)$ flow environment.

Fig. 8

Non-dimensional mean particle X positions versus non-dimensional time, steady flow, neutrally buoyant particles ($\beta =0$)

The concentration of particles, c(x, y, t), can be defined as the number of particles per unit area. The non-dimensional particle concentration, C, can then be defined as:

$$\begin{aligned} C( {X,Y,T})=c\frac{h^2}{N} \end{aligned}$$

(34)

Using Aris’ moment method, Sumer (1974) found that the zeroth moment of concentration, defined by

$$\begin{aligned} C_0 =\int \limits _{-\infty }^\infty {C\mathrm{d}\xi } \end{aligned}$$

(35)

in which

$$\begin{aligned} \xi =\frac{x-\overline{U} \;t}{h} \end{aligned}$$

(36)

is, for large times, given as:

$$\begin{aligned} C_0 ( {Y,T})=\frac{\sin \pi \beta }{\pi \beta }\left( {\frac{1-Y}{Y}}\right) ^\beta ,\quad T>7.5 \end{aligned}$$

(37)

For neutrally buoyant particles $C_0 ( {Y,T})$, for large times, becomes unity. In Fig. 9, the zeroth concentration moment calculated from the present simulation (corresponding to $T=20$) is plotted. The concentration value is not 1 but 0.5 since the thickness of the flow domain is 2h. As can be seen, model results are consistent with the analytical expression.

Fig. 9

Zeroth moment of concentration ($C_0$) for large times compared with analytical solution, steady flow, neutrally buoyant particles ($\beta =0$)

Figure 10 shows the variance of the particle X positions plotted as function of time. Initially the variance increases with time like $\overline{(X-\overline{X} )^2} \propto t^2$, but after approximately $T \approx$ 3, it increases linearly, $\overline{(X-\overline{X} )^2} \propto t$, in agreement with the theory.

Fig. 10

Change of non-dimensional particles variance with time, steady flow, neutrally buoyant particles ($\beta =0$)

When the same simulation is conducted for heavy particles with $\beta =0.3$, the picture changes significantly. Figure 11 shows the particle mean X positions with time. It can be seen that the particle cloud travels slower than the mean Eulerian flow velocity and it is delayed more and more with time.

Fig. 11

Non-dimensional mean particle X positions versus non-dimensional time, steady flow, heavy particles ($\beta =0.3$)

In Fig. 12, the zeroth concentration moment at $T=20$ is given for a heavy particle with $\beta =0.3$ along with the analytical profile from Sumer (1974), Eq. (37) above. A good agreement can be seen.

Fig. 12

Zeroth moment of concentration ($C_{0})$ for large times compared with analytical solution, steady flow, heavy particles ($\beta$ = 0.3)

The numerical simulation has been repeated for particles with different settling velocities. The results of these simulations are presented in Fig. 13 in terms of the ratio of dispersion coefficient of heavy particles to that of neutrally buoyant particles (${D_1 ( \beta )}/{D_1 ( 0)})$. The dispersion coefficient is obtained from $D_1 =\frac{1}{2}\frac{\mathrm{d}\overline{(X-\overline{X} )^2} }{\mathrm{d}T}$ for large times (Eq. 4), i.e. half of the slope of the straight line in the variance plot (in Fig. 10). Incidentally, the numerical value of the dispersion coefficient for neutrally buoyant particles was found to be $D_1 ( 0)=5.96$, and this value is in good agreement with that obtained analytically by Elder (1959), $D_1 ( 0)=5.86$. Along with the model results, the analytical solution of Sumer (1974) and numerical solution of Sayre (1968) are also plotted in Fig. 13. As the particles become heavy, the dispersion coefficient increases significantly, since they start to be transported close to bed where the velocity gradient is steeper. The model results generally agree with the analytical solution of Sumer (1974). The slight discrepancies in Fig. 13 may be attributed to the way in which the reflective boundary condition at the bottom is formulated in these three studies.

Fig. 13

Ratio of dispersion coefficient for heavy particles to that of neutrally buoyant particles versus the non-dimensional settling velocity ($\beta$)

The simulation results summarised in the preceding appear to be in good agreement with the existing information. This has provided confidence in the use of the present methods for the oscillatory flow cases studied in the present research.

5.3 Longitudinal dispersion in oscillating tunnel

The flow parameters used in the oscillating-tunnel simulations are summarised in Table 2. The friction velocity ($U_\mathrm{fm})$ and boundary layer thickness ($\delta$) values given in the table were obtained from the results of flow model runs.

Table 2 The summary of oscillating-tunnel dispersion simulations

Figures 14 and 15 display the particle trajectories for two cases, one for the case of the neutrally buoyant particle ($\beta =0$), and the other for a heavy particle with $\beta =0.3$. The particles are released from the middle of the flow domain and tracked for five wave periods. From Fig. 14, it is seen that the particle is moving dominantly by the free stream velocity when it is away from the walls. Yet, the turbulence and steep velocity gradient make the particle path deviate from the free stream whenever the particle comes close to the walls. This particle spends its time more or less uniformly across the thickness of the flow domain.

Fig. 14

Non-dimensional path of a neutrally buoyant ($\beta =0$) particle released in oscillatory flow for five consequent wave periods, $a/h=21.4$ and $Re_\mathrm{w}=6.2\times 10^{6}$

Fig. 15

Non-dimensional path of a heavy ($\beta =0.3$) particle released in oscillatory flow for five consequent wave periods, $a/h=21.4$ and $Re_\mathrm{w}=6.2\times 10^{6}$

Figure 15 shows the path of the heavy particle (with $\beta =0.3$). The aspects mentioned in the preceding paragraph can more vigorously be seen in this figure. Since this particle is a heavy particle, it spends more time closer to the bottom, and thus experiences relatively stronger shear and turbulence. Hence, one would expect that heavy particles would be subjected to higher dispersion, just as in the case of steady flow.

For direct comparison, Fig. 16 presents the velocity gradient experienced by the two particles throughout their paths. In the oscillatory flow case, the mean of the non-dimensional velocity gradients for the neutrally buoyant particle ($\beta =0$) and heavy particle ($\beta =0.3$) is around 5 and 36, respectively. The difference is around a factor 7 and shows how the settling velocity of the particle affects the dispersion process significantly.

Fig. 16

The non-dimensional velocity gradients ($\frac{{\partial U}}{{\partial Y}}$) that the particles in Figs. 14 and 15 experience along their paths

Figure 17 presents the variance in non-dimensional form ($\frac{\,\overline{( {x-\bar{x}})^2} \,}{\delta ^2})$ plotted versus non-dimensional time ($\frac{tU_\mathrm{fm} }{\delta })$ for the neutrally buoyant particle. In this flow case, one wave period corresponds to 16.33 in terms of non-dimensional time. Contrary to the steady flow case, the variance does not increase with a constant rate, but it exhibits an increase at a constant rate with a superimposed periodic variation. Note that this behaviour was also observed by other researchers as well, e.g. Mazumder and Paul (2012). The non-dimensional dispersion coefficient of neutrally buoyant particles ($\frac{D_1 }{\delta U_\mathrm{fm} })$ is found to be around 1.5, a factor 4 smaller than the steady flow case.

Fig. 17

Non-dimensional variance of particle position ($\frac{{\overline{{{\left( {x - \overline{x} } \right) }^2}} }}{{{\delta ^2}}}$) as a function of non-dimensional time ($\frac{{{{{}^{tU}}_\mathrm{fm}}}}{\delta }$) for the oscillatory tunnel case, black line phase-averaged variance. Neutrally buoyant particles ($\beta =0$), $a/h=21.4$, and $Re_\mathrm{w}=6.2 \times 10^{6}$

To demonstrate the effect of a/h on the dispersion process, the following parametric study was conducted. First, the half thickness of the flow domain (h) was kept constant at 14.5 cm and a/h was changed by changing the wave Reynolds number (Re $_\mathrm{w})$ of the oscillatory flow from 2.7 $\times$ 10$^{5}$ to 9.9 $\times$ 10$^{7}$ (Group no. 1 experiments in Table 2). Secondly, Re $_\mathrm{w}$ was kept constant at 6.2 $\times$ 10$^{6}$ and h was changed covering a wide range of a/h (Group no. 2–10 experiments in Table 2). For each case, the non-dimensional dispersion coefficient was calculated as a function of a/h. The results are plotted in Fig. 18. As can clearly be seen here, both the above results collapse generally onto the same curve, which evidently shows that the non-dimensional dispersion coefficient is not dependent on Re $_\mathrm{w}$, but it is a function of a/h. The results presented by Fig. 18 further show that the dispersion coefficient increases with increasing a/h and attains asymptotically to the dispersion coefficient for steady open-channel-flow value. This is expected because, for large values of a/h (very small values of the tunnel height), the parameter a/h will not play any significant role in the oscillatory boundary layer process, as the boundary layer will extend across the entire tunnel height and, therefore, the tunnel height is no longer a controlling parameter for the boundary layer process. The decrease in the dispersion coefficient with decreasing a/h, on the other hand, can be explained as follows. The smaller the value of the parameter a/h, the larger is the tunnel height and, therefore, the particle should experience smaller velocity gradients during its travel in the flow, implying that the dispersion coefficient should decrease with decreasing a/h.

Fig. 18

Variation of non-dimensional dispersion coefficient for neutrally buoyant particles ($\beta =0$) with non-dimensional amplitude (a / h). Circles $Re_\mathrm{w} =6.2\times 10^{6}$, Triangles variable $Re_\mathrm{w}$

The dispersion of heavy particles with $\beta =0.3$ was simulated for the oscillatory tunnel case depicted in Figs. 4, 5 and 6 (i.e. Re $_\mathrm{w}$ = 6.2 $\times$ 10$^{6}$ and a/$h=21.4$). The time variation of variance of particle positions in non-dimensional form is given in Fig. 19. The behaviour is similar to that of the neutrally buoyant particle, but the dispersion is stronger, given that the particle spends more time close to the bed. The dispersion coefficient came out to be around 4, more than a factor two higher than the neutrally buoyant case. Recall that, in the case of the steady flow, the dispersion coefficient for $\beta =0.3$ was found to be approximately 10, a significantly higher value.

Fig. 19

Non-dimensional variance of particle position ($\frac{{\overline{{{\left( {x - \overline{x} } \right) }^2}} }}{{{\delta ^2}}}$) as a function of non-dimensional time ($\frac{{{{{}^{tU}}_\mathrm{fm}}}}{\delta }$) for oscillatory tunnel case, black line phase-averaged variance. Heavy particles ($\beta =0.3$), $a/h=21.4$, and $Re_\mathrm{w}=6.2 \times 10^{6}$

The zeroth concentration moments (i.e. concentration profiles) of heavy particles are given in Fig. 20. On this figure, the concentration profile of the same particle in steady flow is also shown. As can be seen, the particles tend to come closer to the bed in the oscillatory flow case as compared to the steady flow. Turbulence cannot counterbalance the settling tendency of particles as pronounced as in the steady flow case. As a result, fewer particles can be carried to the upper regions of the flow domain.

Fig. 20

The zeroth concentration moment of heavy particles ($\beta =0.3$) for oscillatory flow compared to steady flow. $a/h=21.4$, and $Re_\mathrm{w}=6.2 \times 10^{6}$

With the purpose of clearly demonstrating the non-dimensional dispersion coefficient as a function of $\beta$ and a/h, the model was run for a wide range of the two governing parameters (Table 2). Results are given as a family of curves in Fig. 21. Note that the $\beta =0$ curve in Fig. 21 is identical with the circles in Fig. 18.

Fig. 21

Variation of non-dimensional dispersion coefficient with non-dimensional amplitude (a / h) for different non-dimensional settling velocities ($\beta$)

Figure 21 suggests that the dependency of the dispersion coefficient on non-dimensional amplitude gets weaker as the particles get heavier (i.e. for large $\beta$ values). Likewise, for short amplitudes (a/40) the effect of $\beta$ on the dispersion coefficient is more dominant than longer amplitudes.

In the study of Mazumder and Paul (2012), the bed condition defined for the particles was different than the present study, such that a settlement (bed absorbency) and re-entrainment procedure was employed. Furthermore, their study does not clearly state an increase in the dispersion coefficient with increasing settling velocity. Thus, the two studies are not directly comparable. Nevertheless, the time variation of the variance found by Mazumder and Paul (2012) qualitatively agrees well with the present results (Figs. 17, 19).

Ng (2004) carried out an analytical study to find out the dispersion coefficients of neutrally buoyant particles in oscillatory flow. They used the free stream velocity and amplitude of the horizontal motion to non-dimensionalise the dispersion coefficient. When converted to the terms of the present study, their findings for a/$h=32$ and 96 come out to be $\frac{D_1 }{\delta U_\mathrm{fm} }\approx$ 3.8 and 5, respectively. In the present model, those two values are found as 3.1 and 4.2. When compared, these findings are not radically diverse, given that the methodology followed by Ng (2004) and the present study is quite different.

5.4 Comparison with the oscillating tunnel experiments

As given in Sect. 4, three oscillating tunnel experiments were conducted (Table 1). The non-dimensional settling velocity of the particles in these experiments was relatively high (namely $\beta =3.6, 4.8$ and 7.2). Nevertheless, the bottom conditions were such that the particles did not deposit on the bottom when they came close to the bottom, and they were maintained in the flow all the time. These experiments were simulated in the numerical model and obtained results are compared with those of the experiments.

Figures 22, 23 and 24 compare the results of the conducted experiments (Table 1) and the numerical simulation for the phase-averaged zeroth moment of concentration profiles ($C_{0})$ as a function of non-dimensional vertical distance (Y).

Fig. 22

Phase-averaged zeroth moment of concentration ($C_{0})$ profile, $\beta =3.6$, a/$h=16.4$,$Re_\mathrm{w} =3.3\times 10^6$

Fig. 23

Phase-averaged zeroth moment of concentration ($C_0$) profile, $\beta =4.8$, $a/h=10.7$, and $Re_\mathrm{w}=1.5 \times 10^{6}$

Fig. 24

Phase-averaged zeroth moment of concentration ($C_{0}$) profile, $\beta =7.2$, $a/h=6.7$, and $Re_\mathrm{w}=60 \times 10^{5}$

As is seen from the figures, the numerical simulations exhibit a general agreement with the experimental results for all the cases. For $\beta =3.6$, the numerical model estimates the concentration at the upper regions slightly lower compared to the experimental results. As for $\beta =7.2$, the estimate of the model is fairly more uniform while the experiment gave a more peaked result. It should, furthermore, be noted that this case also corresponds to the lowest wave Reynolds number (6  $\times$  10$^{5})$, which is in the transitional regime. A probable explanation why the model seems to perform better for $\beta =4.8$ is that, the number of particles captured in the experiments with $\beta =4.8$ is almost two times larger than the other two experimental sets (Table 1). As the number of particles is increased, a better resolution of vertical concentration profile of particles could be obtained in the experiments, and thus a better comparison with the numerical results was maintained. These results suggest that the numerical model can capture the suspension behaviour of heavy particles reasonably well.

The variation of the mean particle positions ($\overline{X} )$ recorded as a function of phase ($\omega _\mathrm{w}t)$ during the three experiments is shown in Figs. 25, 26 and 27, respectively, along with the results of the numerical simulation. Recall that $\omega _\mathrm{w}t =0$ corresponds to the zero up-crossing point in the free stream velocity. The lower panes of these figures show $X_0 =-a\cos ( {\omega _\mathrm{w} t})$, the non-dimensional displacement of the fluid in the free stream as a reference signal. In these figures, the phase-averaged positions of particles in numerical and experimental results are matched at $X=0$, since the experiments start from an arbitrary position.

Fig. 25

Mean particle position ($\overline{X} )$ as a function of phase ($\omega _\mathrm{w}t)$, $\beta =3.6$, a/$h=16.4$, $Re_\mathrm{w} =3.3\times 10^6$

Fig. 26

Mean particle position ($\overline{X} )$ as a function of phase ($\omega _\mathrm{w}t)$, $\beta =4.8$, a/$h=10.7$, $Re_\mathrm{w} =1.5\times 10^6$

Fig. 27

Mean particle position ($\overline{X} )$ as a function of ($\omega _\mathrm{w}t)$, $\beta =7.2$, a/$h=6.7$, $Re_\mathrm{w} =6.0\times 10^5$

For all the three tested cases, the displacement of the particles in the oscillating tunnel experiments is approximately 30 % lower compared to the numerical simulation results. This is somehow expected, since the particles in the oscillating tunnel have a finite size and a higher specific mass than water, and thus they experience a certain lag each time the surrounding fluid is accelerated owing to their inertia.

The variance of particle positions in non-dimensional form ($\frac{\,\overline{( {x-\bar{x}})^2} \,}{\delta ^2})$ is presented as a function of the non-dimensional time ($\frac{tU_\mathrm{fm} }{\delta })$ in Figs. 28a, 29a and 30a for the three cases, such that the half of the slope of the straight line in the variance plot gives non-dimensional dispersion ($\frac{D_1 }{\delta U_\mathrm{fm} })$ coefficient. Figures 28b, 29b, and 30b display the numerical simulation results for comparison.

Fig. 28

Non-dimensional variance of particle position ($\frac{\overline{( {x-\bar{x}})^2}}{\delta ^2})$ as a function of non-dimensional time ($\frac{tU_\mathrm{fm} }{\delta })$, a Experiments (black line phase-averaged variance), b Numerical simulations. $\beta =3.6$, a/$h=16.4$, $Re_\mathrm{w} =3.3\times 10^6$

Fig. 29

Non-dimensional variance of particle position ($\frac{\overline{( {x-\bar{x}})^2} }{\delta ^2})$ as a function of non-dimensional time ($\frac{tU_\mathrm{fm} }{\delta })$, a experiments (black line phase-averaged variance), b numerical simulations. $\beta =4.8$, $a/h=10.7$, $Re_\mathrm{w} =1.5\times 10^6$

Fig. 30

Non-dimensional variance of particle position ($\frac{\overline{( {x-\bar{x}})^2} }{\delta ^2})$ as a function of non-dimensional time ($\frac{tU_\mathrm{fm} }{\delta })$, a experiments (black line phase-averaged variance), b numerical simulations. $\beta =7.2$, $a/h=6.7$, $Re_\mathrm{w} =6.0\times 10^5$

It is seen that the dispersion behaviour in the experiments basically composed of an increasing trend of variance superimposed with some periodic variations, as depicted in Sect. 5.3. The time variation of the variance is not captured for large times due to the small sample size. Although the numerical model generally captures the constant increasing trend fairly well, the periodic variation amplitudes calculated by the model are significantly lower than those of the experiments, presumably attributed to the settling velocity of the particles being increased beyond the suspension threshold (i.e. $\frac{w_\mathrm{s}}{U_\mathrm{fm} }>>1)$ .

The results of non-dimensional dispersion coefficients found from the half of the slope of the trend line of the experimental data in Figs. 28a, 29a and 30a are compared in Table 3. The numerical model result is in fairly good agreement with the experimental result for $\beta =3.6$ while it differs as much as more than 35 % for the other two cases. No clear explanation has been found for this discrepancy. Nevertheless, the experimental data suffer from the small sample size, as pointed out earlier. Another possible reason is that the dispersion model mainly deals with fine sediment which has low fall velocity and is kept in suspension in almost all times. As the sediment gets heavier (as the fall velocity increases), the particles presumably tend to be transported as bed sediment, at least during some phases of the flow, which causes a discrepancy between the model and experimental results.

Table 3 Comparison of model results with experiments

6 Longitudinal dispersion in wave boundary layer: practical application

In this section, a practical example will be given on how the presented numerical model can be used to simulate the longitudinal dispersion of heavy particles in wave boundary layers. The following case is considered:

Excavated material (fine sand with $d_{50}=0.10$ mm) from a dredging operation is disposed to a coastal area at a depth of 10 m. The settling velocity of the material is $w_\mathrm{s}=0.65$ cm/s. Some time after disposal, the material disperses under the action of storm waves with a root-mean square wave height of $H_\mathrm{rms}=2.6$ m and zero-crossing period of $T_\mathrm{z}=8$ s. This process will be simulated as follows by means of the presented numerical model.

First, the flow model is run for the determination of the mean flow and turbulence characteristics developing under the wave boundary layer. As the free stream velocity, the horizontal orbital velocity at the bed found from the linear wave theory, $U_\mathrm{m}=1.01$ m/s, was used. The model upper boundary (i.e. frictionless lid) was set to $h=1$ m. It can be shown by linear wave theory that the orbital velocity at the top of model upper boundary is practically equal to the one at the bed ($U_{m (y=1m)} =1.02$ m/s). Thus, the oscillatory flow model results are directly applicable for this wave boundary layer case.

The flow model yielded a friction velocity of $U_\mathrm{fm}=5.5$ cm/s and boundary layer thickness of ${\delta } =2.5$ cm, compatible with the ones that can be found from the classical wave boundary layer formulae (e.g. Fredsøe and Deigaard 1992, chp. 2). For these conditions, the non-dimensional settling velocity of the particles corresponds to $\beta =0.3$. With the output of the flow model, the random-walk dispersion model was run with $N=10{,}000$ particles. The particles are released from the bed and none of the particles were seen to travel higher than 0.32 m from the bed. Figure 31 displays the variance of particle position versus time for the first eight waves at the beginning of the storm. From the variance-versus-time data, the dispersion coefficient was calculated as $\frac{D_1 }{\delta U_\mathrm{fm} }=3.5$ or $D_{1}=4.8\times 10^{-3}$  m$^{2}$/s.

Fig. 31

Non-dimensional variance of particle position ($\frac{\overline{( {x-\bar{x}})^2} }{\delta ^2})$ as a function of non-dimensional time ($\frac{tU_\mathrm{fm} }{\delta })$ for the first eight waves (black line phase-averaged variance). Water depth 10 m, $\beta =0.3$, $H_\mathrm{rms}=2.6$ m, $T_\mathrm{z}=8$ s

To illustrate the spreading of the dumped sediment as a result of the dispersion caused by the storm waves, it is assumed that the initial layer of the dumped sediment (before the storm breaks) covered a width of $B_{0}=20$ m in a uniform manner with a small thickness compared to depth. If the storm lasts about eight hours, the increase of the variance of the sediment longitudinal positions can be calculated by use of Eq. (4):

$$\begin{aligned}&\left( {\overline{( {x-\bar{x}})^2} }\right) _1 -\left( {\overline{( {x-\bar{x}})^2} }\right) _0 =2D_1 \,t_\mathrm{storm}\nonumber \\&\quad =2\times (4.8\times 10^{-3})\times (8\times 3600)\cong 277m^{2} \end{aligned}$$

(38)

The variance of initial positions of the disposed sediment layer, which is composed of N grains uniformly distributed between $-\frac{B_0 }{2}$ and $\frac{B_0 }{2}$, can be calculated as follows:

$$\begin{aligned} \overline{( {x-\bar{x}})^2} _0= & {} \frac{1}{N}\sum \limits _{n=1}^N {( {x-\bar{x}})^2}\nonumber \\= & {} \frac{1}{N}\sum \limits _{n=1}^N {\left( {n\frac{B_0}{N}-\frac{B_0 }{2}}\right) ^2} \approx 33 m^{2} \end{aligned}$$

(39)

The result of Eq. (39) converges to a fixed value for 100. Hence the variance of the longitudinal positions of the dumped sediment layer becomes $( {\overline{( {x-\bar{x}})^2} })_1 =310$ m$^{2}$ after the storm disperses it. With a back calculation employing Eq. (39), the equivalent width of the sediment layer after the storm can be found as $B_{1} \approx$ 61 m. This example has been illustrated in Fig. 32, obtained by the model simulation, which shows the initial and final concentration distribution (i.e. the probability density function) of the particles in the dumped sediment layer.

Fig. 32

Initial (upper panel) and final (lower panel) concentration distribution (pdf of particles) of the disposed sediment to a depth of 10 m, $\beta =0.3$, $H_\mathrm{rms}=2.6$ m, $T_\mathrm{z}=8$ s, $t_\mathrm{storm}=8$ h

7 Conclusions

In the present study, the dispersion process of suspended sediments in an oscillatory tunnel is studied with application on wave boundary layers. The problem is modelled by means of a Lagrangian numerical model that combines a k–$\omega$ RANS model and a random-walk numerical scheme. Model results for dispersion of neutrally buoyant and heavy particles under different oscillatory flow conditions were presented and discussed. Finally, an application example regarding the dispersion of fine sand dumped on seabed under storm waves was presented. The conclusions drawn from the present study can be summarised as follows:

• The results for steady flow case show that the model results are consistent with the existing information in the literature, and well represent the physics of the process for suspended particles unless the settling velocity is significantly high.

• In the case of oscillatory flow, the variance of particle positions does not increase monotonously as in steady flow case, but it exhibits some periodic oscillations accompanied by a continuous trend of dispersion.

• In the oscillatory flow case, the dispersion process is governed by the Rouse parameter as well as an additional non-dimensional parameter representing the amplitude of the motion. Generally, in oscillatory flow the dispersion coefficients are significantly less than that of the steady flow. The primary reason for this decrease appears to be the lower shear and turbulence experienced by the particles. This is demonstrated graphically and numerically in the study.

• For oscillatory flow, the dispersion coefficient increases with the settling velocity in a similar manner with the steady flow case. The reason for this increase is, again, that heavy particles moving closer to bed experience a higher shear (i.e. velocity gradient) and turbulence throughout their paths, which presumably disperses them farther and farther.

• When the amplitude of the oscillatory motion increases, dispersion coefficient of neutrally buoyant particles increases significantly and asymptotically attains to the value of the steady, open channel, flow case. As the particles get heavier (i.e. as the Rouse parameter increases), the dispersion coefficient under oscillatory flow becomes less dependent on the amplitude of the motion.

• The comparison of the present findings with the result of the experimental campaign, conducted in a U-tunnel facility, showed a reasonable agreement although the settling velocity values in the experimental study were beyond the suspension regime.

• Finally, the use of the present model was demonstrated with a practical case example in which the dispersion of fine sand under storm waves was studied.

The next stage of the present research shall be including provisions for the effect of steady streaming as well as co-existing current. Additionally, the settlement and re-suspension processes (see for example Sumer 1977) of the bed sediment may also be included in the model calculations.