D Bonamy

F Daviaud

L Laurent

Experimental study of granular surface flows via a fast camera: a continuous description

Depth averaged conservation equations are written for granular surface flows. Their application to the study of steady surface flows in a rotating drum allows to find experimentally the constitutive relations needed to close these equations from measurements of the velocity profile in the flowing layer at the center of the drum and from the flowing layer thickness and the static/flowing boundary profiles. The velocity varies linearly with depth, with a gradient independent of both the flowing layer thickness and the static/flowing boundary local slope. The first two closure relations relating the flow rate and the momentum flux to the flowing layer thickness and the slope are then deduced. Measurements of the profile of the flowing layer thickness and the static/flowing boundary in the whole drum explicitly give the last relation concerning the force acting on the flowing layer. Finally, these closure relations are compared to existing continuous models of surface flows.

I. INTRODUCTION

Granular media present numbers of interesting and unusual properties [START_REF] Jaeger | Granular solids, liquids and gases[END_REF]: the intrinsic dissipative nature of the interactions between the constituent macroscopic particles sets granular matter apart from conventional gases, liquids or solids. One of the most interesting phenomenon in granular systems is the transition from a static equilibrium to a granular flow: contrary to ordinary fluids, they can remain static even with an inclined free surface. But when the angle of the surface exceeds some threshold value θ d , the pile can not sustain the steep surface and starts to flow until its angle relaxes under a given angle of repose θ r . The motion has the particularity to be a surface flow: most of the pile remains static. The condition governing the transition between the solid state and the liquid state as well as the internal equations of the flowing layer are still today under debate [START_REF] Mills | Model for a stationary dense granular flow along an inclined wall[END_REF]. This makes the description of an avalanche rather difficult.

People have benefited from the particularity of avalanches to be surface flows to apply to them depth averaged conservation laws, the Saint Venant equations [START_REF] De Barre Saint-Venant | Mémoires sur des formules nouvelles pour la solution des problèmes relatifs aux eaux courantes[END_REF], and to use the approximations developed in hydrodynamics for a thin fluid film flowing down an inclined plane. Savage and Hutter [START_REF] Savage | The motion of a finite mass of granular material down a rough incline[END_REF], and more recently Pouliquen [START_REF] Pouliquen | On the shape of granular fronts down rough inclined planes[END_REF] applied these equations to describe the motion of granular material down a rough inclined plane. Khakhar et al. [START_REF] Khakhar | Tranverse flow and mixing of granular materials in a rotating cylinder[END_REF], and T. Elperin and A. Vikhansky [START_REF] Elperin | Granular flow in a rotating cylindrical drum[END_REF] have studied steady flows in a rotating drum using the same theoretical framework. They succeed to report on the different free-surface shapes successively observed when the rotation speed increases. Finally, Douady et al. [START_REF] Douady | On granular surface flow equations[END_REF] and very recently Khakhar et al. [START_REF] Khakhar | Surface flow of granular materials: model and experiments in heap formation[END_REF] have used these equations to describe avalanches at the surface of a heap. Another model, referred to as the BCRE model [START_REF] Bouchaud | A model for the dynamics of sandpile surfaces[END_REF], has also been proposed a few years ago: it deals with thin surface granular flows including avalanches and succeeds for example to report on the existence of the two angles θ d and θ r . This model has then been extended to describe thick surface flows [START_REF] Boutreux | Surface flow of granular materials: a modified picture for thick avalanches[END_REF][START_REF] Aradian | Thick surface flows of granular materials: the effect of the velocity profile on the avalanche amplitude[END_REF].

All these models are based on the assumption that a strict separation between flowing grains and static grains can be made. They end in a system of two coupled equations with the same form, differing only through their closure relations. Those are chosen from theoretical considerations, linking for example the coupling term to the effective force applying on the flowing phase (for Saint-Venant approach) or to the microscopic effects of the flowing grains collisions on the static bed (for BCRE approach). An original approach, proposed by Aranson and Tsimring, should also be mentioned [START_REF] Aranson | Continuous description of avalanches in granular media[END_REF]. Assuming a continuous transition between flowing grains and static grains, they have developped a theory based on hydrodynamic equations coupled with an order parameter equation describing the solid/flowing transition. This theory accounts for most of the behaviours recently observed in experiments on rough inclined plane [START_REF] Daerr | Two types of avalanches behaviour in granular media[END_REF].

Thin rotating drums, i.e. rotating drums whose gap is small compared to the diameter, have been widely used to investigate properties of granular materials: Evesque and Rajchenbach [START_REF] Evesque | Characterization of glass bead avalanches by using the technique of a rotating cylinder[END_REF] and Jaeger et al. [START_REF] Jaeger | Relaxation at the Angle of Repose[END_REF] have studied avalanche size distribution. Rajchenbach has used it in an attempt to find constitutive mechanical laws in 2D beads packing [START_REF] Rajchenbach | Granular flows[END_REF]. Nakagawa et al. [START_REF] Nakagawa | Non Invasive measurements of granular flows by magnetic resonance imaging[END_REF] have used magnetic resonance imaging to measure density and velocity field in a thick rotating drum. We propose here to study steady surface flows in a rotating drum for a quasi-2D packing of monodisperse beads in order to find experimentally the closure relations needed to complete the depth averaged hydrodynamic equations.

In section II, basic conservation laws are presented. It is shown how the closure relations can be deduced from the experimental velocity profile inside the flowing layer at the center of the drum and from the flowing layer profile. The experimental set-up is then described (section III). Beads are studied on the microscopic scale in section IV: velocity and volume fraction profiles are measured on elementary slices via a fast camera. In section V, the flowing layer thickness and the angle between the interface and the horizontal are measured in the whole drum. The closure relations are deduced. Finally, they are compared in section VI to the ones assumed by the different models proposed in the litterature.

II. THEORETICAL FRAME: SAINT-VENANT DESCRIPTION

