Modelling of sediment suspensions in oscillating grid turbulence

A k–ε 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 Rou0 and the product of the bulk Rouse number and the bulk Richardson number Ri0 Rou0, built on oscillating grid parameters. For Rou0 < 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 zm from the grid. The depth of the suspension layer is found to vary as zm/z0∝(Ri0 Rou0)−1/4. For Ri0 Rou0 ≪ 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 Rou0. For Rou0 > 0.01 the mean sediment concentration and turbulence decay simultaneously with increasing distance from the grid, and no sharp interface is observed.


Introduction
In many natural ows, the production of turbulence is due to the shear of the mean ow and is localised mainly at the bottom of the uid 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 Hopÿnger (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 e ect of turbulence di usion.
The former studies by E and Hopÿnger (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 signiÿcant 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 di cult to interpret because of the complexity of measuring the settling velocity of mud ocs. A striking property observed in cohesive sediment suspensions (Mory et al., 2002) was that the ux Richardson number, deÿned as the ratio of the buoyancy ux 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 di usive turbulence in stable stratiÿed uid (Briggs et al., 1998), it is generally considered that most of the mechanical energy input is dissipated viscously in stratiÿed ows (Hopÿnger and Linden, 1982) and, in particular, the ux Richardson number is below 0.25 in turbulent mean-shear ows (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 -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 di usive turbulence. There was a long debate in the past as to whether it is appropriate to use a k -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 -model were determined from experiments with turbulent mean-shear ows. This debate is connected to another very long one about the decay laws for turbulent kinetic energy in oscillating grid experiments (Hopÿnger 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 -model was used before by Sonin (1983) and Matsunaga et al. (1999) to describe purely di usive turbulence. When expressing the decay law in oscillating grid turbulence in the form of a power law (i.e. k˙z n ), one observes that the exponent n predicted by the k -model in clear water (with no particles) is di erent 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 -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 -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 modiÿes the decay of turbulence as compared with the decay in clear water, and to understand how this will produce lutocline formation. The k -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 speciÿcally interpreted in terms of cohesive sediments. The range of application of the present study is much concerned with ÿne-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 di cult for cohesive sediments as the settling velocity of sediment may vary due to occulation e ects. Modelling of occulation e ects 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 Hopÿnger (1987) and Huppert et al. (1995), which deal with non-cohesive particles.

2.
A steady k-model for particle suspensions

Formulation
We consider the case of non-cohesive particles which are maintained in suspension in a stirring tank by an oscillating grid. The conÿguration is sketched in Fig. 1. The vertical axis is oriented upwards, z = 0 is the grid mean position. The grid oscillation generates a turbulent ow. We assume here that there is no mean ow inside the tank. For constant stirring conditions, a steady-state particle suspension is established which may take the forms of two di erent typical vertical distributions of concentration. In the ÿrst case, sketched in Fig. 1a, which is basically obtained when the particle settling velocity is su ciently large, the particle settling ux is large and cannot be balanced by the turbulent buoyancy ux, 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 ux is much smaller than in the ÿrst case. Because the turbulent di usivity is high, the equilibrium of the settling ux with the turbulent buoyancy ux 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 ÿne-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 in the suspension layer is given by the mass balance equation where C(z) is the mean concentration at elevation z, and w and c are the uctuations of vertical velocity and sediment concentration. The averaged upward turbulent solid particle ux is modelled as The eddy di usivity t (z) is The eddy viscosity t (z) is related to the turbulent kinetic energy k(z) and the rate of dissipation of the turbulent kinetic energy (z) as The vertical distributions of k(z) and (z) are determined from the classical k -equations (e.g. Rodi, 1984;Winterwerp, 2002). The steady-state equations without mean ow are written below. The equation of conservation of the turbulent kinetic energy k is The equation of conservation of the dissipation rate of turbulent kinetic energy has a similar form where g is the acceleration of gravity, s = 2:65 kg=l and w = 1 kg=l are taken as the particle and water densities, respectively. We use the usual values for the k -model constants as k = 1; = 1:3; c = 0:09 and c 2 = 1:92: The values of the constants appearing in (3) and (6) to account for buoyancy e ects are not as usual. We take c = 0:7 and c 3 = 0: (8) Winterwerp (2002) suggested that c 3 = 1 for a stable stratiÿcation. 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 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 (z) are the solutions of with G = g( s − w )= s w . Eq. (9) expresses the momentum balance: the upward turbulent buoyancy ux balances the downward settling ux of particles. The spatial distribution of the turbulent kinetic energy is the result of three di erent e ects. We identify in (10) the di usion of the turbulent kinetic energy (D), the dissipation ( ) and the buoyancy e ect (B). A similar analysis holds for the distribution of the turbulent kinetic energy dissipation rate (Eq. (11)).
Introducing the local Rouse number and the local Richardson number we observe from (10) that the buoyancy e ect is negligible for the decay of the turbulent kinetic energy when Ri(z)Rou(z)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 In a similar way, a dimensional analysis of (9) shows that the mean concentration varies vertically over a distance of the order of which is much larger than Z k when Rou(z)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)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.

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 (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)1 for 0 ¡ z ¡ z 1 .
We refer in the present section to a paper by Matsunaga et al. (1999). Using the same k -model as we do, they determined the distributions of k(z) and (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 -model for modelling of purely di usive turbulence. On the other, their solution establishes the dependence of the distribution of k(z) and (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 -model to describe the transportation of turbulence by the sole e ect of di usion is questionable. This model was conceived by considering ows in which the distribution of turbulence inside the domain is governed by the production of turbulence from the mean shear ow and its constants were determined from experimental results obtained for such ows. The case of turbulence produced in an oscillating grid experiment is rather di erent. When the geometrical properties of the grid are suitably chosen, the mean ow 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 e ect of di usion. The ability of a k -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 -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 -model for the turbulent kinetic energy and the dissipation rate are −4:98 and −8:46, respectively, whereas the exponents suggested by several authors (Hopÿnger and Toly, 1976;Ura et al., 1987;DeSilva and Fernando, 1994) from their experiments were −2 and −4, respectively. This is a signiÿcant di erence from the mathematical point of view, but Matsunaga et al. showed that, on the basis of available experimental data, the results of the k -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 -model. Considering that Matsunaga et al. provided clear proof that the k -model is able to address oscillating grid turbulence, at least for the conditions of available experimental data, the k -model with its usual constant is used here. We additionally refer to the quantiÿcation of the turbulent kinetic energy k 0 and dissipation rate 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 = ¿ 5500 they obtained k 0 = 0:6f 2 S 2 (S=M ) 1=4 and 0 = 0:45f 3 S 2 S=M: We only consider this case here. Using (16) the local Rouse number and the local Richardson number 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 modiÿed 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 1. Introducing the following dimensionless quantities: and using a new independent variable Eqs. (9)-(11) of the k -model can be rewritten in a dimensionless form as 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, 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 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.

Numerical method and boundary conditions
A speciÿc 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 m ) = 0; k(z m ) = 0 and (z m ) = 0: In order to ÿnd 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 ÿnite di erence scheme. The integration is performed by starting from the lower boundary toward the lutocline for the boundary conditions (27) Integration of the k equation (5) over the vertical gives The consistent solution is obtained if the ux F(z) tends numerically to zero at z = z m . Close to the grid the decay of turbulence is modiÿed by the particle suspension only to a very limited extent as shown by Eqs. the results of the di erent cases computed showed that k ≈ Rou 0 Ri 0 and ≈ Rou 0 Ri 0 =c 2 . Fig. 2 shows the dependence of k and with Rou 0 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 di erent values of k and lead to non-valid solutions (for instance, C, k or are not decreasing in all the suspension with increasing distance from the grid, or one of the latter quantities becomes negative). The values of k and are adjusted just to ensure the upper boundary condition (28) at z = z m . The quantities k and must be very small because, for the conditions investigated, the turbulence in the vicinity of the grid is hardly modiÿed 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 k and with Rou 0 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 di erent behaviours. For all cases plotted, the grid oscillation parameters are the same. Di erent sediment concentrations C 0 and di erent settling velocities w s distinguish the three conditions. Three di erent bulk Rouse numbers and three di erent bulk Richardson numbers are considered, namely Rou 0 = 10 −3 , 10 −2 , 10 −1 and Ri 0 = 10 −1 , 1, 10. The di erent combinations of values of Rou 0 and Ri 0 makes the quantity Ri 0 Rou 0 vary from 10 −3 to 10 −1 . The vertical variations in the local Richardson number Ri(z), the local Rouse number Rou(z) and the local values of Ri(z) Rou(z) are also shown in the lower part of Fig. 3 in order to help the interpretation. The di erent cases plotted in Fig. 3 clearly show the e ect of the bulk Rouse number on the particles and turbulence distributions inside the suspension. The vertical concentration proÿles change considerably with the Rouse number Rou 0 . For the lowest values of the Rouse number (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 (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 proÿle. The bulk Rouse number Rou 0 therefore appears to be the parameter that determines the shape of the concentration proÿle and the conditions for which a homogeneous concentrated suspension is obtained. The vertical variations in k cw and cw in clear water are also shown in Fig. 3

(dotted line). Because
Ri 0 Rou 0 1 for all cases shown in Fig. 3, we observe that the decay of turbulence is not much a ected by the presence of particles in the lower part of the suspension layer (ẑ = z=z 0 ¡ẑ 1 = 0:5), as discussed in Section 2.2. The variations in k and k cw and the variations in and cw are almost superimposed in the suspension layer, respectively, except in the vicinity of the lutocline.
The vertical proÿles of the eddy viscosity 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 t ¡ 50 . For the cases considered here, t ¿ 50 as far as position z located less than 1 cm below the lutocline. At this level, the rapid decrease in k and associated with the occurrence of the lutocline is already settled. Using low Reynolds damping functions (e.g. Patel et al., 1985) would presumably not signiÿcantly modify the position of the lutocline. Additionally, a k -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 -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 signiÿcantly higher level of turbulence to predict a position of the lutocline signiÿcantly di erent from that given by the k -model.
For the same conditions as those of Figs. 3 and 4, the vertical proÿle of the turbulent kinetic energy ux F(z) (Eq. (29)) is ÿnally displayed in Fig. 5. We verify that this quantity vanishes at 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 ux at level z m was less than 10 −5 of the ux 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 su ciently small. The local ux Richardson number, deÿned as the ratio of the buoyancy e ects to the available turbulent energy in (10), provides an alternative estimate of the e ect of buoyancy on the decay of turbulence. Because D = B + , the ux Richardson number also takes the simple form as or is alternatively expressed in terms of the local Rouse number and of the local Richardson number as For the seven conditions already considered, Fig. 6 shows the vertical variations in the ux Richardson number. Depending on the conditions, di erent behaviours can be noticed. We clearly observe that the ux 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 proÿles of the ux Richardson number for computations with the bulk Rouse number Rou 0 = 10 −1 , which do not show the occurrence of a lutocline, display a smoother behaviour. As discussed in Section 1, the observation that the ux Richardson number may approach the maximum value of 1 in grid turbulence experiments is a striking property, which di ers strongly from observations in turbulence with mean shear ows. Larger values of the ux Richardson number imply that the dissipation rate is smaller as compared to the buoyancy term. As the energy is transmitted by di usive 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 e cient.

The depth of particle suspension layers
Determining the depth of the particle suspension layer as a function of the bulk Rouse number Rou 0 and the bulk Richardson number Ri 0 is of practical signiÿcance by providing a simple and direct quantiÿcation of experimental observations. Fig. 3 has already highlighted di erent behaviours 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 su ciently decreased due to dissipation and buoyancy e ects. The buoyancy e ect suddenly enhances the decay of turbulence and a lutocline is formed (as sketched in Fig. 1b). On the other hand, when Rou 0 is su ciently large, the concentration inside the suspension decreases simultaneously with the decay of turbulence. The buoyancy e ect 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 su ciently 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 e ect. However, when Rou 0 is su ciently 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 di erent 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 di erent data sets are identiÿed. The focus in this paper is on computation conditions with Rou 0 6 0:01. A lutocline is clearly identiÿed 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 deÿned in this case, because there is no sharp transition of density. It was determined, as explained in Section 2.3, as the ÿrst position where the concentration becomes zero or negative. This mathematical deÿnition 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 di erence with the case Rou 0 ¡ 0:01.

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 observed, for the same grid oscillation conditions, depending on the concentration level inside the suspension. Signiÿcant variations in the depth of the particle suspension layer were measured. Cohesive sediment suspensions are particularly di cult 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 ocs having di erent 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 deÿned 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 ux Richardson number reaches a value close to one. This was an unexpected result as it di ers signiÿcantly from previous experiments on suspensions in mean-shear turbulent ows. Ivey and Imberger (1991) pointed out that the ux 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 ux Richardson number is close to one.
Experiments with non-cohesive particle suspensions in an oscillating grid tank have been carried out by E and Hopÿnger (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 Hopÿnger (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 Hopÿnger, 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 ÿgures 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 di erent authors. The experiments by E and Hopÿnger 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-deÿned lutocline. The concentration decreases regularly with increasing distance from the grid (see the dashed lines in Fig. 3). which establishes that z m =z 0˙( Rou 0 Ri 0 ) −n , where the exponent n is 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 e ect 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 modiÿed analytical model, inspired by Huppert et al.'s model, is presented in Appendix B, which shows that the dependency z m =z 0˙( Rou 0 Ri 0 ) −1=4 is quite general. This analytical model di ers basically from the numerical model by the fact that the analytical model assumes that the integral lengthscale is not modiÿed in the suspension as compared to clear water. The comparison demonstrates that this crude assumption has no e ect on the ÿnal result. The more signiÿcant limitation of Huppert et al.'s model is that it can only address the case of suspensions with very low Rouse numbers.

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 di usive 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 deÿned 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 modiÿed version presented in Appendix B, but in Noh and Fernando's model the integral lengthscale is limited by the buoyancy lengthscale l b = √ k=N (where N is the Brunt-V ais al a frequency), which decreases rapidly to zero at a certain level that deÿnes the interface.
Our numerical computations highlight the respective e ects of the bulk Rouse number Rou 0 and the bulk Richardson number Ri 0 through the quantity Rou 0 Ri 0 . On the one hand, concentrated suspensions with a homogeneous mean concentration are observed when Rou 0 is below 0.01. When Rou 0 Ri 0 is small, the decay of turbulence is a ected 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 ux Richardson number is close to 1 in the vicinity of the interface. This conÿrms 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 0 Ri 0 . A dependency in the form z m =z 0( Rou 0 Ri 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 modiÿed version (Appendix B), when the decay laws suggested by Hopÿnger and Toly (1976) are used. The dependence of the analytical model is in the form z m =z 0˙( Rou 0 Ri 0 ) −0:21 for the decay laws proposed by Huppert et al. (1995). We observe that the di erence is small and in agreement with the conclusion of Matsunaga et al. (1999), which points out that the di erent models and decay laws for oscillating grid turbulence are not signiÿcantly di erent for the range of parameters usually considered.
The exponents p, q and r are real positive values because the three quantities k, and C should vanish at the lutocline. Introducing (A.1) into (9) We note that, for very small Z, the ÿrst 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 e ects. However, when using the standard model for fully developed turbulence the viscous di usive term (ÿrst term) dominates the turbulent dissipation (third term) when p ¿ − 1. We thus necessarily require the buoyancy term to balance the viscous di usive term and r = p − 2: (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 equation), there is no physical solution close to the lutocline (k and diverge). The exponents p, q and r are not determined at this stage but the k -equations indicate that buoyancy e ects dominate just below the interface.
Appendix B. Buoyancy ux 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 modiÿed 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 modiÿed 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. The notations used in this appendix are those of Huppert et al. in order to highlight the similarities and di erences 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 -model (Eq. (5)) are the same ( = u 3 =qz). Only the di usion terms in Eqs. (B.1) and (5) slightly di er. The main di erence between our analytical model and Huppert et al.'s model is that we do not assume that the buoyancy ux decreases linearly from the grid to the edge of the layer as stated by