# Modelling of sediment suspensions in oscillating grid turbulence

Hervé Michallet, Mathieu Mory

► To cite this version:

Hervé Michallet, Mathieu Mory. Modelling of sediment suspensions in oscillating grid turbulence. Fluid Dynamics Research, 2004, 35, pp.87-106. 10.1016/j.fluiddyn.2004.04.004. hal-00260290

HAL Id: hal-00260290

https://hal.science/hal-00260290v1

Submitted on 16 Nov 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Modelling of sediment suspensions in oscillating grid turbulence

H. Michalletª, M. Moryᵇ,\*

## Abstract

A $k-\epsilon$ model is used to describe the steady state of fine-grained sediments maintained in suspension by purely diffusive turbulence, as generated in oscillating grid turbulence experiments. The behaviour is shown to depend both on the bulk Rouse number $Rou_0$ and the product of the bulk Rouse number and the bulk Richardson number $Ri_0$. $Rou_0$ is built on oscillating grid parameters. For $Rou_0 < 0.01$, concentrated suspensions are observed with a homogeneous particle concentration in the suspension layer. An interface, called lutocline, separates the suspension layer from the clear water at a distance $z_m$ from the grid. The depth of the suspension layer is found to vary as $z_m/z_0 \propto (Ri_0 Rou_0)^{-1/4}$. For $Ri_0 Rou_0 \ll 1$ the decay of turbulence is affected by the particle concentration only in a region very close to the interface. In this case the flux Richardson number approaches the value of 1 near the interface. The lutocline is seen to vanish for large values of $Rou_0$. For $Rou_0 > 0.01$ the mean sediment concentration and turbulence decay simultaneously with increasing distance from the grid, and no sharp interface is observed.

PACS: 47.27 Eq; 47.27 Qb; 47.55 Hd; 47.55 Kf; 92.10 Wa; 92.10 Sx

Keywords: Turbulence simulation and modelling; Turbulent diffusion; Stratified flows; Multiphase and particle laden flows; Sediment transport; Coastal and estuarine oceanography

<footnote>ªLaboratoire des Ecoulements Géophysiques et Industriels (UJF-INPG-CNRS), BP53, 38041 Grenoble Cedex 9,
France</footnote>

<footnote>ᵇEcole Nationale Supérieure en Génie des Technologies Industrielles, Université de Pau et des Pays de l'Adour,
BP 7511, 64075 Pau, France</footnote>

<footnote>\* Corresponding author: Tel.: +33-05-59-40-78-02; fax: +33-05-59-40-78-01.
E-mail addresses: herve.michallet@hing.inpg.fr (H. Michallet), mathieu.mory@univ-pau.fr (M. Mory).</footnote>
# 1. Introduction

In many natural flows, the production of turbulence is due to the shear of the mean flow and is localised mainly at the bottom of the fluid layer. It is well known that turbulence can maintain high concentrations of particles in a steady-state suspension, independent of the way turbulence is produced. In particular, E and Hopfinger (1987) and Huppert et al. (1995) have shown that a steady two-layer system is formed in oscillating grid turbulence experiments, with a sharp interface separating a lower layer having an approximately constant sediment concentration and an upper layer containing almost no sediment. In such experiments, the turbulence produced by the grid oscillation is transported inside the tank by the sole effect of turbulence diffusion.

The former studies by E and Hopfinger (1987) and Huppert et al. (1995) were concerned with non-cohesive particles. The present work originated from oscillating grid turbulence experiments (Gratiot, 2000; Mory et al., 2002) performed using natural cohesive sediments experiments made of mud. Concentrated benthic suspensions (CBS) were also observed with a sharp interface. In the case of muddy sediments, very significant variations in the depth of the suspension layer were observed depending on the concentration in the suspension. These are the result of large variations in the settling velocity of muddy sediments. Experiments with cohesive sediment are particularly difficult to interpret because of the complexity of measuring the settling velocity of mud flocs. A striking property observed in cohesive sediment suspensions (Mory et al., 2002) was that the flux Richardson number, defined as the ratio of the buoyancy flux to the available turbulence, increases towards a value close to one in the vicinity of the lutocline. Although such estimates have also been obtained from numerical simulations of purely diffusive turbulence in stable stratified fluid (Briggs et al., 1998), it is generally considered that most of the mechanical energy input is dissipated viscously in stratified flows (Hopfinger and Linden, 1982) and, in particular, the flux Richardson number is below 0.25 in turbulent mean-shear flows (Ivey and Imberger, 1991).

In this paper, we investigate the conditions of occurrence of particle suspension layers in an oscillating grid tank experiment using a steady $k-\epsilon$ model. This model is still the simplest and cheapest for engineering applications, although more elaborate models such as Reynolds stress models (Straatman et al., 1998) and direct numerical simulation (Briggs et al., 1998) have also been used recently for diffusive turbulence. There was a long debate in the past as to whether it is appropriate to use a $k-\epsilon$ model to describe zero-mean-shear turbulence, because this model was conceived in the framework of turbulence with mean-shear and the well-known constants of the steady $k-\epsilon$ model were determined from experiments with turbulent mean-shear flows. This debate is connected to another very long one about the decay laws for turbulent kinetic energy in oscillating grid experiments (Hopfinger and Toly, 1976; Nokes, 1988; DeSilva and Fernando, 1994; Huppert et al., 1995; among others). Our purpose is not to add another stone in these debates. The $k-\epsilon$ model was used before by Sonin (1983) and Matsunaga et al. (1999) to describe purely diffusive turbulence. When expressing the decay law in oscillating grid turbulence in the form of a power law (i.e. $k \propto z^n$), one observes that the exponent $n$ predicted by the $k-\epsilon$ model in clear water (with no particles) is different from the values given by various authors in light of experimental results. However, Matsunaga et al. (1999) have shown that the results of the $k-\epsilon$ model in clear water compares well with oscillating grid turbulence experiments, at least in the range of experimental conditions available, and that the usual constants of the $k-\epsilon$ model are acceptable for modelling of oscillating grid turbulence. This is discussed in more detail in Section 2.2 of this paper.
The purpose of the present paper is to investigate how a particle suspension modifies the decay of turbulence as compared with the decay in clear water, and to understand how this will produce lutocline formation. The $k-\epsilon$ model and the computation procedure are presented in Section 2, with an emphasis on its application to oscillating grid experiments. Sections 3 and 4 are devoted to the results of numerical computations. Although this work was initially motivated by muddy concentrated benthic suspensions, and cases of computations were chosen from the experimental conditions published in Gratiot (2000) and Mory et al. (2002), our results are not specifically interpreted in terms of cohesive sediments. The range of application of the present study is much concerned with fine-grained sediments, including cohesive sediments. For the present modelling, however, cohesive and non-cohesive sediments are considered in the same way, because the sediment properties enter only through the settling velocity of sediments. The application of the results of this study will only be more difficult for cohesive sediments as the settling velocity of sediment may vary due to flocculation effects. Modelling of flocculation effects depending on the turbulence level is beyond the scope of the present study. Section 5 compares the results of our computations with a wide range of grid oscillation particle suspensions experiments, in particular those of E and Hopfinger (1987) and Huppert et al. (1995), which deal with non-cohesive particles.