Let us consider the situation depicted fig. 1. The approach is limited to the 2D situation. The material can then be divided into two parts: a thin cascading layer flows on a static bed. The boundary between these two phases is supposed to be sharp. This assumption will be discussed and justified in section IV. As the flowing layer is assumed to be thin, the velocity of grains is essentially parallel to the static/flowing boundary. The most natural frame is consequently the curvilinear frame ( -→ e x , -→ e z ) where -→ e x (resp. -→ e z ) is locally tangent (resp. normal) to this interface. The z origin is set at the free surface. Calling (v x , v z ) the velocity at the (x, z) coordinate, the mass and momentum conservation equations can be written [START_REF] Douady | On granular surface flow equations[END_REF][START_REF] Khakhar | Surface flow of granular materials: model and experiments in heap formation[END_REF]:

∂ ∂t (< ρ > R) + ∂ ∂x (R < ρv x >) + (ρv z )| z=-R = 0 (1) ∂ ∂t (R < ρv x >) + ∂ ∂x (R < ρv x 2 >) + (ρv x v z )| z=-R = F ( 2
)
where R is the flowing layer thickness, ρ the solid fraction and F the x-component of the force acting on the volume of the flowing layer bounded by the two surfaces Σ(x) and Σ(x + dx). The quantities (R < ρv x >) and (R < ρv x 2 >) are respectively the flow rate and the flux of x-component of momentum through the surface Σ(x) (see fig. 1) while (ρv z )| z=-R and (ρv x v z )| z=-R represent respectively the flow rate and the flux of x-component of momentum passing through the static/flowing boundary. < A > denotes the average of the quantity A across the flowing layer :

< A >= 1 R 0 -R Adz (3)
In the following, the volume fraction ρ is supposed to be constant (this assumption will be discussed and justified in section IV) and eq. (1,2) are simplified accordingly. The expressions of < v x >, < v x 2 > and F as a fonction of the flow macroscopic parameters are required to close the two equations [START_REF] Jaeger | Granular solids, liquids and gases[END_REF][START_REF] Mills | Model for a stationary dense granular flow along an inclined wall[END_REF]. Translating the packing upwards should not change its dynamical equations. Consequently, the three closure relations can be written as < v x > (R, θ), < v x 2 > (R, θ) and F (R, θ) where θ(x, t) is the angle made by the solid/flowing interface with the horizontal. Let us apply these two conservation equations to the description of steady surface flows in a rotating drum (see fig. 1). In this case, time derivatives terms in eq. (1,2) vanish and the velocity v| z=-R of the static grains at the bed/layer interface can be written:

v| z=-R = Ωr n ( 4
)
where Ω is the rotation speed of the drum, r the distance between the center of the drum O and the point M of the interface where v| z=-R is calculated, and n the unity vector perpendicular to --→ OM. Then, eq. (1,2) become:

d dx (R < v x >) = -Ω X cos θ + H sin θ (5) d dx (R < v x 2 >) + Ω 2 X cos θ + H sin θ X sin θ -H cos θ .H X cos θ + H sin θ = F ρ (6)
where (X(x), H(x)) is the coordinate of M in the gravity frame whose origin coincides with the center of the drum (see fig. 1) and H the Heaviside function.

The aim of the paper is first to investigate the beads behaviour on the microscopic scale: the profiles of the velocity v x (z) and volume fraction ρ(z) are measured on an elementary slice. The evolution of these profiles with R and θ allows one to deduce the first two closure relations < v x > (R, θ) and < v x 2 > (R, θ). Then, profiles R(x) and θ(x) are measured on the macroscopic scale, in the whole drum. These profiles confirm the form of the first two closure relations, extend their range of validity and allow one to find the last closure relation F (R, θ).

III. EXPERIMENTAL SETUP

The experimental set-up is illustrated in fig. 2. It consists mainly in a duralumin rotating drum of diameter D 0 = 45 cm and of gap e = 7 mm. The drum is half filled by a packing of steel beads of diameter d = 3 ± 0.05 mm. The drum can then be considered as a Hele Shaw cell (D 0 ≫ e ≈ 2d) where the packing has a 3D microscopic structure (especially a disordered structure contrary to 2D packing whose stable structure is crystalline) while keeping a 2D macroscopic geometry. The volume fraction measured just after the filling is 0.57 ± 0.02. The gap e has been chosen slightly larger than 2d to prevent any jamming in beads flow. The use of steel beads of millimetric size allows to have a good control of geometrical and mechanical properties and to limit capillary and electrostatic effects at the expense of a larger scale experiment.

The rotation speed Ω can be varied from 0.001 r.p.m. up to 30 r.p.m. For small Ω -typically smaller than 0.1 r.p.m. -, the mean angle of the free surface θ oscillates between the angle of flow θ d = 27.5 ± 0.7 • and the angle of repose θ r = 23.5 ± 0.6 • . This intermittent behaviour is not studied in this paper. For Ω ≥ 0.1 r.p.m., surface flow becomes steady. The beads are lighted via a continuous halogen lamp and regimes obtained for Ω varying from 1 r.p.m. up to 6 r.p.m. are recorded via a fast camera. On this range of rotation speed, the width of the flowing layer R is well defined and inertial effects are negligible (the Froude number F r = Ω 2 D 0 /2g ≤ 0.01 where g is the gravity constant).

IV. MICROSCOPIC SCALE: VELOCITY AND VOLUME FRACTION PROFILE ON AN ELEMENTARY SLICE

Sequences of frames are recorded via a fast camera at a sampling rate f s = 1 kHz. The shutter speed is of 0.1 ms. The recorded region is located at the point where ( -→ v z | z=-R . -→ n ) = 0 (roughly corresponding to the center of the drum). In this region, θ, R and {v x (z)} are invariant under a small translation along -→ e x . The recorded region can then be considered as a statisfactory elementary slice. R and θ in this area are controlled by changing Ω, but can not been modified independently. The recorded region size in pixel is 480 × 234, one pixel corresponding to 0.227 mm. The digital images are processed to obtain the position of the center of mass of the beads seen through the glass porthole (see fig. 3). The image of a single bead is made up of about 23 ± 3 pixels depending on the distance of the bead to the porthole. The error on the determination of beads location is thus about 50 µm. One can then extract both volume fraction and velocity profiles.

