Vertical grain size sorting in bedload transport on steep slopes with a coupled fluid-discrete element model

In order to study vertical grain size sorting in bedload sediment transport, numerical experiments of two-size particle mixtures were carried out, using a validated coupled fluid-discrete element model developed at Irstea. A 3D 10% steep domain, consisting at initial time of a given number of layers of 4 mm particles deposited on top of a coarser 6 mm particle bed, was sheared by a turbulent and supercritical fluid flow (Shields numbers of 0.1 and 0.3). The elevation of the centre of mass of the infiltrated fine particles is observed to follow the same logarithmic decrease with time, whatever the initial number of fine layers. This decrease is steeper for a higher Shields number. The main result is that this typical behaviour is related at first order to the particle shear rate depth profile.


Introduction
Bedload transport, the coarser bed material shear-driven by the turbulent water flow, has major consequences for public safety through flood alleviation, for the management of water resources, and for environmental sustainability.Bedload transport research has mainly been focused on mass flux and its correlation with the water flow.Yet after more than a century of work we have no satisfactory theory for sediment transport.Empirical sediment transport formulas often poorly compare with field measurements, by sometimes order of magnitudes [1].
An important reason for our limited ability to predict sediment flux is the absence of general understanding of grain-grain interactions in stream channels [2,3].
In mountains, steep slopes drive an intense transport of a wide range of grain sizes leading to grain size sorting also named size segregation [4].This phenomenon largely modifies fluxes and results in patterns that can be seen ubiquitously in nature.
Focusing on grain size sorting, the infiltration of fine sediment into an immobile coarse bed has been extensively studied in fluvial geomorphology because of ecological significance especially related to the reproduction of salmonids [5].However when the coarse bed is moving, another type of infiltration named 'kinetic sieving' can take place even with a coarse to fine diameter size ratio close to unity.Smaller grains may then move downward in the bed through voids that are opened by the general motion.This phenomenon which is important during floods has rarely been investigated in stream channels [6,7].
In industrial processes, usually quite energetic, kinetic sieving-type segregation takes place, and has been more largely studied.Indeed, to be able to perfectly mix input materials and thus to counteract the natural segregation remains a challenge.
At Irstea, a coupled fluid-discrete element model has been developed and validated to study bedload transport [8].In this paper we present a study of vertical segregation and kinetic sieving using the model which has been extended to simulate two grain sizes.Following classical theories such as [9], segregation rates are analysed in term of shear rates.

Coupled fluid-discrete element model
A two-phase Eulerian-Lagrangian flow model developed at Irstea is used in this research.It is composed of a three dimensional discrete element model (DEM), describing each individual particle coupled with a one dimensional Reynolds average Navier Stokes (RANS) model [10,11].
The 3D DEM uses the open source code Yade [12] with a classical spring-dashpot contact law.Hydrodynamic forces are restricted to the buoyancy and drag forces.The drag coefficient depends on the local particle Reynolds number taking into account hindrance effects.
The water flow is considered to be stationary and unidirectional in the streamwise direction so that the fluid profile depends only on the depth (1D).The RANS 1D model uses an integral mixing length closure [10,11].The effect of the particles on the fluid phase is taken into account through the drag force and the solid volume fraction.Both the effect of the particle on the fluid and of the fluid on the particle are considered (two-way coupling).The model has been validated with turbulent bedload transport experiments, in terms of both the sediment transport rate and granular depth profiles.More details can be found in [8,13].This model has also been used in turbulent bedload transport to analyse the dense granular rheology [14] and to revisit the effect of the slope on the vertical flow structure and the transport rate scaling [15].The idea is to contrast the case of one full layer of small beads (N=1) with the 'diluted' case of only a few scattered small particles The elevation of the water free-surface is imposed and creates a fluid flow on top of the granular bed, which is turbulent (bulk Reynolds number of about 10 4 ), hydraulically rough (particle Reynolds of order 10 3 ) and supercritical (Froude greater than 1).
Bi-periodic boundary conditions are imposed for the DEM, and only steady uniform conditions are analyzed.The unidirectional character of the fluid flow allows averaging spatially over the streamwise and spanwise directions, and analyzing the results only over the depth.
Both a low and a high dimensionless bottom fluid shear stress (Shields numbers Sh of 0.1 and 0.3) have been considered.The Shields number is defined as ( ). .with  s and  the glass and water density respectively.b  is unambiguously defined as the maximum of the fluid shear stress [8].Even for the Shields number of 0.3, no suspension occurred.
Due to kinetic sieving, the fine particles infiltrate through the coarser bed (Fig 1b).The centre of mass z c of the infiltrated particle cloud was calculated over time (output every 0.5s).It was necessary to simulate up to 10 4 seconds (3 hours) to be able to follow the infiltration of the fine particles down to the lowermost coarser layers over several time decades.Note that while the height of the coarse layers is taken as 8d for Sh=0.1, it is chosen as 12.d for Sh=0.3 because in the latter case, fine particles segregate and reach the bottom more rapidly.