## 2. A steady $k-\epsilon$ model for particle suspensions

### 2.1. Formulation

We consider the case of non-cohesive particles which are maintained in suspension in a stirring tank by an oscillating grid. The configuration is sketched in Fig. 1. The vertical axis is oriented upwards, $z=0$ is the grid mean position. The grid oscillation generates a turbulent flow. We assume here that there is no mean flow inside the tank. For constant stirring conditions, a steady-state particle suspension is established which may take the forms of two different typical vertical distributions of concentration. In the first case, sketched in Fig. 1a, which is basically obtained when the particle settling velocity is sufficiently large, the particle settling flux is large and cannot be balanced by the turbulent buoyancy flux, unless the mean concentration decreases gradually with increasing distance from the grid. In the second case (Fig. 1b), which is obtained for particles having a small settling velocity, the particle settling flux is much smaller than in the first case. Because the turbulent diffusivity is high, the equilibrium of the settling flux with the turbulent buoyancy flux is obtained in the suspension with the mean concentration in the suspension decreasing very slowly with increasing distance from the grid. This suspension is established up to the lutocline position, denoted $z_m$, which is constant in time. The lutocline is the interface in a two-layer system where the lower layer is the particle suspension and the upper layer is clear water containing no particles. For fine-grained sediments (Fig. 1b), it will appear later that the averaged concentration in the suspension is almost homogeneous in the lower layer except in the vicinity of the lutocline, where the particle concentration suddenly drops. The equilibrium of the suspension layer is characterised by a relationship between $z_m$, $C_0$ (the mean concentration at $z=0$), the particle settling velocity $w_s$ and the grid oscillation parameters (frequency $f$, stroke $S$, mesh $M$).

We assume in our model that all averaged quantities vary only with the position along the vertical axis and that they are homogeneous in a plane perpendicular to the vertical axis. The equilibrium
<figure><img src="image_2.png" /><figcaption>Fig. 1. Sketch of the grid experiment. The oscillating grid maintains a steady sediment suspension: forming a lutocline at the position $z_m$ (b) or not (a). Circular arrows represent turbulence and the density of dots represents the concentration.</figcaption></figure>

in the suspension layer is given by the mass balance equation

$$0 = \frac{d}{dz} [w_s C - \overline{w'c'}], \qquad (1)$$

where $C(z)$ is the mean concentration at elevation $z$, and $w'$ and $c'$ are the fluctuations of vertical velocity and sediment concentration. The averaged upward turbulent solid particle flux is modelled as

$$\overline{w'c'} = -\Gamma_t \frac{dC}{dz}. \qquad (2)$$

The eddy diffusivity $\Gamma_t(z)$ is

$$\Gamma_t = \frac{v_t}{\sigma_c}. \qquad (3)$$

The eddy viscosity $v_t(z)$ is related to the turbulent kinetic energy $k(z)$ and the rate of dissipation of the turbulent kinetic energy $\varepsilon(z)$ as

$$v_t = c_\mu \frac{k^2}{\varepsilon}. \qquad (4)$$

The vertical distributions of $k(z)$ and $\varepsilon(z)$ are determined from the classical $k-\varepsilon$ equations (e.g. Rodi, 1984; Winterwerp, 2002). The steady-state equations without mean flow are written below. The equation of conservation of the turbulent kinetic energy $k$ is

$$0 = \frac{d}{dz} \left[ \left( v + \frac{v_t}{\sigma_k} \right) \frac{dk}{dz} \right] - \varepsilon + \frac{\rho_s - \rho_w}{\rho_s \rho_w} g \Gamma_t \frac{dC}{dz}. \qquad (5)$$

The equation of conservation of the dissipation rate $\varepsilon$ of turbulent kinetic energy has a similar form

$$0 = \frac{d}{dz} \left[ \left( v + \frac{v_t}{\sigma_\varepsilon} \right) \frac{d\varepsilon}{dz} \right] - c_2 \frac{\varepsilon^2}{k} + (1 - c_3) \frac{\varepsilon}{k} \frac{\rho_s - \rho_w}{\rho_s \rho_w} g \Gamma_t \frac{dC}{dz}, \qquad (6)$$

where $g$ is the acceleration of gravity, $\rho_s = 2.65$ kg/l and $\rho_w = 1$ kg/l are taken as the particle and water densities, respectively. We use the usual values for the $k-\varepsilon$ model constants as

$$\sigma_k = 1, \quad \sigma_\varepsilon = 1.3, \quad c_\mu = 0.09 \quad \text{and} \quad c_2 = 1.92. \qquad (7)$$
The values of the constants appearing in (3) and (6) to account for buoyancy effects are not as usual. We take

$$ \sigma_c = 0.7 \quad \text{and} \quad c_3 = 0. \qquad (8) $$

Winterwerp (2002) suggested that $c_3 = 1$ for a stable stratification. We instead show in Appendix A that the buoyancy term in (6) cannot be zero for a solution to exist in the vicinity of the lutocline in the framework of the standard model for fully developed turbulence. We therefore simply take $c_3 = 0$. Note that the numerical results are not very sensitive to the values of $\sigma_c$ and $c_3$ and their choice do not preclude the generality of the conclusions of the present paper. This was shown by additional computations not presented in this paper for clarity.

Eqs. (1)–(6) can be reduced to a set of three equations, which is more tractable. The three quantities $C(z)$, $k(z)$ and $\varepsilon(z)$ are the solutions of

$$ 0 = w_s C + \frac{c_\mu}{\sigma_c} \frac{k^2}{\varepsilon} \frac{dC}{dz}, \qquad (9) $$

$$ 0 = D - \varepsilon - B = \frac{d}{dz} \left[ \left( v + \frac{c_\mu}{\sigma_k} \frac{k^2}{\varepsilon} \right) \frac{dk}{dz} \right] - \varepsilon - G w_s C, \qquad (10) $$

$$ 0 = \frac{d}{dz} \left[ \left( v + \frac{c_\mu}{\sigma_\varepsilon} \frac{k^2}{\varepsilon} \right) \frac{d\varepsilon}{dz} \right] - c_2 \frac{\varepsilon^2}{k} - (1 - c_3) \frac{\varepsilon}{k} G w_s C \qquad (11) $$

with $G = g(\rho_s - \rho_w)/\rho_s\rho_w$. Eq. (9) expresses the momentum balance: the upward turbulent buoyancy flux balances the downward settling flux of particles. The spatial distribution of the turbulent kinetic energy is the result of three different effects. We identify in (10) the diffusion of the turbulent kinetic energy ($D$), the dissipation ($\varepsilon$) and the buoyancy effect ($B$). A similar analysis holds for the distribution of the turbulent kinetic energy dissipation rate (Eq. (11)).