Volume fraction profile can not be measured directly in our experiment since all beads are not seen through the porthole. One has thus calculated from numerical simulated packing the relationship between the real volume fraction and the volume fraction estimated from the beads seen through the porthole (see appendix A). The experimental volume fraction profiles can then be deduced. The errors made on these profiles is mainly dominated by the statistical processing of the data: profiles have been averaged over 1000 frames to limit this error. A typical volume fraction profile is represented fig. 4a: ρ varies typically from 0.3 to 0.4 in the flowing layer, except at the free surface in a small region (about three beads diameters) where ρ drops quickly down to zero. This last region is assumed to be sharp and corresponds to the free surface boundary. The assumption ρ ≃ constant can then be questioned: relaxing this assumption amounts to replace < v x > and < v x 2 > respectively by < ρv x > / < ρ > and < ρ(z)v x 2 > / < ρ >. As ρ variations are negligible compared to v x variations, ρ can be assumed to be constant.

Frames processing allows to evaluate beads velocity (see appendix B). A typical velocity profile is represented fig. 4b. The velocity profile is linear in the whole flowing layer. The inset, representing the variation of the velocity vs. depth in semilogarithmic scale, emphasizes the creep motion of the static phase. In this phase, v x decreases exponentially with the depth with a characteristic decay length λ creep = 7.6 mm ≃ 2.5d in agreement with [START_REF] Komatsu | Creep motion in a granular pile exhibiting steady surface flow[END_REF]. The creep motion is thus localised in a narrow layer at the solid/liquid interface whose thickness can be considered as zero at macroscopic scale. The creep motion is not studied in this paper.

The rotation speed Ω has been varied from 1 r.p.m. to 6 r.p.m. in order to study variations of the velocity profile with R and θ. The corresponding velocity profiles have been represented fig. 5a. In fig. 5b, all profiles have been translated along the z-axis to make the static/flowing interface coincide. The velocity profile remains linear in the whole flowing layer and the velocity gradient is independent of both R and θ.

Velocity profiles measured at five different points equally distributed along the free surface for Ω = 6 r.p.m. are represented in fig. 6. At these locations, d x R and d x θ are no more equal to zero. The thickness of the recorded region has been decreased to minimize these drifts. {v x (z)} does not depend on the location in the drum, which means that {v x (z)} does not depend on d x R and d x θ. The velocity profile can be then written inside the flowing layer as:

v x (z) = γ(z + R) (7)
with γ = 34 ± 0.5 s -1 , or γ ≃ 0.6 g/d. Such order of magnitude can be easily understood [START_REF] Rajchenbach | Granular flows[END_REF]: when a bead collides with an underlying grain, all the kinetic energy is absorbed by multiple collisions. Then, the balance between potential and kinetic energy on a mean free path of order of d leads to a limiting velocity between two adjacent beads layers equal to √ gd, and consequently to a velocity gradient of order of g/d. The knowledge of the profile {v x (z)}(R, θ) enables us to calculate the first two closure relations:

< v x > (R, θ) = 1 2 γR (8) < v x 2 > (R, θ) = 1 3 γ2 R 2 (9)
This scaling law is different from the one observed in granular flows down rough inclined plane for 2D packing [START_REF] Azanza | Ecoulements granulaires bidimmensionnels sur un plan incliné[END_REF] and for 3D packing [START_REF] Pouliquen | Scaling laws in granular flows down rough inclined plane[END_REF]: In these experiments, < v x > was found to scale as R 3/2 . A linear profile with a velocity gradient independent of R has also been observed by Rajchenbach at the surface of a 2D packing in a rotating drum [START_REF] Rajchenbach | Granular flows[END_REF]. This tends to confirm that boundaries conditions can strongly affect the velocity profile and questions the ability of a classical constitutive law relating the shear stress to the strain rate [START_REF] Mills | Model for a stationary dense granular flow along an inclined wall[END_REF][START_REF] Rajchenbach | Granular flows[END_REF].

Let us note that velocity gradient observed by Rajchenbach in 2D packing scales as √ sin θ contrary to ours. However, in our series of experiments, θ varies only from 31 • up to 48 • , which corresponds to variations of √ sin θ from 0.72 to 0.86. Consequently, the difference between γ behaving as √ sin θ and γ constant does not differ much from the error bars and it is not possible to conclude at this point which of the two scaling is the right one. This will be tested on a more significant range of θ, from 0 • up to 48 • , in the next section.

V. MACROSCOPIC SCALE: PROFILES OF THE FLOWING PHASE HEIGHT AND OF THE STATIC/FLOWING BOUNDARY

The determination of the last closure relation F (R, θ) requires the determination of R and θ in the whole drum. The experimental protocol is the following: sequences of 512 frames are recorded via a fast camera at a sampling rate f s = 250 Hz with a pixel resolution of 480 × 400. In order to keep a sufficient resolution on the profiles, the rotating drum is visualized by halves, a pixel corresponding to 0.468 mm. Image processing allows to isolate the flowing layer and both the interface profile H(X) and the free surface profile S(X), expressed in the gravity frame ( -→ e X , -→ e Z ), can be extracted (cf fig. 7). The fluctuations on each point of these profiles are around 1.5 mm. The profiles are then averaged over 20 frames. Consequently, the resulting errors on both H and S are around 0.3 mm. The curvilinear coordinate x(X), the local angle θ(X), and the flowing layer thickness R(X) are then deduced using the procedure given in appendix C. Resulting errors on R is about 1 mm (see fig. 8b). On the upper part of the profile, the local slope is constant on a range whose size decreases when Ω is increased (see fig. 8a). One can study on this range L Ω the variation of the different variables with R for a fixed angle θ Ω depending on the rotation speed Ω.

Let us first confirm the first closure relation (eq. 8) established at the end of the last section: the drift velocity < v x > can be calculated from the experimental x(X), H(X) and R(X) profiles by integrating numerically eq. 5. Fig. 9a shows the variation of < v x > vs. 1 2 γR with γ = 34 s -1 as suggested in the previous section for different θ Ω . This confirms the assumption of a linear profile {v x (z)} with a velocity gradient independent of R.