Results
Both the fluid and particles starting from rest, the system follows a transient behaviour before reaching a characteristic steady state with a constant sediment transport rate.As one would expect, the infiltration evolution depends on the Shields number.From 50s to 5,000s the cloud of fine particles infiltrated through about 4.d coarse layers for Sh=0.1, and through about 6.d for the Sh=0.3 case.
At the beginning, the behaviour is transient, both the number of fine particle layer deposited N and the Shields number influence the evolution of the fine particles centre of mass (Fig. 2).During this period, the size-segregation phenomenon is coupled with the evolution of the system to the steady state, making a clear analysis not possible here.
Therefore, we rather focus on the steady state, which corresponds to the infiltration of fine particles through the dense coarse granular bed in quasi-static motion.In this region, the centre of mass position z c decreases linearly with ) log(t , provided the cloud is not yet affected by the bottom (Fig. 2).
This linear trend was fitted on all curves (Fig. 2).For a given Shields number, the slope was found to be in practice the same whatever N.For Sh=0.1 (Fig. 2a), the mean slope (out of the three N cases) amounts to -0.89 (min of -0.91 and max of -0.87).For Sh=0.3 (Fig. 2b), the mean is -1.21 (min of -1.22 and max of -1.20).Surprisingly the segregation rate does not depend of the concentration of fines, even for N= 0.05 corresponding to a very low number of fine particles initially scattered on top of the coarse bed, though the evolution is a bit less smooth than for higher N.As already mentioned, the segregation rate (given by the absolute slope) increases with the Shields number.Following theories such as Savage and Lun [9] (see also a simplified description in [16]), we will consider herein that the particle shear rate plays the major role.The segregation rate c dZ dt is taken proportional to the particle shear rate (upper case z and z c are normalized by the coarse diameter d).Streamwise coarse-only particle velocity depth profiles appear in Fig. 3 for both Shields numbers.An exponential decrease is observed in the quasi-static part of the bed [13] (see also inset on Fig. 3): According to this simple calculation, the fitted slope a appearing on Fig. 2 should be the same as that of the exponential decay of the particle velocity profile (precisely, it should be the inverse but note that it is z/d which is plotted in function of v p x as usually done to show the elevation profile).Velocity profiles were fitted over the exponential domain in Fig. 3a (Sh=0.1)and Fig. 3b (Sh=0.3)yielding mean values of respectively 0.80 (min of 0.79 and max of 0.84) and 1.10 (min of 1.08, max of 1.12).Those values compare well with the segregation evolution values (Fig. 2).For both Sh=0.1 and Sh=0.3, the difference of fitted slopes amounts to about 10%, giving credit to our assumption that the particle shear rate plays the major role.

Conclusion
After the initial transient period, the elevation of the centre of mass of the infiltrated fine particles is observed to follow the same logarithmic decrease with time, whatever the initial number of fine layers, from a few scattered fine particles up to one full fine layer.This decrease is steeper for a higher Shields number.
The main finding is that this typical behaviour is related at first order to the particle shear rate exponential depth profile present in the quasi-static region.These results confirm that the shear rate is the main driver of the kinetic sieving-type observed segregation in our bedload experiments.
While this is consistent with established theoretical framework for the full layer of fine particle, one would expect a different type of behaviour for isolated particles [17].Future work will therefore be devoted to the study of single fine particle behaviour in a similar system.
Advancing in our understanding of grain size sorting and upscaling from discrete element models to shallow water continuous models should ultimately improve bedload transport and river morphology modelling.

Fig. 1 .
Fig. 1.Infiltration of 1 layer of fine particle into the coarse bed a. t=0s b. t= 188s (Sh=0.1).Our 3D numerical domain is made of a 10% (5.71°) inclined channel filled with water, in which spherical particles are deposited.The coarse bed is made of particles of diameter d=6 mm, 15 times d wide and 15times d long.At initial time, a given number N of layers of 4 mm diameter particles (N = 0.05, 0.5, 1) are deposited on top of the coarse bed (Fig 1a).

Fig. 2 .
Fig. 2. Time evolution of the fine particles centre of mass for a number of layers N of 0.05, 0.5 and 1 , from light yellow (N=0.05bottom) to dark black (N=1 top).See text for the slope of the fit a. Shields of 0.1 b.Shields of 0.3 Fig. 2 shows the evolution over time of the centre of mass of the fine particles z c for both Shields numbers of 0.1 and 0.3 (the origin of z c is at the bottom of the domain).The initial centre of mass position z c is dependent of the number N of fine particle layers since the number of coarse layers is kept constant for a given Shields number.

Fig. 3 .
Fig. 3. Depth profiles of the streamwise particle velocity (m/s) with fitting in the exponential domain(See text for fit values) from light yellow (N=0.05bottom) to dark black (N=1 top).a. Shields of 0.1 b.Shields of 0.3.Inset shows the velocity profile with linear axes (N=0.05)