Introducing the local Rouse number

$$ Rou(z) = \frac{w_s}{\sqrt{k(z)}} \qquad (12) $$

and the local Richardson number

$$ Ri(z) = \frac{G \sqrt{k(z)} C(z)}{\varepsilon(z)}, \qquad (13) $$

we observe from (10) that the buoyancy effect is negligible for the decay of the turbulent kinetic energy when $Ri(z)Rou(z) \ll 1$. The decay of the turbulent kinetic energy is almost the same as in clear water. A dimensional analysis of (10) shows that the characteristic lengthscale of the decay of the turbulent kinetic energy is the integral lengthscale

$$ Z_k \approx \frac{k^{3/2}}{\varepsilon} = l. \qquad (14) $$

In a similar way, a dimensional analysis of (9) shows that the mean concentration varies vertically over a distance of the order of

$$ Z_C \approx \frac{k^{1/2}}{w_s} l = \frac{l}{Rou}, \qquad (15) $$

which is much larger than $Z_k$ when $Rou(z) \ll 1$. In this case, we will observe that the concentration is nearly homogeneous in the lower part of the suspension layer where $Rou(z)Ri(z) \ll 1$, and that
the decay of turbulence is almost the same as in clear water. The quantity $Ri(z)Rou(z)$ increases with increasing distance from the grid. When it becomes of the order of 1, the decay of turbulence becomes much more rapid. A sudden decrease in the concentration occurs while the Rouse number increases rapidly. The lutocline settles down.

## 2.2. The solution in the lower part of the suspension layer in terms of grid oscillation parameters

We determine here the distribution of sediment concentration $C(z)$, turbulent kinetic energy $k(z)$ and dissipation rate $\varepsilon(z)$ in the lower part of the suspension layer $0 < z < z_1$. The position $z_1$ is currently undetermined but we assume that $Ri(z)Rou(z) \ll 1$ for $0 < z < z_1$.

We refer in the present section to a paper by Matsunaga et al. (1999). Using the same $k-\varepsilon$ model as we do, they determined the distributions of $k(z)$ and $\varepsilon(z)$ produced in a grid oscillation experiment in clear water. On the one hand, Matsunaga et al.'s paper contains a useful discussion on the validity of the $k-\varepsilon$ model for modelling of purely diffusive turbulence. On the other, their solution establishes the dependence of the distribution of $k(z)$ and $\varepsilon(z)$ with the grid oscillation parameters (grid mesh $M$, frequency $f$ and stroke $S$ of grid oscillation). We shall observe that the solution in the lower part of the suspension layer is almost the same as in clear water.

The validity of the $k-\varepsilon$ model to describe the transportation of turbulence by the sole effect of diffusion is questionable. This model was conceived by considering flows in which the distribution of turbulence inside the domain is governed by the production of turbulence from the mean shear flow and its constants were determined from experimental results obtained for such flows. The case of turbulence produced in an oscillating grid experiment is rather different. When the geometrical properties of the grid are suitably chosen, the mean flow is very weak in the tank, except in a very limited region close to the grid. Turbulence is transported from the grid by the single effect of diffusion. The ability of a $k-\varepsilon$ model to describe the decay of turbulence with increasing distance from the grid was considered before by Sonin (1983). We observe that the solution obtained by Matsunaga et al. is basically the same as Sonin's. Both authors showed that the $k-\varepsilon$ model predicts decay laws for the turbulent kinetic energy and its dissipation rate in the form of power laws. The power-law exponents they obtained are the same. In addition, Matsunaga et al. compared their solutions with a set of various laboratory experiments, and related their numerical solutions to the parameters of grid oscillation. The exponents of the decay laws predicted by the $k-\varepsilon$ model for the turbulent kinetic energy and the dissipation rate are -4.98 and -8.46, respectively, whereas the exponents suggested by several authors (Hopfinger and Toly, 1976; Ura et al., 1987; DeSilva and Fernando, 1994) from their experiments were -2 and -4, respectively. This is a significant difference from the mathematical point of view, but Matsunaga et al. showed that, on the basis of available experimental data, the results of the $k-\varepsilon$ model with its widely accepted model constants are also applicable to oscillating grid turbulence.

It is not the purpose of the present paper to address again the decay of turbulence produced by an oscillating grid. Our aim is to investigate the conditions of existence of a lutocline separating a suspension layer and a clear water layer using the $k-\varepsilon$ model. Considering that Matsunaga et al. provided clear proof that the $k-\varepsilon$ model is able to address oscillating grid turbulence, at least for the conditions of available experimental data, the $k-\varepsilon$ model with its usual constant is used here. We additionally refer to the quantification of the turbulent kinetic energy $k_0$ and dissipation rate $\varepsilon_0$ at the mean level position determined by Matsunaga et al. in terms of the grid oscillation parameters
$M$, $S$ and $f$. For $Re = fS^2/v \ge 5500$ they obtained

$$k_0 = 0.6 f^2 S^2 (S/M)^{1/4} \quad \text{and} \quad \varepsilon_0 = 0.45 f^3 S^2 S/M. \tag{16}$$

We only consider this case here. Using (16) the local Rouse number

$$Rou_0 = \frac{w_s}{\sqrt{k_0}} \tag{17}$$

and the local Richardson number

$$Ri_0 = \frac{G\sqrt{k_0}C_0}{\varepsilon_0} \tag{18}$$

at the level of the grid mid position can be expressed as bulk parameters of the experiment which compare the properties of the suspension to the grid oscillation parameters.

Following Matsunaga et al. (1999), it is straightforward to show that the decay of turbulence is modified in a negligible way in the lower part of the suspension layer as compared to the solution in clear water when $Rou_0 Ri_0 \ll 1$. Introducing the following dimensionless quantities:

$$\hat{z} = \frac{z}{z_0} \quad \text{with} \quad z_0 = \frac{k_0^{3/2}}{\varepsilon_0}, \quad \hat{k} = \frac{k}{k_0}, \quad \hat{\varepsilon} = \frac{\varepsilon}{\varepsilon_0}, \quad \hat{v}_t = \frac{v_t \varepsilon_0}{k_0^2} \quad \text{and} \quad \hat{C} = C/C_0, \tag{19}$$

and using a new independent variable

$$\frac{d\hat{\zeta}}{d\hat{z}} = \frac{\sigma_k}{\hat{v}_t} \quad \text{for } 0 < \hat{z} < \infty, \quad \text{with } \hat{\zeta}(\hat{z}=0) = \sqrt{\frac{6\sigma_k}{c_\mu}}, \tag{20}$$

Eqs. (9)–(11) of the $k-\varepsilon$ model can be rewritten in a dimensionless form as

$$\frac{d\hat{C}}{d\hat{\zeta}} = -\frac{\sigma_c}{\sigma_k} Rou_0 \hat{C}, \tag{21}$$