Let us answer now to the question about the dependence of γ with θ. To test the different assumptions on this scaling, namely if γ scales as √ sin θ or is independent of θ, profiles H(X) and R(X) have to be taken into account in their whole: Typically, for Ω = 6 r.p.m., θ(X) varies from 0 • to 48 • . The procedure will be the following: from the profile H exp (X) obtained experimentally, one deduces X(x), H(x), θ(x) and, by integrating the eq. 5, the flow rate {R < v x >}(x). The theoretical profiles for the flowing thickness R (a) R < v x >= pR 2 corresponding to a velocity gradient independent of θ

(b) R < v x >= pR 2 √
sin θ corresponding to a velocity gradient γ scaling as √ sin θ as observed by Rajchenbach [START_REF] Rajchenbach | Granular flows[END_REF] These two profiles are then compared to the experimental profile R exp (x) (see fig. 9b). The fit parameter is p. The scenario (a) leads to better fits, with a fit parameter p = 14 s -1 about 20% smaller than γ/2. This difference may come from the boundaries: the velocity gradient γ has been determined from the velocity of only the beads observed through the porthole while p is related to the mean velocity of all the beads.

At this point, one can evaluate the force F . Numerical calculation of the second left-hand term of eq. 6 shows that it is negligible compared to the first left-hand term and can thus be neglected. From eq. 5, the simplified eq. 6 and the closure relations 8 and 9, one finds:

F ρ = γRΩ H sin θ + X cos θ (10)
The last closure relation can then be determined. From the experimental X(x), H(x) and R(x) profiles, one can deduce F (x) using eq. 10. Fig. 10a presents the variation of F/ρ γR vs. R for different θ. F/ρ γR seems to decrease as a parabola with R:

F ρ γR = F 1 (θ) -F 2 (θ)R 2 (11)
A minimal form of the force F should include:

• the weight F W = ρgR sin θ

• a pressure term F P = -kρg dR dx cos θ [START_REF] Pouliquen | On the shape of granular fronts down rough inclined planes[END_REF]. Since R is negligible compared to D 0 , this term can be neglected. • the friction on the static phase: F S.f r = -µρgR cos θ. A priori, µ depends on R and θ. Experiments shows however that µ reach quickly a constant value close to tan θ r independent of R and θ as soon as R is larger to few beads diameters [START_REF] Daerr | Dynamique des Avalanches[END_REF]. In our case, µ can be considered as constant.

• finally, the friction on the lateral boundaries F B.f r = R 0 σ B.f r dz. This last contribution is the one responsible for the decrease of F/ρ γR with R as a parabola. It can be obtained directly from the value of F and the first two contributions F W and F S.f r . Fig 10b represents the variation of F B.f r /ρ vs. R for different θ Ω . Its variation with θ Ω appears to be weak. The form of F B.f r can be understood through the following process:

F B.f r = R 0 σ B.f r dz where σ B.f r (z)
is the stress exerted by the boundary on the layer between z and z + dz. To model σ f r.B , let us apply the following process: σ f r.B = νδp where ν is the frequency of beadboundary collision and δp the x-momentum lost by each collision. As ρ is constant, the mean free path is constant in the layer and ν ∝ ρv x (z). Each collision with the boundary is inelastic. Consequently, δp ∝ v x . Finally one can write σ f r.B = ηρv x 2 with η constant and, after integration over the width of the flowing layer:

F f r.B = -ρ η γ2 3 R 3 (12)
This closure relation indeed allows one to reproduce the experimental F B.f r (R) profiles obtained for different θ Ω (cf. fig 10b). In our experiment, η ≃ 10 m -1 . The order of magnitude of η can be understood in the following way: the momentum lost by each collision is of order (1 -e) where e is the restitution coefficient. For our beads, e ≃ 0.95. The mean free path is of the same order of the bead diameter d. Consequently, one expects η ≃ (1 -e)/d ≃ 10 m -1 .

The last closure relation that follows from the above discussion thus reads:

F = ρgR sin θ -µρgR cos θ - 1 3 η γ2 R 3 (13)

VI. DISCUSSION

Using the closing relations eq. 8, eq. 9 and eq. 13, the conservation equations (1,2) become:

∂R ∂t + V ∂R ∂x = Γ (14)
∂H ∂t = -Γ (15)
where:

V = 2 < v x >= γR (16)
and:

Γ = F ρ γR = g γ (sin θ -µ cos θ) - 1 3 η γR 2 (17)
This system of equations contains three parameters: the velocity gradient γ, the friction coefficient µ and a third parameter η characterizing the friction on lateral boundaries. These three parameters depend on microscopic features of the beads such as their diameter, their coefficient of restitution and their roughness and of the gap of the drum. It can be compared to the different existent continuous models proposed in [START_REF] Douady | On granular surface flow equations[END_REF][START_REF] Khakhar | Surface flow of granular materials: model and experiments in heap formation[END_REF][START_REF] Bouchaud | A model for the dynamics of sandpile surfaces[END_REF][START_REF] Boutreux | Surface flow of granular materials: a modified picture for thick avalanches[END_REF][START_REF] Aradian | Thick surface flows of granular materials: the effect of the velocity profile on the avalanche amplitude[END_REF]. These two equations are very close to the one proposed by Douady et al. [START_REF] Douady | On granular surface flow equations[END_REF]. They have supposed a linear velocity profile with a constant velocity gradient as observed in our experiment. Moreover, for large R, variation of the force F with R and θ is qualitatively similar to ours, but they attribute this variation to the friction of the flowing layer to static grains instead of invoking friction effects on the lateral boundaries as we have done.

Khakhar et al. [START_REF] Khakhar | Surface flow of granular materials: model and experiments in heap formation[END_REF] have proposed a derivation of depth averaged conservation equations very close to that of [START_REF] Douady | On granular surface flow equations[END_REF]. They also have assumed a linear velocity, but with a gradient dependent of θ. As for the closing relations proposed by [START_REF] Douady | On granular surface flow equations[END_REF], variation of F with R is attributed to a dependency of the shear stress at the solid/flowing interface with R rather than to boundary effects.

Another model has also been proposed a few years ago by Bouchaud et al [START_REF] Bouchaud | A model for the dynamics of sandpile surfaces[END_REF]. The original model, devoted to describe avalanche triggering, deals with thin moving layers of R ≤ d and can not be compared to ours. But this model has later been extended to describe thick surface flows by Boutreux et al. [START_REF] Boutreux | Surface flow of granular materials: a modified picture for thick avalanches[END_REF]. These models assume a constant advection velocity V , small variations of θ around θ r , and relate the coupling term Γ to the effect of collisions between the moving grains on the static bed. For thin flowing layers [START_REF] Bouchaud | A model for the dynamics of sandpile surfaces[END_REF], all the moving grains are assumed to be in contact with the static bed and to have the same probability to interact with it: Γ = aR(θ -θ r ) whith a constant. For thick granular flows [START_REF] Boutreux | Surface flow of granular materials: a modified picture for thick avalanches[END_REF], only grains belonging to the first layers can interact with the static bed. The influence of the others is shielded by these first layers. Γ is then more naturally assumed to be independent of R: Γ = a(θ -θ r ) with a constant. This last form is compatible with our form of Γ expanded to first order in θ -θ r . The constant a can be evaluated [START_REF] Douady | On granular surface flow equations[END_REF][START_REF] Khakhar | Surface flow of granular materials: model and experiments in heap formation[END_REF]: a = g cos θ r / γ. Thanks to their simplicity, these two last model allow analytical exact solutions [START_REF] Mahadevan | Propagating fronts on sandpile surfaces[END_REF][START_REF] Dorogovtsev | How sandpiles spill: sandpiles problems in a thick flow regim[END_REF][START_REF] Thorsten | Exact solutions of a model for granular avalanches[END_REF]. But the assumption of a constant velocity V prevents them for describing quantitatively avalanches front [START_REF] Douady | The four fronts and the two avalanches[END_REF] or avalanches amplitudes [START_REF] Aradian | Thick surface flows of granular materials: the effect of the velocity profile on the avalanche amplitude[END_REF].

VII. CONCLUSION

Depth averaged conservation equations have been written for granular surface flows. Their applications to the study of steady surface flows in a rotating drum allow one to find experimentally the closure relations of these equations from measurements of the velocity profile in the flowing layer at the center of the drum and from the flowing layer thickness R and the interface slope θ profiles. These measurements have been performed for a quasi 2D packing of steel beads. The velocity profile has been found linear with a gradient γ independent of both R and θ. The closure relations have been deduced: the mean drift velocity < v x > and the mean kinetic energy < v x 2 > are given respectively by < v x >= 1 2 γR and R < v x 2 >= 1 3 γ2 R 2 . The last relation relating the force F acting on the flowing layer to R and θ has also been explicitly given and reveals the importance of lateral boundaries. Two parameters, namely the friction coefficient µ of the flowing layer on the static bed (close to tan θ r ) and a coefficient η characterizing the friction on lateral boundaries are needed to describe F . Finally, the closed equations have been compared to existent continous models. Their differences have been discussed: a new term whose form is given should be included to take into account lateral boundaries effects. The precise study of the dependence of the three parameters γ, µ and η with the beads diameters, the beads coefficient of restitution, the beads surface roughness and the cell gap represents interesting topic for a future investigation.

VIII. APPENDIX

A. Determination of the volume fraction profile 3D disordered packing of different volume fractions confined between two plates have been simulated. Numerical packing have been constructed by dropping randomly beads in a cell of gap e. Beads are randomly located at a given altitude. Their altitudes are then minimized under the two following constrains: the beads should stay in the cell and the distance between their center of mass and the center of mass of beads already in the cell should be superior to d. Confined packing of volume fraction 0.5 can then be obtained. Confined packing of lower volume fraction are then obtained either by contracting beads and cell dimensions without changing the position of beads center of mass or by removing randomly a given number of beads. One can then determine the beads actually observed through the porthole, i.e. the beads whose center of mass is not hidden by the the presence of the other beads closer to the porthole. Both real volume fraction ρ real and the volume fraction ρ view evaluated from the observed beads are calculated (see fig. 11). ρ real varies linearly with ρ view . For d = 3 mm and e = 7 mm:

ρ real = 1.61ρ view (18)
Direct experimental measurements of ρ real and ρ view just after the filling the drum, i.e. for ρ real = 0.57 (see section III), confirm the value of the ratio ρ real /ρ view .

B. Determination of the velocity of beads

The processing to determine beads velocity is the following: each bead is tracked over ten successive frames. As soon as the displacement between two successive frames is larger than d, there is an ambiguity and the bead is considered as lost. Then, only beads tracked over at least five successive frames are kept. They represent more than 99% of all the observed beads. The largest accessible velocity is thus f s d/2 = 1500 mm/s and the error on these velocities is about 10 mm/s. The local angle θ of the surface flow with the horizontal axis is then determined from the mean velocity of the flow calculated by averaging the velocity of all beads in all frames. All frames are then rotated by θ and divided into horizontal layers one bead diameter wide. The velocity profile is then calculated by averaging the velocity of the beads in each layer on 200 frames. The errors made on the velocity profile is then mainly dominated by the statistical processing of the data. It is less than 20 mm/s for a 99% confidence interval.

C. Determination of x, θ and R

From experimental H(X) and S(X) profiles, x(X) is calculated via x = XP X 1 + dH dX (y) 2 dy where P is the upper boundary of the flowing layer (H(X P ) = S(X P )). The local angle θ(X) is defined as: θ = atan(dH/dX). The procedure to determine R(x) is the following: for a given x, one first finds the corresponding static/flowing boundary point P H of coordinate (X(x), H(x)) in gravity frame. Then one looks for the point P S of coordinate (X S , S(x S )) at the free surface so that ---→ P S P H be orthogonal to -→ e x (numerically, this consists of minimising (X(x) -X S ) cos θ(x) + (H(x) -S(X S )) sin θ(x) with respect to X S ). R(x) is then defined as |P S P H |.

  th (x) are then calculated assuming respectively the relations:

FIG. 11 .

 11 FIG. 1. Bonamy, Phys. Fluids

ACKNOWLEDGMENTS