$$\frac{d^2\hat{k}}{d\hat{\zeta}^2} = \frac{c_\mu}{\sigma_k} \left[ \hat{k}^2 + Ri_0 Rou_0 \frac{\hat{k}^2}{\hat{\varepsilon}} \hat{C} \right], \tag{22}$$

$$\frac{d^2\hat{\varepsilon}}{d\hat{\zeta}^2} = \frac{c_\mu \sigma_\varepsilon}{\sigma_k^2} \left[ c_2 \hat{\varepsilon} \hat{k} + (1 - c_3) Ri_0 Rou_0 \hat{k} \hat{C} \right]. \tag{23}$$

The kinematic viscosity is neglected in (22) and (23) as compared to the eddy viscosity because this set of equations will only be considered in the lower part of the suspension layer. In clear water, $\hat{C} = 0$ and Eqs. (22) and (23) reduce to a simple set of equations, whose solution was given by Matsunaga et al. (1999) to be

$$k_{cw}(z) = k_0 \left( 1 + \frac{s}{6} \sqrt{\frac{6\sigma_k}{c_\mu}} \frac{z}{z_0} \right)^{-2/s}, \tag{24}$$

$$\varepsilon_{cw}(z) = \varepsilon_0 \left( 1 + \frac{s}{6} \sqrt{\frac{6\sigma_k}{c_\mu}} \frac{z}{z_0} \right)^{-(3+s)/s} \tag{25}$$
with

$$s = -\frac{7}{2} + \sqrt{\frac{1}{4} + \frac{6c_2\sigma_{\varepsilon}}{\sigma_k}} \qquad (26)$$

The numerical procedure used to solve the system of Eqs. (21) and (22), and in particular the values of the boundary conditions, is described in the next section.

## 2.3. Numerical method and boundary conditions

A specific procedure must be used to determine the solution in the suspension layer because the position $z_m$ of the lutocline is not known. The position $z_m$ must be determined as part of the solution. The solution of the system of Eqs. (21)–(23) in the domain $0 < z < z_1$, and of the system of Eqs. (9)–(11) in the domain $z_1 < z < z_m$ has to verify the following boundary conditions:

$$C(z=0) = C_0, \quad k(z=0) = k_0 \quad \text{and} \quad \varepsilon(z=0) = \varepsilon_0, \qquad (27)$$

$$C(z_m) = 0, \quad k(z_m) = 0 \quad \text{and} \quad \varepsilon(z_m) = 0. \qquad (28)$$

In order to find the solution and simultaneously determine $z_m$ systems (9)–(11) and (21)–(23) is solved using a shooting technique to determine the lutocline position. The equations are digitised by a Runge-Kutta finite difference scheme. The integration is performed by starting from the lower boundary toward the lutocline for the boundary conditions (27) and by assuming empirical values for the derivatives $[dk/dz]_0 = dk/dz(z=0)$ and $[d\varepsilon/dz]_0 = d\varepsilon/dz(z=0)$. For any choice of boundary conditions, the computation is stopped at a position $z$ if one of the values $C(z)$, $k(z)$ or $\varepsilon(z)$ becomes negative or zero, whereas the other two remain positive. Using a ‘trial and error’ technique, new values are given for the derivatives of $k$ and $\varepsilon$ at $z=0$, and a subsequent integration of the equations is performed. The procedure is repeated until the values of the derivatives $[dk/dz]_0$ and $[d\varepsilon/dz]_0$ have been found, for which integration provides a solution where $C(z)$, $k(z)$ or $\varepsilon(z)$ vanish simultaneously within a grid mesh at $z=z_m$. The solution obtained is valid if $k$, $\varepsilon$ and $C$ decrease monotonically in the interval $0 < z < z_m$. A subsequent check of the validity of the solution is made by considering the flux of the turbulent kinetic energy at level $z$,

$$F(z) = - \left( v + \frac{v_t}{\sigma_k} \right) \frac{dk}{dz}. \qquad (29)$$

Integration of the $k$ equation (5) over the vertical gives

$$F(z) = F(0) - \int_{0}^{z} \varepsilon \, dz - G w_s \int_{0}^{z} C \, dz. \qquad (30)$$

The consistent solution is obtained if the flux $F(z)$ tends numerically to zero at $z = z_m$.

Close to the grid the decay of turbulence is modified by the particle suspension only to a very limited extent as shown by Eqs. (22) and (23). For the cases computed in this study the bulk Rouse number $Rou_0$ and the bulk Richardson number $Ri_0$ ranged between $10^{-4}$ and $10^{-1}$, and between $10^{-2}$ and $10^2$, respectively. The quantity $Rou_0 Ri_0$ varied from $10^{-4}$ to $10^{-1}$. The derivatives $[dk/dz]_0$ and $[d\varepsilon/dz]_0$ at the grid mid position $z=0$ were almost the same as in clear water. Writing

$$\left[\frac{dk}{dz}\right]_0 = \left[\frac{dk_{cw}}{dz}\right]_0 (1 + \delta_k) \quad \text{and} \quad \left[\frac{d\varepsilon}{dz}\right]_0 = \left[\frac{d\varepsilon_{cw}}{dz}\right]_0 (1 + \delta_\varepsilon), \qquad (31)$$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Computed values of $\delta_k$ (○) and $\delta_\epsilon$ (□) for the numerical solutions.</figcaption></figure>

the results of the different cases computed showed that $\delta_k \approx \text{Rou}_0 \text{Ri}_0$ and $\delta_\epsilon \approx \text{Rou}_0 \text{Ri}_0/c_2$. Fig. 2 shows the dependence of $\delta_k$ and $\delta_\epsilon$ with $\text{Rou}_0 \text{Ri}_0$ obtained from our computations when the position $z_m$ of the lutocline was found, where the concentration, turbulent kinetic energy and dissipation rate all vanish. Slightly different values of $\delta_k$ and $\delta_\epsilon$ lead to non-valid solutions (for instance, $C, k$ or $\epsilon$ are not decreasing in all the suspension with increasing distance from the grid, or one of the latter quantities becomes negative). The values of $\delta_k$ and $\delta_\epsilon$ are adjusted just to ensure the upper boundary condition (28) at $z = z_m$. The quantities $\delta_k$ and $\delta_\epsilon$ must be very small because, for the conditions investigated, the turbulence in the vicinity of the grid is hardly modified by the presence of particles as compared to the decay of turbulence in clear water. Because of this, we have no simple understanding of the dependence of $\delta_k$ and $\delta_\epsilon$ with $\text{Rou}_0 \text{Ri}_0$ displayed in Fig. 2.

### 3. The distributions of particle concentration and turbulence in the suspension layer