We thanks J.-Ph Bouchaud for a critical reading of this manuscript. We are grateful to B. Faucherand and M. Planelles for their participation to the data collection. We also acknowledge discussions with M. Bonetti, P. Claudin, A. Daerr, S. Douady, P. Evesque, D. Khakhar and O. Pouliquen. We thank C. Gasquet and P. Meininger for technical support.

• Fig. 1: steady surface flow in a drum at rotation speed Ω. A thin cascading layer of thickness R (in grey) flows on a static bed whose local slope with the horizontal is θ. The transport theorem is applied on the volume limited by the vertical slices Σ(x) and Σ(x + dx) and resulting equations are written in the curvilinear frame ( e x , e z ). The frame ( e X , e Z ) is the gravity frame located at the center of the drum O. The solid rotation at a point M located at the solid/flowing interface is Ωr n where r is the distance between O and M and n the unit vector perpendicular to --→ OM.

• Fig. 2: schematic drawing of the experimental set-up

• Fig. 3: determination of beads velocities. Left: one of the 200 raw frames for Ω = 3 r.p.m.. Frames are binarised and beads position is then defined as the center of mass of white areas. Right: corresponding calculated velocity field.

• Fig. 4: left: solid fraction profile (averaged over 1000 frames). Right: velocity profile (averaged over 200 frames) at the center of the rotating drum for Ω = 6 r.p.m.. v x is linear with z in the flowing layer. Inset:

vs. the depth showing the exponential creep deformation in the static phase. The straigth line corresponds to an exponential fit. In both figures, the errorbars correspond to a 99% confidence interval.

• Fig. 5: left: velocity profile at the center of the drum for different Ω: for all these (R, θ) couples, v x is linear with z. Right: all these profiles have been translated along the z axis to make their static/flowing interface coincide. The velocity gradient γ is independent of both R and θ.

• Fig. 6: velocity profile for Ω = 6 r.p.m. at five different locations equally distributed along the free surface. The velocity gradient depends weakly on the location in the drum. All these profiles have been tranlated along the z axis to make their static/flowing interface coincide. The velocity gradient γ is consequently independent of both d x R and d x θ.

• Fig. 7: frame processing. Left: one of the 512 raw frames (Ω = 3 r.p.m.). Frames are subtracted 2 by 2, binarised and gathered by bursts of 25 frames. The resulting frame is smoothed via a dilation-erosion filter. Right: one of the 20 resulting frames showing the flowing layer.

• Fig. 8: typical profiles obtained for Ω = 3 r.p.m.. Left : profiles H(X) and S(X). Dashed circle correspond to rotating drum limits. One can see a large range L Ω in the right part of the drum where the layer thickness varies while the static/flowing boundary slope almost does not change. Right : profile R(x).

• Fig. 9: left: mean drift velocity < v x > vs. R. The solid line corresponds to eq.12 (translated vertically of a constant value.): η ≃ 10 m -1 .

• Fig. 11: variation of the real volume fraction profile ρ real vs. the volume fraction profile ρ view taking into account only the beads effectively seen through the porthole. Data come from numerical simulations.


D Bonamy

F Daviaud

L Laurent

M Bonetti

J P Bouchaud

Multiscale Clustering in Granular Surface Flows

Keywords: numbers: 45.70.Mg, 05.45.Df, 83.50. -v, 45.70.Qj

We investigate steady granular surface flows in a rotating drum and demonstrate the existence of rigid clusters of grains embedded in the flowing layer. These clusters appear to be fractal and their size is power law distributed from the grain size scale up to the thickness of the flowing layer. The implications of the absence of a characteristic length scale on available theoretical models of dense granular flows are discussed. Finally, we suggest a possible explanation of the difference between velocity profiles observed in surface flows and in flows down a rough inclined plane.

Granular materials share properties with both usual liquids and solids. They can form an inclined free surface without flowing, but when the angle of the free surface exceeds some threshold value, an avalanche occurs. Global behaviors can be described by models derived from fluid mechanics [1,2] or nonlinear physics [3,4]. However, some experimental results remain far from being understood. Detailed measurements of the mean density profile and the mean tangential velocity profile obtained in threedimensional flow (3D) by NMR [5] and in two-dimensional flows (2D) by direct image analysis [6,7] show strong evidence that the relation between stress and strain is nonlocal: (i) the velocity gradient is found to be constant in the flowing layer whereas momentum balance predicts a linear variation of the shear stress with depth [8]; (ii) the velocity gradient presents a different scaling with the depth for dense granular flows down a rough inclined plane [9] indicating the nonlocal influence of boundary conditions on internal rheology inside the flowing layer; (iii) the velocity gradient does not vanish at the free surface at variance with typical fluids. Very recent models have tried to account for some nonlocal effects [6,8,10]. We report here the first experimental evidence of rigid clusters of grains embedded in the flow and characterize their geometrical and statistical properties. Although clustering instabilities driven by the inelasticity of grain collision are well known in granular gases [START_REF] Goldhirsch | Physics of Dry Granular Media[END_REF], they have never been observed in dense surface flows. We find that these clusters are fractal and their size is power law distributed from the microscopic scale -the diameter of a grain -up to the macroscopic scale -the flowing layer thickness. Therefore, no characteristic correlation length can be defined in the flowing layer. Almost 50% of flowing beads belong to these objects. We discuss recently proposed nonlocal models in the light of these experimental observations and propose an argument to explain the differences observed between granular surface flows and granular flow down a rough inclined plane.

Experimental setup.-The experimental setup is illustrated in Fig. 1a. It consists in a duralumin rotating drum of diameter D 0 45 cm and variable gap, half filled with steel beads of diameter d 3 6 0.05 mm. Two drum thicknesses are used:

(A) A gap of 7 mm so that a quasi-2D packing is obtained but with a local 3D microscopic disordered structure. The fast camera allows us then to track the beads actually seen through the transparent sidewall of the tumbler: Around 40% of the beads are hidden and cannot be tracked.

(B) A gap of 3 mm so that a pure 2D packing is obtained and all the beads can be tracked. Additional metallic obstacles (bound to the inside of the drum) have been added in the static part of the packing to prevent any 2D ordering effects which would induce nongeneric effects. The main drawback is that the measurements can be made only during the time when the obstacles are in the bottom part of the drum.