Fig. 3 shows the computed vertical variations in concentration, turbulent kinetic energy and dissipation rate in the suspension layer for seven particular cases, which are representative of different behaviours. For all cases plotted, the grid oscillation parameters are the same. Different sediment concentrations $C_0$ and different settling velocities $w_s$ distinguish the three conditions. Three different bulk Rouse numbers and three different bulk Richardson numbers are considered, namely $\text{Rou}_0 = 10^{-3}$, $10^{-2}$, $10^{-1}$ and $\text{Ri}_0 = 10^{-1}$, 1, 10. The different combinations of values of $\text{Rou}_0$ and $\text{Ri}_0$ makes the quantity $\text{Ri}_0 \text{Rou}_0$ vary from $10^{-3}$ to $10^{-1}$. The vertical variations in the local Richardson number $\text{Ri}(z)$, the local Rouse number $\text{Rou}(z)$ and the local values of $\text{Ri}(z) \text{Rou}(z)$ are also shown in the lower part of Fig. 3 in order to help the interpretation. The different cases plotted in Fig. 3 clearly show the effect of the bulk Rouse number on the particles and turbulence distributions inside the suspension. The vertical concentration profiles change considerably with the Rouse number $\text{Rou}_0$. For the lowest values of the Rouse number ($\text{Rou}_0 = 10^{-3}$ and $10^{-2}$) the main observation is that a lutocline is obtained and the concentration is almost uniform in the suspension layer. For the largest Rouse number ($\text{Rou}_0 = 10^{-1}$) the concentration decreases regularly inside the suspension layer with increasing distance from the grid and no clear lutocline is observed in the concentration profile.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Variations in (a) the dimensionless concentration $C(z)$, (b) turbulent kinetic energy $k(z)$ and (c) dissipation rate $\varepsilon(z)$ with distance $z$ from the grid mean position for different properties of particles. The vertical variations in (d) the local Richardson number $Ri(z)$ (Eq. (13)), (e) local Rouse number $Rou(z)$ (Eq. (12)) and (f) local values of $Ri(z)$ $Rou(z)$ are shown in addition. The same grid oscillation conditions are used for all curves plotted: $f=3$ Hz, $S=4.5$ cm, $M=7.5$ cm, so that $k_0 = 96 \text{ cm}^2 \text{ s}^{-2}$, $\varepsilon_0 = 148 \text{ cm}^2 \text{ s}^{-3}$ and $z_0 = 6.4$ cm. Seven different cases with particles are plotted, where the type of lines indicates the different values of $Rou_0$ while the thickness indicates the different values of $Ri_0$: $Rou_0 = 10^{-3}$ (dashed-dotted lines), $10^{-2}$ (solid lines), $10^{-1}$ (dashed lines), and $Ri_0 = 10^{-1}$ (doubly thick lines), 1 (thin lines), 10 (thick lines). The vertical variations in $k_{cw}$ and $\varepsilon_{cw}$ in clear water are also shown in b and c (dotted line).</figcaption></figure>

The bulk Rouse number $Rou_0$ therefore appears to be the parameter that determines the shape of the concentration profile and the conditions for which a homogeneous concentrated suspension is obtained.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Vertical profile of the eddy viscosity divided by the kinematic viscosity for the conditions and symbols of Fig. 3.</figcaption></figure>

The vertical variations in $k_{cw}$ and $\varepsilon_{cw}$ in clear water are also shown in Fig. 3 (dotted line). Because $Ri_0 Rou_0 \ll 1$ for all cases shown in Fig. 3, we observe that the decay of turbulence is not much affected by the presence of particles in the lower part of the suspension layer ($\hat{z} = z/z_0 < \hat{z}_1 = 0.5$), as discussed in Section 2.2. The variations in $k$ and $k_{cw}$ and the variations in $\varepsilon$ and $\varepsilon_{cw}$ are almost superimposed in the suspension layer, respectively, except in the vicinity of the lutocline.

The vertical profiles of the eddy viscosity $v_t$ are plotted in Fig. 4 for the seven cases shown in Fig. 3. This quantity decreases rapidly when the lutocline is approached. Winterwerp (1999) discussed that low Reynolds damping functions are needed when $v_t < 50v$. For the cases considered here, $v_t > 50v$ as far as position $z$ located less than 1 cm below the lutocline. At this level, the rapid decrease in $k$ and $\varepsilon$ associated with the occurrence of the lutocline is already settled. Using low Reynolds damping functions (e.g. Patel et al., 1985) would presumably not significantly modify the position of the lutocline. Additionally, a $k-\varepsilon$ model assumes that turbulence is isotropic, and this might also be a shortcoming for modelling the region very close to the lutocline. Anyhow, we do not pretend to describe the interface in all its complexity but rather the conditions leading to its formation. Because the $k-\varepsilon$ model shows that concentration and turbulence decrease very suddenly in the vicinity of the lutocline, a more relevant and sophisticated model for describing low Reynolds number and anisotropic turbulence should be able to sustain significantly higher level of turbulence to predict a position of the lutocline significantly different from that given by the $k-\varepsilon$ model.

For the same conditions as those of Figs. 3 and 4, the vertical profile of the turbulent kinetic energy flux $F(z)$ (Eq. (29)) is finally displayed in Fig. 5. We verify that this quantity vanishes at
<figure><img src="image_6.png" /><figcaption>Fig. 5. Vertical profile of the dimensionless flux of turbulent kinetic energy for the conditions and symbols of Fig. 3.</figcaption></figure>

the lutocline when such an interface is observed. As discussed in Section 2.3, computing $F(z)$ is a check of numerical convergence. We considered that numerical convergence was obtained when the energy flux at level $z_m$ was less than $10^{-5}$ of the flux $F_0$ at the grid mean position.

For the cases considered in Figs. 3–5, the lutocline is a very sudden transition between the suspension layer and the clear water layer, and the concentration is almost uniform inside the suspension layer when the Rouse number $Rou_0$ is sufficiently small. The local flux Richardson number, defined as the ratio of the buoyancy effects to the available turbulent energy in (10),

$$Rf = \frac{B}{D} = \frac{Gw_s C(z)}{\frac{d}{dz}[(v + (c_\mu/\sigma_k)(k^2/\epsilon))\frac{dk}{dz}]}, \quad (32)$$

provides an alternative estimate of the effect of buoyancy on the decay of turbulence. Because $D = B + \epsilon$, the flux Richardson number also takes the simple form as

$$Rf = \frac{B}{D} = \frac{Gw_s C(z)}{\epsilon(z) + Gw_s C(z)}, \quad (33)$$

or is alternatively expressed in terms of the local Rouse number and of the local Richardson number as

$$Rf = \frac{B}{D} = \frac{Rou(z)Ri(z)}{1 + Rou(z)Ri(z)}. \quad (34)$$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Vertical variations in the flux Richardson number for the conditions and symbols of Fig. 3.</figcaption></figure>