In the recorded region, located at the center of the drum, the flowing layer thickness R and the mean angle u of the flow can be changed by varying the rotation speed V. Regimes obtained for V varying from 1 to 8 rpm are investigated. On this range of rotating speeds, surface flows are steady and inertial effects are negligible (the Froude number Fr V 2 D 0 ͞2g, g being the gravity constant, varies from 2.5 3 10 24 to 1.6 3 10 22 when V varies from 1 to 8 rpm).

The beads are lighted via a continuous halogen lamp. Sequences of 200 frames are recorded via a fast camera at a sampling rate of 1 kHz. The recorded region size in pixels is 480 3 234, one pixel corresponding to 0.227 3 0.227 mm. Frame processing allows us to obtain the position of the center of mass of the beads seen through the transparent side wall of the tumbler. Since the image of a single bead is made up of about 23 6 3 pixels (depending on the distance to the sidewall), the errors on the determination of bead location is about 150 mm [START_REF] Moukarzel | For a pure 2D packing made up of weakly polydisperse, rigid smooth beads, the coordination of each bead of the bulk is four[END_REF]. Tracking each bead on ten successive frames allows us to evaluate its velocity averaged on 10 ms (see Fig. 1b). The mean angle u of the flow is calculated from the averaging of the velocity of all the beads on the whole sequence of 200 frames. The images are then divided into layers velocity V m ͑z͒ of the layer at depth z is then defined as the average of the velocity of all the beads of the sequence of 200 frames whose center of mass is inside the layer.

For both geometries (A) and (B), two phases can be observed: a solid phase experiencing a creep motion where V m decreases exponentially with z (see also [13]) and a flowing phase exhibiting a clear linear velocity profile (see Fig. 1c for the quasi-2D case) with a velocity gradient ᠨ g dV m ͞dz independent of R, u, the coefficient of restitution of the beads and drum gap [7]. The creep motion is not studied in this Letter. In the following, we focus on the flowing layer and neglect beads with velocities less than 100 mm͞s since these are assumed to belong to the static bed.

Velocity fluctuations in quasi-2D and pure 2D flows.-In order to investigate possible collective effects, we study first the spatial correlations of the instantaneous velocity field for both geometry (A) and (B). Calling V͑x, z, t͒ the velocity of a bead of the frame t located at the ͑x, z͒ coordinate in ͑e x , e z ͒, where the unit vector e x (respectively, e z ) is parallel (respectively, perpendicular) to the mean flow, the bead fluctuation velocity is defined as Ṽ͑x, z, t͒ V͑x, z, t͒ 2 V m ͑z͒e x . Since the orientation of the fluctuation ñ Ṽ͞ Ṽ is not correlated to the depth z (whereas the amplitude Ṽ of the fluctuation decreases with the depth z), we focused on the instantaneous field made up of these orientations.

One of these instantaneous fields is represented in Fig. 2a and reveals aggregates of strongly correlated beads. These correlations can be quantified by looking at the orientation-orientation correlation functions: Given

V 1 rpm, V 2 rpm, V 3 rpm, V 4 rpm, V 6
rpm, and V 8 rpm. For these different V, i.e., for different u and R, the distribution decreases as a power law (straight line) with the same exponent a ͑A͒ 2.9 6 0.1.

two particles, labeled by 1 and 2, let k be the unit vector pointing from the center of 1 to the center of 2. The vector ñ1 has component ñk 1 parallel to and ñЌ 1 perpendicular to k; it is likewise for particle 2. Following [14], we define two correlation functions C k ͑r͒ P ñk 1 ñk 2 and C Ќ ͑r͒ P ñЌ 1 ñЌ 2 where the sums are over all particles such that the distance between the two particles is within d͞4 of r. These two correlation functions reveal strong long range correlations (of the order of 30%) that extend on a range of the same order as the flowing layer thickness R (see Fig. 2b). To isolate the correlated aggregates, we have used the following criteria: Two beads i and j belong to the same cluster whenever (i) they are in contact [START_REF]difficulties: two beads are in contact as soon as the distance separating their center of mass is smaller than 1.1d. The factor 1.1 accounts for the precision in the experimental location of the beads[END_REF] and (ii) the angle d ij between ñi and ñj is smaller than a given value d c arbitrarily chosen equal to 60 ± . This 034301-2 034301-2 value corresponds to the point where the experimental orientation distribution crosses the value of the uniform distribution (see Fig. 2c). Values ranging from 30 ± up to 90 ± have also been tried and do not modify the following conclusions. Figures 2d and2e show typical frames of the individual clusters for geometry (A), respectively, for V 8 rpm and V 1 rpm. Note that the size of the largest ones is comparable to the flowing layer thickness R: The distribution of the number of beads N in a cluster is a power-law N 2a in both geometries (A) and (B), with an exponent equal, respectively, to a ͑A͒ 2.9 6 0.1 and a ͑B͒ 2.0 6 0.1 for quasi-2D and 2D flows (see Fig. 2f for the quasi-2D case). In both cases, the exponent is independent of V, i.e., independent of both R and u. This indicates that no typical correlation length scale can be defined to characterize the spatial correlations of the velocity field.