For the seven conditions already considered, Fig. 6 shows the vertical variations in the flux Richardson number. Depending on the conditions, different behaviours can be noticed. We clearly observe that the flux Richardson number peaks to one at a position very close to the lutocline for all cases where one is observed, implying that the turbulent kinetic energy dissipation rate becomes very low and the particle suspension suddenly drops to zero. On the other hand, the profiles of the flux Richardson number for computations with the bulk Rouse number $R_{u0} = 10^{-1}$, which do not show the occurrence of a lutocline, display a smoother behaviour. As discussed in Section 1, the observation that the flux Richardson number may approach the maximum value of 1 in grid turbulence experiments is a striking property, which differs strongly from observations in turbulence with mean shear flows. Larger values of the flux Richardson number imply that the dissipation rate is smaller as compared to the buoyancy term. As the energy is transmitted by diffusive turbulence at larger scales far from the grid, we expect that the dissipation rate is reduced for a same level of the turbulent energy, whereas mixing is more efficient.

## 4. The depth of particle suspension layers

Determining the depth of the particle suspension layer as a function of the bulk Rouse number $R_{u0}$ and the bulk Richardson number $Ri_0$ is of practical significance by providing a simple and direct quantification of experimental observations. Fig. 3 has already highlighted different behaviours
<figure><img src="image_8.png" /><figcaption>Fig. 7. Dimensionless position of the lutocline $z_m/z_0$ as a function of $Rou_0 Ri_0$. △: computations for $Rou_0 > 0.01$; ▽: computations for $Rou_0 < 0.01$.</figcaption></figure>

depending on the value of $Rou_0$. On the one hand, when $Rou_0$ is small, the particle concentration cannot decrease below a minimum distance from the grid, where turbulence level has sufficiently decreased due to dissipation and buoyancy effects. The buoyancy effect suddenly enhances the decay of turbulence and a lutocline is formed (as sketched in Fig. 1b). On the other hand, when $Rou_0$ is sufficiently large, the concentration inside the suspension decreases simultaneously with the decay of turbulence. The buoyancy effect is limited by the decrease in concentration, and no lutocline is formed (Fig. 1a). Some consequences are visible on Fig. 3 and are interpreted in terms of the depth of the suspension layer. When $Rou_0$ is sufficiently small ($Rou_0 = 10^{-3}$ and $10^{-2}$ in Fig. 3) and a lutocline is observed, the depth of the suspension layer decreases with increasing values of $Rou_0 Ri_0$ as the decay of turbulence is enhanced by the buoyancy effect. However, when $Rou_0$ is sufficiently large ($Rou_0 = 10^{-1}$ in Fig. 3), the particle concentration decreases regularly with increasing distance from the grid. The turbulence can maintain particles in suspension at a higher distance from the grid as compared to the depth of the suspension layer obtained for the same value $Rou_0 Ri_0$ but for a smaller $Rou_0$.

The different behaviours obtained from our numerical computations are summarised in Fig. 7, which presents the depth of the suspension layer as a function of $Rou_0 Ri_0$. Two different data sets are identified. The focus in this paper is on computation conditions with $Rou_0 \le 0.01$. A lutocline is clearly identified and the results of the computations shown in Fig. 7 indicate a regular decrease in $z_m/z_0$ with increasing values of $Rou_0 Ri_0$. The computational results are compared with the experimental results in Section 5. For $Rou_0 > 0.01$, on the other hand, higher values are obtained for the quantity $z_m/z_0$. The lutocline position is not well defined in this case, because there is no sharp transition of density. It was determined, as explained in Section 2.3, as the first position where the concentration becomes zero or negative. This mathematical definition lacks a physical basis and the values of $z_m/z_0$ shown in Fig. 7 for $Rou_0 > 0.01$ are somewhat indicative. Only a few computations were performed for $Rou_0 > 0.01$; their purpose was to identify the difference with the case $Rou_0 < 0.01$.
5. Comparison with experiments

Our numerical study was motivated by observations made in recent oscillating grid experiments
with mud suspensions (Gratiot, 2000; Mory et al., 2002). A variety of suspension layers were ob-
served, for the same grid oscillation conditions, depending on the concentration level inside the sus-
pension. Significant variations in the depth of the particle suspension layer were measured. Cohesive
sediment suspensions are particularly difficult to characterize in the laboratory because measurement
of the particle settling velocity is tedious, as it varies with the concentration. The determination of an
averaged settling velocity certainly does not fully account for the deposition of flocs having differ-
ent size and buoyancy properties. Because of the uncertainty of settling velocity measurements, the
results of our k-ε model are not compared directly with the experimental results presented by Mory
et al. (2002), as this quantitative comparison does not allow clear conclusions to be drawn. However,
the range of conditions for which numerical computations were performed has been defined in view
of the latter experiments, considering large variations in the concentration and settling velocity. A
striking result of mud suspension experiments carried out by Mory et al. (2002) was the observation
that in the vicinity of the lutocline the flux Richardson number reaches a value close to one. This
was an unexpected result as it differs significantly from previous experiments on suspensions in
mean-shear turbulent flows. Ivey and Imberger (1991) pointed out that the flux Richardson number
is always below 0.25. Actually, our numerical computations give some support to the argument that
grid turbulence can generate particle suspensions where the flux Richardson number is close to one.

Experiments with non-cohesive particle suspensions in an oscillating grid tank have been carried
out by E and Hopfinger (1987) and Huppert et al. (1995). It is more straightforward to compare
their results with our computations as the settling velocity of the particles is more easily determined
than in the case of cohesive sediments. Moreover, E and Hopfinger (1987); Gratiot (2000) and Mory
et al. (2002) used the same grid tank set-up. The position of the lutocline obtained in the experiments
of E and Hopfinger, and in those of Huppert et al. is compared in Fig. 8 with our numerical results
presented in Fig. 7. In Figs. 8a and b, the values of the dimensionless ratio z_m/z_0 are plotted as
a function of Ri_0 and of Rou_0, respectively. The two figures do not provide any insight into the
dependence of z_m/z_0 on Ri_0 and Rou_0, but they indicate the range of conditions considered by the
different authors. The experiments by E and Hopfinger had a high Rou_0, and Rou_0 Ri_0 was also
rather large. For such conditions our model predicts that the decay of turbulence is rapidly enhanced
by the buoyancy, and there is no well-defined lutocline. The concentration decreases regularly with
increasing distance from the grid (see the dashed lines in Fig. 3). Fig. 8c shows the variations in
z_m/z_0 with Rou_0 Ri_0. E and Hopfinger's data are clearly above the other data. Figs. 8a and b indicate
that Huppert et al. investigated a wider range of conditions that coincide with our computation cases,
except when Rou_0 Ri_0 is of the order of one. The data of Huppert et al. show a similar dependence
of z_m/z_0 on Rou_0 Ri_0 to that obtained in our computations. Examining the data in more detail, we
find that Huppert et al.'s experimental data for Ri_0 Rou_0 < 2 × 10⁻³ overlap satisfactorily the results
of our computations with Rou_0 > 0.01, whereas the data for Ri_0 Rou_0 > 2 × 10⁻³ presents a better
agreement with our computations with Rou_0 < 0.01, although Fig. 8b shows that Huppert et al.'s
results with Ri_0 Rou_0 < 2 × 10⁻³ correspond to experiments with Rou_0 < 0.01.

Fig. 8 shows that both our computations with *Rou*₀ *<* 0.01 and the experimental results predict
a variation in the depth of the suspension layer with *Rou*₀*Ri*₀ of the form *z*<sub>m</sub>/*z*₀ ∝ (*Rou*₀*Ri*₀)<sup>-1/4</sup>.
This result seems to indicate a general behaviour. Huppert et al. proposed a simple analytical model
<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison of the dimensionless position of the lutocline obtained from $k-\epsilon$ computations with the experimental data by E and Hopfinger (1987) (○) and by Huppert et al. (1995). Huppert et al.'s results were obtained for different particle diameters (+: $\lambda$ = 0.055 mm; ×: $\lambda$ = 0.058 mm; ∴: $\lambda$ = 0.108 mm; *: $\lambda$ = 0.176 mm). △: computations for $Rou_0$ > 0.01; ▽: computations for $Rou_0$ < 0.01.</figcaption></figure>

which establishes that $z_m/z_0 \propto (Rou_0 Ri_0)^{-n}$, where the exponent $n$ is $\frac{1}{4}$ or 0.21, depending on the turbulence decay law adopted. Huppert et al.'s model solves the turbulent kinetic energy equation, assuming that the integral lengthscale varies linearly with distance from the grid (as in clear water) and that the concentration is homogeneous inside the suspension layer. This model does not consider the effect of the Rouse number, but it is partly equivalent to our numerical model when $Rou_0 < 0.01$. Huppert et al. made additional assumptions, which actually do not appear to be required. A modified analytical model, inspired by Huppert et al.'s model, is presented in Appendix B, which shows that the dependency $z_m/z_0 \propto (Rou_0 Ri_0)^{-1/4}$ is quite general. This analytical model differs basically from the numerical model by the fact that the analytical model assumes that the integral lengthscale is not modified in the suspension as compared to clear water. The comparison demonstrates that this crude assumption has no effect on the final result. The more significant limitation of Huppert et al.'s model is that it can only address the case of suspensions with very low Rouse numbers.

# 6. Conclusion

In the comparison of our model with studies carried out by other authors, we have left aside the work by Noh and Fernando (1991), who derived two equations for modelling sediment suspen-
sions at low concentrations in diffusive turbulence. They computed unsteady solutions in the range $0.003 < Rou_0 < 0.1$ and $Ri_0 < 0.01$. This low $Ri_0$ regime does not correspond to standard experimental conditions such as those considered in the present paper. Noh and Fernando's model consists of one equation for $C$ similar to (9) and one equation for $k$. The dissipation rate is defined using the integral lengthscale $l$, which is assumed to vary linearly with increasing distance from the grid. This model is very similar to the model of Huppert et al. (1995) and its modified version presented in Appendix B, but in Noh and Fernando's model the integral lengthscale is limited by the buoyancy lengthscale $l_b = \sqrt{k/N}$ (where $N$ is the Brunt-Väisälä frequency), which decreases rapidly to zero at a certain level that defines the interface.

Our numerical computations highlight the respective effects of the bulk Rouse number $Rou_0$ and the bulk Richardson number $Ri_0$ through the quantity $Rou_0Ri_0$. On the one hand, concentrated suspensions with a homogeneous mean concentration are observed when $Rou_0$ is below 0.01. When $Rou_0Ri_0$ is small, the decay of turbulence is affected by the particle suspension only within a very short distance from the lutocline. The transition leading to the lutocline formation is very sudden and the flux Richardson number is close to 1 in the vicinity of the interface. This confirms the observations by Gratiot (2000) and Mory et al. (2002). On the other hand, a simultaneous decay of sediment concentration and turbulence is observed when $Rou_0$ is greater than 0.01.

The present paper focuses on concentrated sediment suspensions. For $Rou_0 < 0.01$, the depth of the suspension layer decreases with increasing values of $Rou_0Ri_0$. A dependency in the form $z_m/z_0 \propto (Rou_0Ri_0)^{-1/4}$ was obtained from our computation. An interesting observation is that this dependence is the same as that predicted by the analytical model proposed by Huppert et al. (1995), and also obtained using its modified version (Appendix B), when the decay laws suggested by Hopfinger and Toly (1976) are used. The dependence of the analytical model is in the form $z_m/z_0 \propto (Rou_0Ri_0)^{-0.21}$ for the decay laws proposed by Huppert et al. (1995). We observe that the difference is small and in agreement with the conclusion of Matsunaga et al. (1999), which points out that the different models and decay laws for oscillating grid turbulence are not significantly different for the range of parameters usually considered.

## Acknowledgements

This work has been partially funded by the European Commission, Directorate General XII for Science, Research and Development under Contract No. MAS3-CT97-0082 (COSINUS Project). Many thanks to our COSINUS partners, especially E. Toorman, O. Petersen and N. Gratiot for constructive criticisms.

## Appendix A. On the buoyancy term in the $\epsilon$-equation

We consider the suspension case, i.e. the system of Eqs. (9)-(11). If a lutocline exists at $z=z_m$, let us assume that the solution very close to the lutocline (for $Z=z_m-z$ very small) has the following form:

$$k = aZ^p, \quad \epsilon = bZ^q \quad \text{and} \quad C = cZ^r. \tag{A.1}$$
The exponents $p$, $q$ and $r$ are real positive values because the three quantities $k$, $\epsilon$ and $C$ should vanish at the lutocline. Introducing (A.1) into (9), it may be deduced that

$$q = 2p - 1. \tag{A.2}$$

Introducing (A.1) into (10) with (A.2), then gives

$$vap(p-1)Z^{p-2} + \frac{c_{\mu}}{\sigma_k} \frac{a^3}{b} p^2 Z^{p-1} - bZ^{2p-1} - Gw_s c Z^r = 0. \tag{A.3}$$

We note that, for very small $Z$, the first term dominates the second one ($p-2 < p-1$), implying that the kinematic viscosity is stronger than the eddy viscosity very close to the lutocline. Close to the interface, the eddy viscosity model (4) should be corrected due to low Reynolds number effects. However, when using the standard model for fully developed turbulence the viscous diffusive term (first term) dominates the turbulent dissipation (third term) when $p > 1$. We thus necessarily require the buoyancy term to balance the viscous diffusive term and

$$r = p - 2. \tag{A.4}$$

One may check that introducing (A.1) into (11) gives the same result. This leads to an interesting remark: if $c_3$ is set equal to one in (11) (i.e. considering no buoyancy term in the $\epsilon$ equation), there is no physical solution close to the lutocline ($k$ and $\epsilon$ diverge). The exponents $p$, $q$ and $r$ are not determined at this stage but the $k-\epsilon$ equations indicate that buoyancy effects dominate just below the interface.

## Appendix B. Buoyancy flux model applied to suspensions of particles

Huppert et al. (1995) proposed a model predicting the depth of a suspension layer of particles in a grid turbulence experiment as a function of the particle properties and grid oscillation parameters. A modified version of this model is described below, which is more general as it does not require some of the assumptions made by Huppert et al.. This modified model and Huppert et al.'s model share the following properties:

(i) The particle concentration inside the suspension layer is assumed to be homogeneous.

(ii) The integral lengthscale in the suspension layer varies linearly with increasing distance from the grid, i.e. $l = qz$.

(iii) The turbulent kinetic energy equation is written as Eq. (12) in Huppert et al.'s paper,

$$r \frac{d(u^3)}{dz} = -\frac{u^3}{qz} - \frac{g}{\rho_0} \overline{w' \rho'}. \tag{B.1}$$

The notations used in this appendix are those of Huppert et al. in order to highlight the similarities and differences between the two models. $u$ denotes the turbulent rms velocity ($k=u^2$). The buoyancy term and the dissipation term in the turbulent kinetic energy equation (B.1) and in the turbulent kinetic energy equation of the $k-\epsilon$ model (Eq. (5)) are the same ($\epsilon=u^3/qz$). Only the diffusion terms in Eqs. (B.1) and (5) slightly differ.

The main difference between our analytical model and Huppert et al.'s model is that we do not assume that the buoyancy flux decreases linearly from the grid to the edge of the layer as stated by
Huppert et al. (Eq. (13)). Because of the mass balance equation (Eq. (1) in the present paper) the buoyancy flux has to be homogeneous when the particle concentration is homogeneous. Using the mass balance equation, (B.1) is therefore rewritten as

$$r \frac{d(u^3)}{dz} = -\frac{u^3}{qz} - B_0 \quad \text{with} \quad B_0 = \frac{g}{\rho_0} \frac{\rho_p - \rho_0}{\rho_p} w_s C. \qquad (B.2)$$

Searching for a solution of the form $u^3 = \beta z + F(z)$, we simply find that the solution of (B.2) is

$$u^3(z) = B_0 \frac{q}{rq-1} z + \left( u_0^3 - B_0 \frac{q}{rq-1} z_0 \right) \left( \frac{z_0}{z} \right)^{1/rq} \qquad (B.3)$$

when the turbulent velocity at $z = z_0$ is taken to be $u(z_0) = u_0$. This integration of (B.2) does not require any assumption of the form $z/z_0 \gg 1$ as used by Huppert et al. The position of the interface is given by the location where the turbulent velocity vanishes, and therefore

$$\frac{D}{z_0} = \left( 1 + \frac{1-rq}{q} \frac{u_0^3}{B_0 z_0} \right)^{\frac{rq}{1+rq}}. \qquad (B.4)$$

This law takes different forms depending on the value of the quantity $rq$. For $rq = \frac{1}{3}$, as taken by Hopfinger and Toly (1976), we deduce

$$\frac{D}{z_0} = \left( 1 + 2r \frac{u_0^3}{B_0 z_0} \right)^{\frac{1}{4}}, \qquad (B.5)$$

whereas the value $rq = 1/3.84$, suggested by Huppert et al. leads to

$$\frac{D}{z_0} = \left( 1 + 2.84r \frac{u_0^3}{B_0 z_0} \right)^{\frac{0.21}{2}}. \qquad (B.6)$$

We find the same power laws as those found by Huppert et al. The only quantitative difference is the proportionality coefficient. For conditions producing a suspension layer in which the particle concentration is homogeneous, it is reasonable that $u_0^3/B_0z_0 \gg 1$, so that a simple power law dependence of $D/z_0$ versus the grid parameters and particle properties is obtained.

## References

* Briggs, D.A., Ferziger, J.H., Koseff, J.R., Monismith, S.G., 1998. Turbulent mixing in a shear-free stably stratified two-layer fluid. J. Fluid Mech. 354, 175–208.

* DeSilva, I.P.D., Fernando, H.J.S., 1994. Oscillating grids as a source of nearly isotropic turbulence. Phys. Fluids 6, 2455–2464.

* E, X., Hopfinger, E.J., 1987. Stratification by solid particle suspensions. Third International Symposium on Stratified Flows, Caltech, Pasadena, pp. 1–8.

* Gratiot, N., 2000. Etude expérimentale de la formation des couches de crème de vase turbulentes. Ph.D. Thesis, Université Joseph Fourier, Grenoble.

* Hopfinger, E.J., Linden, P.F., 1982. Thermoclines in zero-mean-shear turbulence. J. Fluid Mech. 114, 157–173.

* Hopfinger, E.J., Toly, J.A., 1976. Spatially decaying turbulence and its relation to mixing across density interfaces. J. Fluid Mech. 78, 155–175.

* Huppert, H.E., Turner, J.S., Hallworth, M.A., 1995. Sedimentation and entrainment in dense layers of suspended particles stirred by an oscillating grid. J. Fluid Mech. 289, 263–293.
Ivey, G.N., Imberger, J., 1991. On the nature of turbulence in a stratified fluid. Part 1: the energetics of mixing. J. Phys. Oceanogr. 21, 650-658.

Matsunaga, N., Sugihara, Y., Komatsu, T., Masuda, A., 1999. Quantitative properties of oscillating-grid turbulence in a homogeneous fluid. Fluid Dyn. Res. 25, 147-165.

Mory, M., Gratiot, N., Manning, A., Michallet, H., 2002. CBS layers in a diffusive turbulence grid oscillation experiments. In: Winterwerp, J.C., Kranenburg, C. (Eds.), Proceedings in Marine Science, Vol. 5. Fine Sediment Dynamics in the Marine Environment, Elsevier Science B.V., Amsterdam, pp. 139-154.

Noh, Y., Fernando, H.J.S., 1991. Dispersion of suspended particles in turbulent flow. Phys. Fluids A 3 (7), 1730-1740.

Nokes, R.I., 1988. On the entrainment rate across a density interface. J. Fluid Mech. 188, 185-204.

Patel, V.C., Rodi, W., Scheuerer, G., 1985. Turbulence models for near-wall and low Reynolds number flows: a review. AAIA J. 23 (9), 1308-1319.

Rodi, W., 1984. Turbulence models and their applications in hydraulics—a state-of-the-art review. IAHR monograph, Delft.

Sonin, A.A., 1983. A calibration of the model for the diffusion of turbulence. Phys. Fluids 26 (10), 2769.

Straatman, A.G., Stubley, G.D., Raithby, G.D., 1998. Examination of diffusion modeling using zero-mean-shear turbulence. AAIA J. 36 (6), 929-935.

Ura, M., Komatsu, T., Matsunaga, N., 1987. Entrainment due to oscillating-grid turbulence in two-layered fluid. In: Chen, C.J., Chen, L.-D., Holly Jr., F.M. (Eds.), Turbulence Measurements and Flow Modeling. Hemisphere, New York, pp. 109-118.

Winterwerp, J.C., 1999. On the dynamics of high-concentrated mud suspensions. Ph.D. Thesis, Technische Universiteit, Delft.

Winterwerp, J.C., 2002. On the flocculation and settling velocity of estuarine mud. Cont. Shelf Res. 22, 1339-1360.