Volume fraction fluctuations in pure 2D flows.-In order to find out the physical origin of these correlated clusters, measurements on the local volume fraction fluctuations were performed on the 2D stacking [geometry (B)]. A Voronoi tesselation can then be used to define the local volume fraction n associated with each bead. Let us note that this technique assumes that all beads' positions are known and consequently cannot be applied for the quasi-2D packing of geometry (A). The relative fluctuations of n are very small, a few percent (see Fig. 3a), but sufficient to alter significantly the beads' behavior. According to the "Reynolds" dilatancy concept [START_REF] Reynolds | [END_REF], an assembly of rigid particles cannot deform whenever its volume fraction is higher than a given threshold n c . For pure 2D packing of monodisperse beads, n c p͞4. Clusters of beads in contact with n $ n c have been isolated in the flowing layer (see Figs. 3d-3f). Almost 50% of beads of the flowing layer belong to one of these solid clusters. As was found for the clusters defined using velocity correlation, the size of largest ones is of the order of R. The size distribution is again a power law with an exponent a ͑B͒ n Ӎ 1.5, weakly dependent of V (see Fig. 3c). This exponent is slightly smaller than the one obtained from the analysis of velocity fluctuations. It should be noted that although there is some overlap between both types of clusters, they are not identical. We have also plotted the number of beads N in a cluster against the cluster radius of gyration R g (Fig. 3b). The power-law scaling is indicative of a fractal structure with an apparent fractal dimension equal to d ͑B͒ n 1.7 6 0.1. The mean aspect ratio of the clusters, calculated as the ratio of the two principal axes of their inertia tensor, does not depend on their size and is found to be close to 2. They preferentially orient along the flow direction. Their lifetime, defined as the duration during which their size is larger than the half of their maximal size, is of order of p d͞g Ӎ 10 22 s, where g is the gravity. Discussion.-These solid clusters can be compared to the collective objects postulated in different nonlocal models for granular flows. Mills et al. [8] postulate the exis- tence of 1D transient solid chains in the flowing layer, well separated from each other. In their approach, which describes granular flow down a rough inclined plane, a grain colliding with one of these chains transmits its momentum throughout the whole chain. On the other hand, for surface flows, Rajchenbach [6] suggests that each shock impact is imparted to the whole substrate beneath and thus that the momentum is "dissipated" into the whole packing. This description has been recently extended by Andreotti and Douady [10] who account for the possible trapping of flowing grains in the bumps of the static bed. The localization of the flow within a layer of finite thickness and the linear velocity profile are well reproduced in their approach. However, their model does not account for the invariance of the velocity gradient when the thickness R or the local angle u are changed. The collective objects that we observe are neither chainlike nor the whole packing, but rather widely distributed fractal clusters involving 50% of the beads of the flowing layer. No typical length 034301-3 034301-3 scale can be defined. This last point may be responsible for the selection of a velocity gradient independent of both R and u, i.e., independent of the shear stress. The observation of the successive clusterized frames shows that the largest clusters are emitted by the static phase and die either by fragmenting into smaller ones or by sticking back to the static phase. Differences between granular flows down a rough inclined plane and granular surface flows can then be rationalized as follows: clustering results from the competition between inelastic multiple collisions, which tends to aggregate grains together [START_REF] Goldhirsch | Physics of Dry Granular Media[END_REF], and shear, which erodes clusters. For flows down an inclined plane, these two effects lead to clusters of a typical size, whereas in surface flows, the static bed plays the role of a cluster reservoir: its erosion by the flowing grains can generate very large size clusters that then split and cascade into smaller and smaller ones (this mechanism might explain why we observe power-law distributions). In this view, the cluster size distribution, and a fortiori the velocity gradient, depend crucially on the boundary conditions.

The absence of a characteristic correlation length to describe the locally "jammed" clusters of beads is quite interesting. First, this contrasts with an assumption made in most local and nonlocal models [4,6,8,10], where the transition between the "solid" phase and the "liquid" phase is supposed to occur over a well defined length scale. The existence of multiscale rigid clusters indicates that the flowing phase is actually critical and suggests the proximity of a continuous "jamming" transition, of the type proposed in [17] (see also [18] for a related discussion). In this respect, it is interesting to note that similar power-law distributed clusters of strongly correlated motion have recently been observed in a colloid close to the glass transition [19].

FIG. 1 .

 1 FIG. 1. (a) Sketch of the experimental setup. (b) A typical instanteous velocity field obtained for quasi-2D flows (geometry A) with V 8 rpm in the rectangular box of (a). (c) Velocity profile measured at the center of the drum for three different rotation speeds: V 1 rpm, V 3 rpm, and V 8 rpm. The arrows indicate the end of each profile. The precision on each point is better than the point size.

FIG. 2 (

 2 FIG. 2 (color). Clusters of beads with correlated velocity fluctuation orientation in quasi-2D flows [geometry (A)]. (a) Typical instantaneous field of the orientation of the beads velocity fluctuations obtained for a quasi-2D flow [geometry (A)] with V 8 rpm. The original instantaneous velocity field is represented in gray. (b) The corresponding orientation correlation functions as a function of particle separation (see text for details). (c) Distribution function P͑d͒ of the angle d between the orientation of the velocity fluctuations of two beads in contact. The straight line corresponds to a uniform function. (d),(e) Typical frames of the isolated clusters, respectively, for V 8 rpm and V 1 rpm (see text for details). The position of the beads belonging to the flowing layer are represented with a black circle.The gray arrows represent the raw velocity field. Beads belonging to the same cluster are drawn in the same color (f ) Probability density function of the size (in number of beads N) of these clusters. Data have been shifted for clarity. From bottom to top, V 1 rpm, V 2 rpm, V 3 rpm, V 4 rpm, V 6 rpm, and V 8 rpm. For these different V, i.e., for different u and R, the distribution decreases as a power law (straight line) with the same exponent a ͑A͒ 2.9 6 0.1.

FIG. 3 (

 3 FIG. 3 (color). Presence of rigid clusters in the flowing layer of pure 2D flows [geometry (B)]. (a) Averaged volume fraction profile n m ͑z͒. The straight line corresponds to n n c . (b) Cluster size N plotted against their radius of gyration R g . The error bar indicates the standard deviation on R g of clusters with a given N. A power-law behavior (straight line) indicates fractal structures of dimension d ͑A͒ f 1.6 6 0.1. (c) Probability density function of the size (in number of beads N) of the "rigid" clusters (see text for details). Data have been shifted for clarity.From bottom to top, V 2 rpm, V 7 rpm, V 10 rpm, and V 13 rpm. They are power law distributed with an exponent close to a ͑B͒ n Ӎ 1.5 weakly dependent of V. Frames (d), (e), and (f) show a typical sequence of these clusters. Beads belonging to the static phase appear in gray. Beads belonging to the same cluster are drawn in the same color. The clusters are tracked on the successive frames. Clusters keep their color when they remain in successive frames.

We thank B. Dubrulle, O. Pouliquen, and D. Salin for fruitful discussions and C. Gasquet, V. Padilla, and P. Meininger for technical assistance.