Hugo A Martin

Anne Mangeney

Xiaoping Jia

Bertrand Antti Maury

Aline Lefebvre-Lepot

Yvon Maday

Paul Dérand

Ultrasound-induced dense granular flows: a two-time scale modeling

come   L'archive ouverte pluridisciplinaire

Introduction

The shear instability in a granular medium is involved in many natural hazards such as seismic fault slips and landslides. One of important and challenging issue in seismic hazard investigation is to understand how small-amplitude seismic waves (of the order of microstrain) can trigger remotely earthquakes from hundreds of km [START_REF] Gomberg | Earthquake triggering by seismic waves following the Landers and Hector Mine earthquakes[END_REF][START_REF] Hill | Seismicity Remotely Triggered by the Magnitude 7.3 Landers, California, Earthquake[END_REF] or more locally landslides around a volcano [START_REF] Keefer | Investigating landslides caused by earthquakes -A historical review[END_REF]. Recent observations showed that perturbations from local foreshock activity are probably a part of the earthquake nucleation process [START_REF] Bouchon | The long precursory phase of most large interplate earthquakes[END_REF] and that large rockfall events and avalanches can be triggered by volcanic seismicity [START_REF] Durand | On the Link Between External Forcings and Slope Instabilities in the Piton de la Fournaise Summit Crater, Reunion Island[END_REF]. Dynamic stress from seismic waves can destabilize granular solids and force failure earlier in time relative to an unperturbed fault or slope. Indeed, static and dynamic properties of dense granular media are determined by inhomogeneous contact force networks, exhibiting multiple metastable configurations. Sound waves propagating from grain to grain provide not only a unique probe of such optically opaque networks [START_REF] Liu | Sound in sand[END_REF][START_REF] Jia | Ultrasound propagation in externally stressed granular media[END_REF]) but also a controlled perturbation via vibration-induced softening and dissipation [START_REF] Johnson | Nonlinear dynamics, granular media and dynamic earthquake triggering[END_REF][START_REF] Jia | Elastic weakening of a dense granular pack by acoustic fluidization: Slipping, compaction, and aging[END_REF]. Granular media undergo a transition from a jammed solid state to a flowing liquid state when the external shear exceeds the static yield stress (Fig. 1(a,b)).

Previous works pointed this transition as bifurcation phenomenon [START_REF] Jaeger | Friction in Granular Flows[END_REF][START_REF] Quartier | Dynamics of a grain on a sandpile model[END_REF][START_REF] Baldassarri | Brownian Forces in Sheared Granular Matter[END_REF][START_REF] Dijksman | Jamming, Yielding, and Rheology of Weakly Vibrated Granular Media[END_REF], similar to solid friction at multicontact interfaces [START_REF] Baumberger | Solid friction from stick-slip down to pinning and aging[END_REF] described by the rate and state friction law [START_REF] Marone | LABORATORY-DERIVED FRICTION LAWS AND THEIR APPLICATION TO SEISMIC FAULTING[END_REF][START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF]. Here, the friction coefficient is defined as 𝜇 = 𝜏/𝜎 𝑛 , as the ratio of the shear stress normalized by the normal stress from which the static and dynamic coefficients of friction 𝜇 𝑠,𝑑 = 𝜏 𝑠,𝑑 /𝜎 𝑛 follow, where 𝜏 𝑠 is the static friction stress at yield while 𝜏 𝑑 is the dynamic friction stress. In the inclined plane geometry considered in this article (Fig. 1(a)), we have 𝜇 = tan 𝜃 and 𝜇 𝑠 = tan 𝜃 𝑚 , with 𝜃 𝑚 the (maximum) angle of avalanche. The angle of repose 𝜃 𝑟 being a few percent lower than 𝜃 𝑚 [START_REF] Pouliquen | Onset of Granular Flows on an Inclined Rough Surface: Dilatancy Effects[END_REF][START_REF] Daerr | Two types of avalanche behaviour in granular media[END_REF][START_REF] Coussot | Avalanche Behavior in Yield Stress Fluids[END_REF][START_REF] Wyart | On the dependence of the avalanche angle on the granular layer thickness[END_REF]) corresponds to the dynamic friction 𝜇 𝑑 = tan 𝜃 𝑟 at the minimum shear load (Fig. 1(b)). It has been shown [START_REF] Baumberger | Solid friction from stick-slip down to pinning and aging[END_REF][START_REF] Nasuno | Friction in Granular Layers: Hysteresis and Precursors[END_REF][START_REF] Bureau | Jamming creep of a frictional interface[END_REF][START_REF] Parteli | Nontrivial temporal scaling in a Galilean stick-slip dynamics[END_REF]) that for shear forces far below the static threshold 𝜇 ≪ 𝜇 𝑠 , both granular layers and rough solid interfaces respond elastically as shown in Fig. 1(b), in the jammed state (region I). For 𝜇 ≲ 𝜇 𝑠 a shear lower than the threshold, nonlinear response occurs with creep-like irreversible motion. For 𝜇 ≳ 𝜇 𝑠 , the system yields and starts to slide or flow over a transient characteristic distance [START_REF] Marone | LABORATORY-DERIVED FRICTION LAWS AND THEIR APPLICATION TO SEISMIC FAULTING[END_REF][START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF][START_REF] Baumberger | Solid friction from stick-slip down to pinning and aging[END_REF]) before reaching the steady flow region II/III at a velocity 𝑉 or shear rate 𝛾 imposed by the load.

The possible failure of a granular medium, such as a sand pile, caused by external vibrations, has been known for a long time in engineering and geophysical applications, however a unified physical description still lacks. The vibrations considered are most of the time of large amplitude 𝑈 0 ≳ 𝑑 with 𝑑 the grain size and low frequency 𝑓 < 𝑓 0 , where 𝑓 0 is a characteristic frequency determined by the stiffness of interfacial contacts [START_REF] Bureau | Jamming creep of a frictional interface[END_REF]. The amount of shaking is usually estimated by the reduced peak acceleration of the grain Γ = 𝑎/𝑔 with 𝑎 the instantaneous acceleration and 𝑔 the gravity. When Γ > 1, vertical vibrations cancel almost normal forces exerting on the grain (confined under gravity) and modify consequently the spatial arrangement of grains, resulting in phenomena such as compaction, convection, shear banding, to mention a few [START_REF] Jaeger | Relaxation at the Angle of Repose[END_REF][START_REF] Clement | Fluidization of a Bidimensional Powder[END_REF][START_REF] D'anna | Observing brownian motion in vibrationfluidized granular matter[END_REF]. This is similar to the oscillation effect on the normal stress facilitating sliding [START_REF] Zaloj | Modifying Friction by Manipulating Normal Response to Lateral Motion[END_REF][START_REF] Cochard | Stabilization of Frictional Sliding by Normal Load Modulation[END_REF] and also to the scenario of the acoustic lubrication in a confined continuous medium. In this scenario, the acoustic pressure 𝑝 𝑎 = (𝜌𝑐)𝑣 𝑎 , with 𝑐 the sound speed and 𝑣 𝑎 the vibration velocity, is expected to The inertial flow triggered by ultrasonic vibration is mostly uniform and continuous.

(b) Sketch of the normalized stress-strain rate relation in a sheared granular medium: static state (V = 0), unstable flow (velocity-weakening) and stable flow (velocity strengthening). Under a shear 𝜇 between 𝜇 𝑠 and 𝜇 𝑑 , the metastable state can be switched abruptly to a flowing state by acoustic perturbation (lubrication).

temporally relieve the pressure of the overburden, thereby decreasing the yield stress [START_REF] Melosh | Dynamical weakening of faults by acoustic fluidization[END_REF]. However, the above scenarios, involving large-amplitude vibrations, cannot explain the dynamic earthquake triggering by seismic waves at micro-and nano-strain amplitude [START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF][START_REF] Gomberg | Earthquake triggering by seismic waves following the Landers and Hector Mine earthquakes[END_REF], nor the laboratory experiments using nanometer-amplitude ultrasound to soften the material modulus by 30 % via nonlinear dynamics [START_REF] Johnson | Nonlinear dynamics, granular media and dynamic earthquake triggering[END_REF][START_REF] Jia | Elastic weakening of a dense granular pack by acoustic fluidization: Slipping, compaction, and aging[END_REF]. Also, some modifications of the stick-slip cycle by ultrasound remain unexplained [START_REF] Johnson | Effects of acoustic waves on stick-slip in granular media and implications for earthquakes[END_REF]. In these situations, the oscillation frequency of ultrasound 𝑓 ⩾ 40 kHz is high compared to the characteristic frequency of 𝑓 0 ≃ 5 kHz in millimeter-thick granular layers [START_REF] Baumberger | Solid friction from stick-slip down to pinning and aging[END_REF][START_REF] Bureau | Jamming creep of a frictional interface[END_REF]) so that grains cannot have enough inertial normal motion to suppress the weight of the overburden. On the other hand, for a nanometer ultrasonic vibration, the collision-like pressure estimated as 𝑝 𝑐 ≃ 10 -4 Pa ≪ 𝜎 𝑛 (≃ 10 Pa) is too small to be considered.

Recently, it was evidenced by lab experiments (Leópoldès et al. 2020) that the triggering of granular instabilities by small-amplitude (𝑈 0 /𝑑 ≃ 10 -5 ) and high-frequency ( 𝑓 = 70 kHz) ultrasound waves can be explained by acoustic shear lubrication of grain contacts. This interparticle friction reduction consequently lowers the effective friction coefficient on the macroscopic scale and triggers the granular flow at an inclined angle below 𝜃 𝑚 (Fig. 1(b)) without causing the rearrangement of grain positions, . However, it remains unclear how such local effects could give rise to collective motion and also their possible delay of response [START_REF] Durand | On the Link Between External Forcings and Slope Instabilities in the Piton de la Fournaise Summit Crater, Reunion Island[END_REF]. In particular, the micro processes of triggering remain poorly understood, as small-amplitude ultrasound does not induce grain displacement per se at the relevant length scale (i.e., grain diameter 𝑑) during avalanches.

In this article, we address these questions by presenting numerical simulations of the onset of flows of a granular layer initially static on an inclined plane, triggered by ultrasound applied from the basal plane. To do so, we propose an original coupling model, relying on a timescale separation, with one time scale representing the grain-scale dynamics (called the grainmotion time scale) and the other representing the ultrasonic vibration in the granular medium (called the vibration time scale). At the vibration time scale, the granular packing (particles) is considered quasi-static. Therefore, our methodology involves coupling a Discrete Element Model (DEM) for grain flows with a mass-spring model (network) for wave vibration. The latter is investigated in the steady (harmonic) regime through vibrational (eigen) modes. Once the amplitudes of the vibration field are computed, the acoustic lubrication effect on the (entire) granular medium is taken into account through modifications of the interparticle friction coefficients, according to the Mindlin model (Leópoldès et al. 2020). To our knowledge, this is the first time such a coupling numerical methodology has been employed to study the acoustic lubrication effect on the granular flow.

The advantages of this time-scale separation are twofold. Firstly, it simplifies (isolates) the simulation of acoustic propagation through a "frozen" (or a snapshot of) granular network at a given flowing instance. Here the interparticle friction coefficients may be modified solely by the irreversible ultrasound-matter interaction [START_REF] Johnson | Nonlinear dynamics, granular media and dynamic earthquake triggering[END_REF][START_REF] Jia | Elastic weakening of a dense granular pack by acoustic fluidization: Slipping, compaction, and aging[END_REF]Leópoldès et al. 2020) and are not affected by particle collisions or other phenomena occurring at the grain-motion time scale. Additionally, this approach offers high computational efficiency since we can use relatively large time steps (Δ𝑡 𝑔 = 1 ms at the grain-motion time scale) without losing the ultrasound-induced effects, thereby allowing for the consideration of large assemblies of grains that can be compared to experimental data in laboratory settings (see typical times of simulations in Martin et al. (2023a)). The computational time efficiency is discussed in Section 2.4. The model utilized at the grain motion time scale is the convex optimization contact dynamics (cocd) discrete element method, which has been proposed and validated in Martin et al. (2023a), and compared with other granular models in Martin et al. (2023b).

The numerical method is presented at various time scales in Section 2. Section 3 contains all the results we obtain, both with and without destabilization due to basal vibrations. In Section 4, we then discuss the physical interparticle mechanisms responsible for the destabilization we have identified, the changes that occur when we alter the vibration parameters, and compare our results with experiments.

Numerical models for granular flows and ultrasonic vibrations

To model the acoustic triggering of granular flows induced by basal ultrasounds at the laboratory scale, our strategy is to take into account the vibration-induced change of the interparticle friction coefficient 𝜇 𝑝 . This is done by considering the very different time scales of grain motion and ultrasounds, respectively. We then consider a grain-motion time scale that is of the order of small grain motions, and a vibration time scale.

Concerning the grain-motion time scale, for a granular flow on a slope, the horizontal velocity scale is given by 𝑈 =

√

𝑔𝐻 cos 𝜃, where 𝑔 is the gravitational constant, 𝐻 is the thickness of the granular mass, and 𝜃 is the slope angle. In our configuration, 𝐻 lies between 3𝑑 and 14.4𝑑, where 𝑑 is the mean diameter of the grains (0.7 mm), resulting in 0.002 ⩽ 𝐻 ⩽ 0.011 m, and thus 0.14 ≲ 𝑈 ≲ 0.32 m s -1 . We use a time scale where a grain moving at velocity 𝑈 covers its own radius, i.e., 𝐿 = 0.35 mm. This time scale is then given by 𝑇 = 𝐿/𝑈, and it falls within the range of approximately 10 -3 ≲ 𝑇 ≲ 2.5 × 10 -3 s. We adopt the smallest value as the common time scale for all our grain displacement simulations, i.e., the grain-motion time scale is Δ𝑡 𝑔 = 1 ms.

Regarding the vibration time scale, experimental measurements of the speed of sound in a granular assembly at rest and maintained by a low confined stress (for instance by gravity) have been made where small but finite values of sound speed 𝑐 = 10 m s -1 to 100 m s -1 were measured [START_REF] Liu | Sound in sand[END_REF][START_REF] Bonneau | Evidence of rayleigh-hertz surface waves and shear stiffness anomaly in granular media[END_REF][START_REF] Van Den Wildenberg | Evolution of granular packings by nonlinear acoustic waves[END_REF][START_REF] Brum | Drastic slowdown of the rayleigh-like wave in unjammed granular suspensions[END_REF]). These values are much larger than the sound speed predicted by the effective medium theory based on the simplified normal contact force (i.e., the Hertz model), 𝑐 ≃ 𝑝 1/6 in a granular pack loaded by a hydrostatic pressure 𝑝 = 𝜌𝑔ℎ. This would give rise to 𝑐 ≃ 0.1 m s -1 for ℎ = 5 mm to 10 mm. Such a discrepancy may stem from the interlocking effect (i.e., arching) due to frictional contact (tangential) forces and heterogeneous anisotropic stress networks (i.e., force-chains), which depend on the sample configuration and loading history (memory effect) as well as the confining boundary [START_REF] Jaeger | Friction in Granular Flows[END_REF]). As shown in [START_REF] Khidas | Anisotropic nonlinear elasticity in a spherical-bead pack: Influence of the fabric anisotropy[END_REF] and [START_REF] Johnson | Linear and Nonlinear Elasticity of Granular Media: Stress-Induced Anisotropy of a Random Sphere Pack[END_REF], elastic wave velocities (longitudinal and transversal) do remain finite in a compacted granular packing after the removal of the applied stress, likely due to the tight wedging or interlocking of grains which results in a residual stress network. Using again the mean grain radius as a characteristic distance, we get a time scale approximately from 𝑇 = 𝐿/𝑐, i.e., 𝑇 ≃ 12 µs, estimated with 𝑐 = 30 m s -1 (Leópoldès et al. 2020). Hence, we adopt this value as the time scale for sound propagation Δ𝑡 𝑤 = 10 µs. Note that it is smaller than the grain displacement time scale by a factor of 100, justifying our approach that separates these two time scales. Consequently, at the vibration time scale Δ𝑡 𝑤 , the granular assembly can be considered as frozen. We thus consider two different models at these two different time scales. At the grainmotion time scale Δ𝑡 𝑔 , grain motion is described based on the discrete element method cocd (Martin et al. 2023a), which represents the macroscopic motion of each particle in a granular assembly. Then, at the vibration time scale, the vibration model computes the infinitesimal perturbations of each grain position, around an equilibrium configuration provided by the grain-motion model. We first briefly introduce the grain-motion model cocd developed by Martin et al. (2023a) (Section 2.1). Then we present the the wave equation and the vibration model (Section 2.2), and finally describe the Mindlin model (Section 2.3) that deals with the vibration-induced reduction of interparticle friction coefficients. We close this modeling section by briefly presenting the computational time efficiency (Section 2.4).

Grain-motion model cocd

In our grain-motion model cocd, the granular media is represented by a collection of rigid particles (see Fig. 1(a)) like glass spheres (Martin et al. 2023a). The equations of motion are solved for each particle at every time step to determine their respective contact forces. These interactions are traditionally described using the Hertz theory, utilizing non-linear damped springs as commonly seen in the molecular dynamics (MD) framework (see Cundall's pioneering work in [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]). However, in this article, we adopt an alternative method called contact dynamics (CD), originally introduced by Moreau and Jean in the 1990s [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics : application to granular materials[END_REF][START_REF] Moreau | Numerical aspects of the sweeping process[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Moreau | An introduction to Unilateral Dynamics[END_REF]). In contrast to MD, where contact forces are modeled using functions derived from the Hertz theory, CD employs linear impulses. These impulses adhere to contact laws governing both normal repulsion and tangential friction.

Numerous numerical techniques have been put forth for CD models [START_REF] Maury | A time-stepping scheme for inelastic collisions: Numerical handling of the nonoverlapping constraint[END_REF][START_REF] Staron | Study of the collapse of granular columns using two-dimensional discretegrain simulation[END_REF][START_REF] Anitescu | Optimization-based simulation of nonsmooth rigid multibody dynamics[END_REF][START_REF] Tasora | Large-scale parallel multi-body dynamics with frictional contact on the graphical processing unit[END_REF][START_REF] Radjai | Contact dynamics as a nonsmooth discrete element method[END_REF][START_REF] Acary | A formulation of the linear discrete Coulomb friction problem via convex optimization[END_REF][START_REF] Seguin | Clustering and flow around a sphere moving into a grain cloud[END_REF]Martin et al. 2023a;[START_REF] Bloch | On convex numerical schemes for inelastic contacts with friction[END_REF]. Close to the framework of the SCoPI software (2022), in cocd's approach, particle velocities and positions are computed simultaneously through an implicit scheme, necessitating the solution of a convex minimization problem during each time integration step.

Within the CD framework, two contact laws are validated at each contact and time step. The first law establishes a complementarity problem between the normal distance and the intensity of the normal contact force. This implies that grains cannot overlap or engage in interaction unless they are in contact, and the force between any two grains is inherently repulsive. In mathematical terms, this relationship is represented as

𝑓 𝑖 𝑗 𝑛 ⩾ 0 , D 𝑖 𝑗 ⩾ 0 , 𝑓 𝑖 𝑗 𝑛 D 𝑖 𝑗 = 0,
where 𝑖 and 𝑗 denote two particles, D 𝑖 𝑗 is their normal distance, and 𝑓 𝑖 𝑗 𝑛 represents the intensity of the normal force between them (see Fig. 2(a) and Fig. 2(b)). The second validated contact law pertains to Coulomb friction law (Fig. 2(a) and Fig. 2(c)), encompassing both the tangential and normal components of the contact force belonging to Coulomb's cone, is expressed as

𝒇 𝑖 𝑗 𝑡 = -𝜇 𝑖 𝑗 𝑓 𝑖 𝑗 𝑛 𝒗 𝑖 𝑗 𝑡 /∥𝒗 𝑖 𝑗 𝑡 ∥, if ∥𝒗 𝑖 𝑗 𝑡 ∥ > 0, ∥ 𝒇 𝑖 𝑗 𝑡 ∥ ⩽ 𝜇 𝑖 𝑗 𝑓 𝑖 𝑗 𝑛 , if ∥𝒗 𝑖 𝑗
𝑡 ∥ = 0, where ∥ • ∥ denotes the euclidean norm, 𝜇 𝑖 𝑗 is the interparticle friction coefficient (𝜇 𝑝 ), 𝒇 𝑖 𝑗 𝑡 ∈ R 3 denotes the tangential force vector, and 𝒗 𝑖 𝑗 𝑡 ∈ R 3 is the tangential relative velocity vector between the two spheres 𝑖 and 𝑗. To be as general as possible, we will subsequently introduce the equations of motion in 3D despite the simulations that are presented in this article are in 2D. Consider a mechanical system in R 3 , consisting of 𝑁 rigid spheres capable of rotation, each with specified fixed radii 𝑟 𝑖 > 0 and masses 𝑚 𝑖 > 0, 𝑖 = 1, . . . , 𝑁. The center of sphere 𝑖 is represented by 𝒄 𝑖 ∈ R 3 , and its instantaneous velocity by 𝒗 𝑖 ∈ R 3 . As we are exclusively dealing with spheres, we do not track the orientation of bodies; instead, we focus solely on the instantaneous rotation vector, denoted as 𝝎 𝑖 ∈ R 3 . We denote by 𝒄 = (𝒄 1 , . . . , 𝒄 𝑁 ) ∈ R 3𝑁 and 𝒖 = (𝒗 1 , 𝝎 1 , . . . , 𝒗 𝑁 , 𝝎 𝑁 ) ∈ R 6𝑁 the generalized position and velocity field vectors.

The signed distance between spheres 𝑖 and 𝑗 is defined by :

D 𝑖 𝑗 (𝒄) = ∥𝒄 𝑖 -𝒄 𝑗 ∥ -𝑟 𝑖 + 𝑟 𝑗 , (2.1)
so that the non-overlapping condition writes D 𝑖 𝑗 ⩾ 0.

For any pair of grains 𝑖 and 𝑗, with respective centers represented by 𝒄 𝑖 and 𝒄 𝑗 , we use 𝑪 𝑖 and 𝑪 𝑗 to indicate the points that establish the distance (with 𝑪 𝑖 = 𝑪 𝑗 when the spheres are in contact; refer to Fig. 2(d)). We define the corresponding position vectors as 𝒓 𝑖 = 𝑪 𝑖 -𝒄 𝑖 , 𝒓 𝑗 = 𝑪 𝑗 -𝒄 𝑗 .

We define the direction perpendicular to the surfaces of the particles at points 𝑪 𝑖 and 𝑪 𝑗 , a direction common to both particles. We introduce the unit vector 𝒏 𝑖 𝑗 ∈ R 3 , which is defined as the corresponding normal vector pointing towards particle 𝑖. Given that we are dealing with spherical particles, we have:

𝒏 𝑖 𝑗 = 𝒄 𝑖 -𝒄 𝑗 ∥𝒄 𝑖 -𝒄 𝑗 ∥ .
We denote by P 𝑖 𝑗 𝒗 = 𝒗 -(𝒗 • 𝒏 𝑖 𝑗 )𝒏 𝑖 𝑗 ∈ R 3 the projection of 𝒗 on Π 𝑖 𝑗 , the plane that is orthogonal to 𝒏 𝑖 𝑗 and thus parallel to the tangent planes in 𝑪 𝑖 and 𝑪 𝑗 .

We also define A 𝑖 𝑗 from R 6𝑁 to R 3 as the linear operator which maps the generalized velocity field 𝒖 ∈ R 6𝑁 to the relative velocity between the points 𝑪 𝑖 and 𝑪 𝑗 at which the distance between spheres 𝑖 and 𝑗 is attained, i.e.

A 𝑖 𝑗 𝒖 = 𝒗 𝑖 + 𝝎 𝑖 ∧ 𝒓 𝑖 -(𝒗 𝑗 + 𝝎 𝑗 ∧ 𝒓 𝑗 ) ∈ R 3 .
Direct calculations demonstrate that for any generalized velocity 𝒖 ∈ R 6𝑁 and any vector

𝒇 ∈ R 3 , we obtain A 𝑖 𝑗 𝒖 • 𝒇 = 𝒖 • A 𝑇 𝑖 𝑗 𝒇 with A 𝑇 𝑖 𝑗 𝒇 = (0, . . . , 0, 𝒇 , 𝒓 𝑖 ∧ 𝒇 position 𝑖 , 0, . . . , 0, -𝒇 , -𝒓 𝑗 ∧ 𝒇 position 𝑗 , 0, . . . , 0) ∈ R 6𝑁 ,
so that A 𝑇 𝑖 𝑗 maps a vector 𝒇 ∈ R 3 to the generalized force/moment vector corresponding to the force 𝒇 exerted on particle 𝑖 at point 𝑪 𝑖 and the opposite force -𝒇 exerted on particle 𝑗 at point 𝑪 𝑗 .

The vector P 𝑖 𝑗 A 𝑖 𝑗 𝒖 = 𝒗 𝑖 𝑗 𝑡 ∈ R 3 represents the tangential relative velocity. Consequently, when two spheres are in contact without any relative normal motion (i.e., 𝒏 𝑖 𝑗 • A 𝑖 𝑗 𝒖 = 0), the expression ∥P 𝑖 𝑗 A 𝑖 𝑗 𝒖∥ = 0 indicates a rolling motion with no slip, while ∥P 𝑖 𝑗 A 𝑖 𝑗 𝒖∥ > 0 corresponds to a sliding motion.

At any time, we shall denote by 𝐼 𝑐 the set of all possible pairs of contacts: 𝐼 𝑐 = {(𝑖, 𝑗) 1 ⩽ 𝑖 < 𝑗 ⩽ 𝑁 }. Note that the pair of grains 𝑖 and 𝑗 is represented only once in 𝐼 𝑐 through the couple (𝑖, 𝑗) if 𝑖 < 𝑗 and ( 𝑗, 𝑖) if 𝑗 < 𝑖. We denote by 𝑁 𝑐 the total number of potential pairs of spheres, which is also the cardinal of set 𝐼 𝑐 , i.e., 𝑁 𝑐 = card(𝐼 𝑐 ) = 𝑁 (𝑁 -1)/2.

We consider that no external torque is exerted on the grains. If 𝒇 𝑒𝑥𝑡 𝑖 ∈ R 3 is the external force exerted on particle 𝑖, we define the generalized force vector as

𝒇 𝑒𝑥𝑡 = ( 𝒇 𝑒𝑥𝑡 1 , 0, . . . , 𝒇 𝑒𝑥𝑡 𝑛 , 0) ∈ R 6𝑁 .
We then define the 6𝑁 × 6𝑁 generalized mass matrix (masses and moments of inertia) as

M = diag (𝑚 1 , 𝑚 1 , 𝑚 1 , 𝐽 1 , 𝐽 1 , 𝐽 1 , 𝑚 2 , . . . , 𝐽 𝑁 , 𝐽 𝑁 , 𝐽 𝑁 ) .
Finally, the equations of motion write:

M d𝒖 d𝑡 = 𝒇 𝑒𝑥𝑡 + ∑ 𝛼∈ 𝐼 𝑐 A 𝑇 𝛼 𝑓 𝛼 𝑛 𝒏 𝛼 + 𝒇 𝛼 𝑡 , (2.2) 𝑓 𝛼 𝑛 ⩾ 0 , D 𝛼 ⩾ 0 , 𝑓 𝛼 𝑛 D 𝛼 = 0, 𝛼 ∈ 𝐼 𝑐 , (2.3) If D 𝛼 (𝒄) = 0 then (A 𝛼 𝒖 + ) • 𝒏 𝛼 = 0, 𝛼 ∈ 𝐼 𝑐 , (2.4) If ∥P 𝛼 A 𝛼 𝒖 + ∥ > 0 (sliding motion) , 𝒇 𝛼 𝑡 = -𝜇 𝛼 𝑓 𝛼 𝑛 P 𝛼 A 𝛼 𝒖 + ∥P 𝛼 A 𝛼 𝒖 + ∥ , 𝛼 ∈ 𝐼 𝑐 , (2.5) If ∥P 𝛼 A 𝛼 𝒖 + ∥ = 0 (no slip) , ∥ 𝒇 𝛼 𝑡 ∥ ⩽ 𝜇 𝛼 𝑓 𝛼 𝑛 , 𝛼 ∈ 𝐼 𝑐 . (2.6) Equation (2.4
) is added to the normal (Eq. ( 2.3)) and tangential contact laws (Eqs. (2.5)-

(2.6)), specifying an inelastic collision law. Observe that translational and rotational velocities are prone to being non-smooth, experiencing instantaneous jumps during collisions. Specifically, the post-collision velocity 𝒖 + may differ from the pre-collision velocity 𝒖 -. Consequently, the mentioned evolution is to be interpreted in a weak, distributional sense. Let's provide some additional remarks on the preceding equations. For a pair of grains 𝛼 = (𝑖, 𝑗) ∈ 𝐼 𝑐 , where 𝐼 𝑐 denotes the set of contacts, the corresponding vector 𝑓 𝑖 𝑗 𝑛 𝒏 𝑖 𝑗 + 𝒇 𝑖 𝑗 𝑡 ∈ R 3 is transmitted to both particles 𝑖 and 𝑗 through the transpose of A 𝑖 𝑗 . To elaborate, let's introduce the following definitions:

𝑓 𝑗𝑖 𝑛 = 𝑓 𝑖 𝑗 𝑛 , 𝒇 𝑗𝑖 𝑡 = -𝒇 𝑖 𝑗 𝑡 , ∀𝛼 = (𝑖, 𝑗) ∈ 𝐼 𝑐 .
Then, utilizing the expression for A 𝑇 𝑖 𝑗 , Eq. ( 2.2) can be reformulated as follows:

𝑚 𝑖 𝒗 𝑖 = 𝒇 𝑒𝑥𝑡 𝑖 + ∑ 𝑗, 𝑗≠𝑖 ( 𝑓 𝑖 𝑗 𝑛 𝒏 𝑖 𝑗 + 𝒇 𝑖 𝑗 𝑡 ), ∀𝑖 = 1 . . . 𝑁, 𝐽 𝑖 𝝎 𝑖 = ∑ 𝑗, 𝑗≠𝑖 (𝒓 𝑖 ∧ 𝒇 𝑖 𝑗 𝑡 ), ∀𝑖 = 1 . . . 𝑁.
This corresponds to Newton's second law, where the contact between particles 𝑖 and 𝑗 induces the force 𝑓 From Eq. ( 2.3), we deduce 𝑓 𝑖 𝑗 𝑛 = 𝑓 𝛼 𝑛 ⩾ 0. This, combined with the orientation of 𝒏 𝑖 𝑗 from particle 𝑗 to particle 𝑖, ensures that this force is repulsive, as anticipated. Equation (2.3) also guarantees that the distances between the particles remain positive, and the normal force is null whenever the distance is strictly positive (i.e., the particles are not in contact).

The mechanical characteristics of such granular media arise from a synergy of geometrical particle rearrangements and interparticle friction forces. To be precise, the macroscopic static friction coefficient 𝜇 𝑠 = tan(𝜃 𝑚 ) -with 𝜃 𝑚 the (maximum) angle of avalanche -can be understood as a composite of the interparticle friction coefficient 𝜇 𝑝 and the geometric confinement effect (dilatancy) 𝜇 𝑔 (Leópoldès et al. 2020). The coefficient 𝜇 𝑔 is influenced by factors such as grain shapes, masses, or inertia, while 𝜇 𝑝 serves as a defined parameter within the model, integral to the classical Coulomb law of friction for all grain-to-grain or grain-to-wall interactions.

The value of 𝜇 𝑝 employed in this article is 𝜇 𝑝 = 0.25. This choice falls within a comparable range to the friction coefficient 𝜇 𝑝 = 0.3, calibrated through 3D simulations and experiments (Martin et al. 2023a,b), and the friction coefficient measured for an ideal glass-to-glass contact, 𝜇 𝑝 = 0.4 (www.engineeringtoolbox.com 2022), as well as that determined using DEM, 𝜇 𝑝 = 0.16 [START_REF] Tang | Measurement of restitution and friction coefficients for granular particles and discrete element simulation for the tests of glass beads[END_REF].

Wave equation and vibrational modes 2.2.1. Wave equation

At the wave propagation or vibration time scale, the grains motions may be supposed quasi-static or frozen. As mentioned above, however, the quantitative description of sound propagation in such weakly confined amorphous-like granular media is not available [START_REF] Makse | Granular packings: Nonlinear elasticity, sound propagation, and collective relaxation dynamics[END_REF]). To capture qualitatively the interaction between and ultrasound granular flow (see Section 2.3), we model the vibration of grains, for a first approximation, as in a 2D network of mass-spring (mimicking a normal contact stiffness). Here the tangential force and rotational motion are neglected , i.e., the particles are considered as frictionless. Nevertheless, the friction forces are accounted for in the grain-motion model cocd and for investigating the lubrication effect at grain contacts induced by shear acoustic waves (see below section 2.3). The sound propagation is characterized by the perturbed positions of the centers of the masses 𝒄 𝑖 ∈ R 3 ; see .e.g., a similar model proposed in [START_REF] Somfai | Elastic wave propagation in confined granular systems[END_REF]. Thanks to an expansion around an assumed equilibrium configuration, we establish a wave equation for infinitesimal perturbations from this equilibrium position. The full description of the equations derivation can be found in Appendix A.

Similarly to the operator P 𝑖 𝑗 A 𝑖 𝑗 and since we do not account for rotational motion, we define N 𝑖 𝑗 as the linear operator from R 3𝑁 to R which maps the generalized grain velocity vector of translation ū = (𝒗 1 , 𝒗 2 , . . . , 𝒗 𝑁 ) ∈ R 3𝑁 to the relative normal velocity between the spheres 𝑖 and 𝑗 (grains), projected on the line generated by the normal vector 𝒏 𝑖 𝑗 , i.e.,

N 𝑖 𝑗 ū = (𝒗 𝑖 -𝒗 𝑗 ) • 𝒏 𝑖 𝑗 ∈ R.
Straightforward computations show that for any generalized velocity ū ∈ R 3𝑁 and any scalar

𝑓 𝑛 ∈ R, we have 𝑓 𝑛 N 𝑖 𝑗 ū = ū • N 𝑇 𝑖 𝑗 𝑓 𝑛 with N 𝑇 𝑖 𝑗 𝑓 𝑛 = (0, . . . , 0, 𝑓 𝑛 𝒏 𝑖 𝑗 position 𝑖 , 0, . . . , 0, -𝑓 𝑛 𝒏 𝑖 𝑗 position 𝑗 , 0, . . . , 0) ∈ R 3𝑁 , (2.7)
so that N 𝑇 𝑖 𝑗 maps a scalar 𝑓 𝑛 ∈ R to the generalized force vector corresponding to the force 𝑓 𝑛 𝒏 𝑖 𝑗 exerted on particle 𝑖 at point 𝒄 𝑖 and the opposite force -𝑓 𝑛 𝒏 𝑖 𝑗 exerted on particle 𝑗 at point 𝒄 𝑗 .

We then define the linear operator N from R 3𝑁 into R 𝑁 𝑐 corresponding to the combination of all maps N 𝛼 , 𝛼 ∈ 𝐼 𝑐 , i.e., for any 𝒗 ∈ R 3𝑁 , we have N 𝒗 = (N 1,2 𝒗, N 1,3 𝒗, . . . , N 𝑁 -1, 𝑁 𝒗) ∈ R 𝑁 𝑐 , and for any 𝒇 ∈ R 𝑁 𝑐 we have the equality N 𝒗 • 𝒇 = 𝒗 • N 𝑇 𝒇 .

We now define the 𝑁 𝑐 × 𝑁 𝑐 diagonal square matrix K , which contains the elastic properties of the system by

K = 3 2 diag (𝜅 1,2 ) 2/3 ( 𝑓 1,2 𝑛 ) 1/3 , . . . . . . , (𝜅 𝑖 𝑗 ) 2/3 ( 𝑓 𝑖 𝑗 𝑛 ) 1/3 , . . . . . . , (𝜅 𝑁 -1, 𝑁 ) 2/3 ( 𝑓 𝑁 -1, 𝑁 𝑛 ) 1/3 ∈ R 𝑁 𝑐 × 𝑁 𝑐 ,
where 𝑓 𝑖 𝑗 𝑛 still represents the intensity of the normal force between particles 𝑖 and 𝑗 and where 𝜅 𝑖 𝑗 > 0 is a constant depending on grains properties; see Appendix A.1.

We finally define the 3𝑁 × 3𝑁 generalized mass matrix (masses only) as

M = diag (𝑚 1 , 𝑚 1 , 𝑚 1 , 𝑚 2 , 𝑚 2 , 𝑚 2 . . . , 𝑚 𝑁 , 𝑚 𝑁 , 𝑚 𝑁 ) . (2.8)
At the end, a wave equation is defined for any 𝒆 ∈ R 3𝑁 by

M d 2 𝒆 d𝑡 2 + 𝚲𝒆 = 0, (2.9)
where the linear map defined by the matrix

𝚲 = N 𝑇 KN ∈ R 3𝑁 ×3𝑁 ,
can be seen as a kind of discrete Laplace operator, and where 𝜀 indicates that 𝒆 is only an infinitesimal perturbation of configuration 𝒄. The operator 𝚲 depends on the contact network through operator N, formed by the generalized position vector at equilibrium 𝒄 and embeds the elastic properties of the granular assembly, characterized by the constant normal force intensities exerted between any particles 𝑖 and 𝑗, i.e., 𝑓 In the classical form of the wave equation (𝜕 2 𝑡 𝑢 -𝑐 2 Δ𝑢 = 0), the square of the sound speed (constant) 𝑐 2 is before the Laplace operator Δ. In our framework Eq. (2.9), the (local) sound speeds are then proportional to ((𝜅 𝑖 𝑗 ) 2/3 ( 𝑓 𝑖 𝑗 𝑛 ) 1/3 /𝑚 𝑖 ) 1/2 so to the 1/6th power of the normal forces computed at the grain-motion time scale 𝑓 𝑖 𝑗 𝑛 , i.e., 𝑐 ∝ ( 𝑓 𝑖 𝑗 𝑛 ) 1/6 (see discussions above and in Appendix A.1) .

Vibration model

We now present the vibration model which describes the asymptotic limit of the ultrasound vibrations (perturbations) induced in a quasi-static (frozen) granular packing. The generalized position vector of a grain at equilibrium, denoted 𝒄 ∈ R 3𝑁 is excited at a mono-frequency by an external vibration applied from the basal plane. More precisely, the grains in contact with the bottom are submitted to a sinusoidal motion, at a given frequency 𝑓 = 𝜔/2𝜋 and amplitude 𝑈 0 .

At this vibration time scale, we seek for the vibrational (eigen) modes of the granular assembly, on the line of [START_REF] Leibig | Model for the propagation of sound in granular materials[END_REF] and [START_REF] Somfai | Elastic wave propagation in confined granular systems[END_REF], and consider the asymptotic limit regime (i.e., steady harmonic vibrations) of the perturbed positions of all grains at the forced external vibration frequency 𝑓 . Accordingly, the solution to the wave equation (2.9) 𝒆 ∈ R 3𝑁 is separable, and can be written as 𝒆(𝑡) = 𝑈 0 exp(𝑖𝜔𝑡)𝒒, where 𝒒 ∈ R 3𝑁 is a constant vector that does not depend on time, solution to the vibration model (or Helmoltz equation [START_REF] Somfai | Elastic wave propagation in confined granular systems[END_REF][START_REF] Couto | Green's functions for the wave, Helmholtz and Poisson equations in a two-dimensional boundless domain[END_REF])). This equation writes

𝜔 2 1 3𝑁 + 𝚲 𝒒 = 0 ∈ R 3𝑁 ,
(2.10) associated to the Dirichlet kind of boundary condition

𝑞 𝑖 𝑦 = 1, ∀𝑖 ∈ J , (2.11) (a) (b)
Figure 4: Illustration of the vibration-induced perturbations of the interparticle friction coefficients.

where J ⊂ N is the set of grain indexes in contact with the bottom, 𝑞 𝑖 𝑦 ∈ R is the vertical component of the vector 𝒒, and 1 3𝑁 is the 3𝑁 × 3𝑁 identity matrix. Note that since the basal boundary condition is set only on spheres vertical component, the particles belonging to J are free to move on the horizontal plane.

The square matrix 𝚲 is positive semi-definite and the vibrational (eigen) modes of the entire granular mass are directly given by its eigenvalues. Consequently, Eq. (2.10) is ill-posed when 𝜔 2 is one of its eigenvalues generating eventually resonance effects. Nevertheless, except in these situations, one has the existence and the uniqueness of the solution given by 𝒆. Figure 3 shows the histogram of the system normal (eigen) modes computed for a plane slope 𝜃 = 14 • and a layer height 𝐻/𝑑 = 14.4 at 𝑡 = 1.2 s of simulated flow time-a typical set of parameters used for our simulations presented in Section 3.

Note that to mimic the sound speed 𝑐 ≃ 10 m s -1 (and associated effective contact stiffness) observed in the experiments mentioned above, we have adjusted empirically the contact coefficients 𝜅 𝑖 𝑗 (as a kind of fit parameter) in our simplified Hertz model (see Appendix A.1).

Interparticle friction reduction through acoustic lubrication

Let us outline here the coupling between the grain-motion model cocd (see Section 2.1) and the vibration model (Eqs. (2.10)-(2.11)), and the role played by the acoustic lubrication given by the Mindlin model. More precisely, the purpose of the vibration model associated to the Mindlin model is to compute a new coefficient of friction for each contact, which embedded the vibrational perturbations into the dynamics that occurs at the grain-motion time scale. We consider that we are in the framework described in Section 2.2.2, where a mono-frequency sinusoidal vibration is imposed on the basal grains of a pile. In the following, we consider a temporal discretization of the grain-motion model; refer to Martin et al. (2023a) for the full description of the numerical scheme that is used for computing a numerical solution to cocd.

The coupling algorithm that we present in this article consists in running one iteration of the numerical scheme used to compute the solution of the grain-motion model with a typical time step Δ𝑡 𝑔 = 1 ms, providing the current generalized position vector 𝒄 ∈ R 3𝑁 and a set of the normal force intensities 𝑓 𝑖 𝑗 𝑛 , 1 ⩽ 𝑖 < 𝑗 ⩽ 𝑁. We then compute the solution 𝒒 ∈ R 3𝑁 to the vibration model. From the basal frequency 𝑓 and infinitesimal amplitudes 𝑈 0 , the Mindlin model gives us new values of the interparticle friction coefficients 𝜇 𝑖 𝑗 that are used in the Coulomb's law of the next iteration of the grain-motion model; see Eqs. (2.5)-(2.6). In the following, we describe more precisely the way the friction coefficients are modified.

When the vibration model is solved, the vector 𝑈 0 𝒒 ∈ R 3𝑁 corresponds to the ultrasoundinduced perturbations applied on the generalized position vector 𝒄 ∈ R 3𝑁 at each contact int and 𝑡 70 int are the average computational times required to complete one time integration-a simulated time of Δ𝑡 𝑔 -without ( 𝑓 = 0 kHz) and with vibrations ( 𝑓 = 70 kHz). The last column is the computational time (in percentage) spent for the calculation of vibrations and the reduction of the interparticle friction coefficients during one time integration. The MOSEK's tolerance parameter is set at the default value (10 -8 , see details in Martin et al. (2023a)). The simulations were performed on one Intel ® Core tm i7-1065G7 CPU @ 1.30 GHz × 8.

for the basal amplitude 𝑈 0 and frequency 𝑓 . From there, the normal (respectively tangential) infinitesimal displacements at contact is given by 𝑈 0 N 𝑖 𝑗 𝒒 ∈ R 3𝑁 (respectively 𝑈 0 P 𝑖 𝑗 𝒒 ∈ R 3𝑁 with P 𝑖 𝑗 introduced in Section 2.1). The question now arises of how to account for the infinitesimal perturbations of grains onto their macroscopic motion. The observations made in Leópoldès et al. (2020) show that the ultrasounds lead to modifying the static friction coefficient at the contact. That is why we choose to model the feedback of the acoustic waves on grains motion through the modification of the interparticle friction coefficient 𝜇 𝑝 , which is involved only in the tangential contact law (𝜇 𝑖 𝑗 in Coulomb's law) at the grain-motion time scale. Without the vibrations, we consider for the sake of simplicity that there is the same static friction coefficient for every contact but having different values of 𝜇 𝑝 does not change the method. More precisely, the static friction coefficient 𝜇 𝑠 includes both the interparticle friction 𝜇 𝑝 and the geometric trapping 𝜇 𝑔 (dilatancy effect). Because of the small amplitude of ultrasound, we assume that the sound-matter interaction only modifies 𝜇 𝑝 but not 𝜇 𝑔 , hence the new vibration-induced static friction coefficient 𝜇 ★ 𝑠 is resulting from the modification of 𝜇 ★ 𝑝 ≠ 𝜇 𝑝 , only. It means that due to ultrasound propagation, the static friction coefficient 𝜇 𝑖 𝑗 changes depending on perturbation amplitudes between grains. We then denote the vibration-induced interparticle static friction coefficient 𝜇 𝑡 is the microscopic oscillating tangential force intensities. Accordingly, we adopt here this scaling formula for the granular layer to describe the acoustic lubrication of the interparticle friction coefficient

𝜇 ★ 𝑖 𝑗 𝜇 𝑝 = 1 -𝐶 𝜇 𝛿f 𝑖 𝑗 𝑡 𝜇 𝑝 𝑓 𝑖 𝑗 𝑛 , (2.12)
where 𝐶 𝜇 ⩾ 0. The remaining task is then to compute an estimation of ultrasound-induced oscillating tangential force intensities 𝛿f 𝑖 𝑗 𝑡 . We approximate these intensities through a linear elastic law, meaning that they are computed thanks to the microscopic tangential grains displacements,

Computational time efficiency

At each time integration, the numerical scheme used to compute a numerical approximated solution to the grain-motion model requires solving a convex optimization problem. Additionally, the calculation of vibrational (eigen) modes and the modified friction coefficients are incorporated at each iteration. We utilize the MOSEK APS (2010) solver to address the optimization problem-see details in Martin et al. (2023a)-and for computing vibrational modes. Consequently, MOSEK is called twice an iteration. Note that thanks to our two-time scale approach, we do not need to reduce the time step of the numerical method (i.e., Δ𝑡 𝑔 ) when integrating the influence of vibrations in our simulations. We maintain the same time step in cases with or without ultrasound. Computing the modes and new friction coefficients extends the total duration of each time integration. Table 1 provides some statistics on the computational times required for performing our simulations presented in Section 3. In particular, we observe that the proportion of time spent on modes and coefficients calculation decreases as the number of grains 𝑁 increases (see the last column). For example, considering the vibrations increases the computation time of one iteration by 49 % for 𝑁 = 8000 grains, whereas it only increases by 26 % for 32 000 grains. This result is promising because it shows that our method of time-scale separation allows us to account for the rheological modification of the flow by ultrasound with a reasonable additional computation time, especially as the number of grains is large.

Results

We now investigate the transition of the granular flow from the jammed solid state induced either by increasing shear or by basal ultrasound vibrations. The simulations presented in this section were all conducted in 2D, while we introduced the equations of the coupled model in 3D in Section 2. Therefore, the 2D static granular medium consists of layers of spherical grains (diameter 𝑑) of various mean thicknesses (from 𝐻 = 3𝑑 to 14.4𝑑) put on an inclined plane at a slope 𝜃. The lengths of the samples are much larger than the thicknesses with free-stress boundary condition at both edges (Fig. 5). To observe the onset of grain motion, we focus on a specific observation domain, located at the middle of the full domain from 0.0 to 0.4 m (see the dashed rectangle in Fig. 5).

During a simulation, the slope of the basal plane is incrementally inclined (with a step Δ𝜃 = 0.5 • between two consecutive slopes) for increasing the shear. At each new slope, small amounts of grains may rearrange in the granular layer, increasing the kinetic energy 𝐸 𝑘 of the system and and the energy ratio 𝐸 𝑘 /𝐸 𝑡 with 𝐸 𝑡 the total energy. The necessary condition to reach a new stable equilibrium (state) when the slope is increased from 𝜃 to 𝜃 + Δ𝜃 is that this energy ratio satisfies the condition 𝐸 𝑘 /𝐸 𝑡 < 𝜀 𝐸 , during a simulated time of 0.1 s (equivalent to 100 iterations, since the grain motion time step is Δ𝑡 𝑔 = 1 ms). We used the value 𝜀 𝐸 = 3 × 10 -8 , which has been empirically determined during pretests. When a portion of the mass detaches (and then 𝐸 𝑘 /𝐸 𝑡 ≫ 𝜀 𝐸 ), we check whether it is only a local rearrangement, e.g., as seen in Fig. 6(a) around 𝑥 = 0.3 m at 𝑡 = 0.2 s (which keeps moving slightly at 0.4 s but arrest flow at 0.6 s), or it leads to a generalized flow in space (Fig. 6(c)).

In cases where the detachment (rearrangement) arrests on its own (i.e., 𝐸 𝑘 /𝐸 𝑡 < 𝜀 𝐸 during 0.01 s), we continue the simulation by increasing the slope until there is a continuous flow. More specifically, in the case where grain motions triggered the ultrasound spread throughout the domain (𝐸 𝑘 /𝐸 𝑡 remains greater than 𝜀 𝐸 ), we simulate the flow for about 1.2 s without changing the slope.

During the simulated time, the flow depth is considered to be split into two layers of equal depth (index 1 refers to the top surface layer and index 2 to the bottom layer). For each layer, we compute the normalized standard deviation 𝜎 𝑙 , with 𝑙 = 1, 2 of grains' velocity relative to the average velocity (Figs. 6(b,d)), i.e.

𝜎 𝑙 = 1 v𝑥𝑙 1 𝑁 𝑙 𝑁 𝑙 ∑ 𝑖=1 (𝑣 𝑖 𝑥 -v𝑥𝑙 ) 2 , (3.1)
where 𝑁 𝑙 ∈ N is the number of particles belonging to the layer 𝑙 in the observation domain (Fig. 5), 𝑣 𝑖 𝑥 is the downslope velocity component of the disk 𝑖, and v𝑥𝑙 is the mean downslope velocity of the layer 𝑙, i.e., v𝑥𝑙 = (1/𝑁 𝑙 ) 𝑁 𝑙 𝑖=1 𝑣 𝑖 𝑥 . If both 𝜎 𝑙 values converge to a limit below than two specified criteria 𝜀 𝜎 𝑙 , 𝑙 = 1, 2 (determined empirically), the flow is considered as uniform (continuous flow) and we determine the angle as the avalanche angle and denote it by 𝜃 𝑚 (see, e.g., Fig. 6(d) at 𝑡 = 0.8 s). We then assume the steady regime of inertial flow is reached when the two 𝜎 𝑙 are lower roughly than 0.3 for the top layer (𝜀 𝜎 1 = 0.3, vertical dashed line) and 0.9 for the bottom layer (𝜀 𝜎 2 = 0.9, vertical plain line). Once the avalanche angle 𝜃 𝑚 is reached, other simulations are conducted for increasing angles larger than 𝜃 𝑚 , up to 22 • , starting each time from the final configuration of the largest angle where the layer remains static (i.e., 𝜃 𝑚 -Δ𝜃).

We start by considering the case where the flow is solely induced by gravity-driven shear, i.e., by the inclination of the plane (Section 3.1), that is, without any vibration of the basal plane. Secondly, we present the effect of triggering and the dynamics generated by the dual effect of the slope inclination and basal vibration (Section 3.2).

Flow onset induced by gravity-driven shear: delay time to homogeneous flow

Figures 6(a)-(d) depict grain motion in the zoomed regions for 𝐻 = 14.4𝑑, at inclination angles close below (16 • ) and equal to the avalanche angle 𝜃 𝑚 (17 • ). When the inertial flow is initiated, a delay time is clearly observed before all grains are in motions (downslope flow) initiated from the low edge (on the right) (Figs. 6(a,c)). Comparing the grains' velocity norm at 𝜃 = 16 • and 17 • at 𝑡 = 1.2 s, we observe that the moving part of the domain (the low edge in Fig. 6(a) and the full domain in Fig. 6(c)) is intuitively 40 % slower at 16 • than at 17 • . Additionally, Fig. 6(c) shows that in the first moments, the flow is initiated not only from the low edge but also in the middle of the domain; see, e.g., the motion started by the collection of particles centered in 𝑥 = 0.28 m at 0.2 s and the motion initiated from the high edge (on the left) at 0.4 s (Fig. 6(c)). At 0.8 s, these two moving assembly of particles have merged (or percolated; see Section 4.2), finally forming a uniform flowing mass. The normalized standard deviations 𝜎 𝑙 of the grain velocity at the top (blue) and the bottom (red) layers, respectively, are shown in Figs. 6(b,d). These simulations illustrate that the avalanche angle is well around 𝜃 𝑚 = 17 • . Indeed, at 0.8 s, the 𝜎 𝑙 drop below the two criteria (Fig. 6(d)). Conversely, it is clear that at 16 • the normalized standard deviations 𝜎 𝑙 are far above the criteria for uniform flows (Fig. 6(b)).

Even when the whole flow reaches a high level of uniformity (both curves representing 𝜎 𝑙 drop below the criteria), the flow uniformity in the top layer is larger than in the bottom layer, as the blue curve consistently remains below the red curve in Fig. 6(d). This behavior is consistently observed in all our simulations whatever the layer thickness and is consistent with the laboratory experiments of [START_REF] Bachelet | Acoustic emissions of nearly steady and uniform granular flows: A proxy for flow dynamics and velocity fluctuations[END_REF] showing larger velocity fluctuations at the base of the flow.

Figure 7 highlights the avalanche angles for different thicknesses in the configuration where motion is solely generated by the inclination of the plane, i.e., without any basal vibrations. These angles represent the transition from a state where the granular mass is at rest (jammed solid state) to a continuous flow (for inclinations greater than the avalanche angle), sometimes passing from a metastable state that is not represented here. This boundary between the two limit states is depicted by the blue curve in a phase space composed, on the 𝑥-axis, of the inclination of the basal plane and, on the 𝑦-axis, of the average thickness of the granular mass in the initial state.

The shape of this blue curve (monotonic and decreasing) is similar to those reported in the literature [START_REF] Daerr | Two types of avalanche behaviour in granular media[END_REF][START_REF] Forterre | Stability analysis of rapid granular chute flows: formation of longitudinal vortices[END_REF][START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF] showing that as the initial thickness increases, the avalanche angle decreases. The gray area represents the graph zone where the results strongly depend on the boundary conditions because the thickness of the granular layer is very thin (like 𝐻 = 3𝑑).

From angles greater than the avalanche angle 𝜃 𝑚 , we measure the time 𝑡 𝑐𝑜𝑛 it takes for 𝜎 𝑙 in the top and bottom layers to satisfy the criteria characterizing uniform flow. These times are represented in Fig. 8 for the three heights 𝐻 = 6𝑑, 11.4𝑑, and 14.4𝑑. The times 𝑡 𝑐𝑜𝑛 (without basal vibrations) represented by the blue curves are all decreasing. As expected, the larger the slope angle is, the shorter is the time required for the grains reach a uniform flow.

It is worth noting that the triggering duration decreases with the slope angle, but it seems to converge to a limit beyond a certain slope angle, possibly due to the choice of the imposed criteria. This is reflected in Fig. 8 where the curves become almost horizontal for large angles. For example, for 𝐻 = 14.4𝑑, if the flow becomes uniform 0.5 s earlier by increasing the slope angle from 17 • to 18 • , it only becomes uniform 0.02 s earlier by increasing the slope angle from 21 • to 22 • . Finally, one might have expected a clear relationship between the inclination angle 𝜃 and the time 𝑡 𝑐𝑜𝑛 with the thickness of the granular layer 𝐻, but this is not what we observe here. In fact, while the thickness 𝐻 = 6𝑑 has a slower triggering time than the other two at 20 • , indicating that the thinner thickness is slower, this is not the case at 17 • , where the triggering for the largest thickness 𝐻 = 14.4𝑑 is faster than that for 11.4𝑑.

Flow onset induced by vibration: transition to uniform flows via percolation

In this section, we investigate the triggering of granular flows in the presence of basal vibrations at a frequency of 70 kHz. As described before, this triggering is modeled by coupling the cocd model (Section 2.1) with the steady vibration model (Section 2.2.2), through the modification of interparticle friction coefficients 𝜇 𝑝 , which are altered by the vibrations using the Mindlin model (Section 2.3). More specifically, we observe how this basal vibration affects the results obtained without basal vibrations. Note that the basal vibration frequency 𝑓 = 70 kHz is relatively high for a system with normal modes ranging from 1 kHz to 75 kHz, as shown in Fig. 3. We discuss this choice in Section 4.3.

Figures 9(a,b) present the same types of results as Figs. 6(a,b), decomposed into simulation snapshots and curves representing the temporal evolution of the 𝜎 𝑙 in the top and bottom layers, still for a thickness of 𝐻 = 14.4𝑑. The difference here is that the basal vibration is on. Thus, at 𝜃 = 16 • , and contrary to the case without vibrations presented in Figs. 6(a,b), it is observed that a flow of the granular layer is triggered and converges towards a continuous flow after approximately 𝑡 = 0.25 s (Fig. 9(a)). There is therefore a notable difference since only a part of the layer was flowing after 1.2 s in the case without vibrations (Figs. 6(a,b)). As a result, the avalanche angle 𝜃 𝑚 is clearly lower in the case of basal vibrations compared to the case without. Similarly to the 17 • angle without vibrations (Figs. 6(c,d)), the "top" layer is more uniform than the "bottom" layer, as indicated by the blue curve being significantly below the red curve (Fig. 9(b)). Unlike in the case without basal vibrations, we do not observe the initiation of flow occurring primarily at the front of the layer (low edge on the right in Fig. 6(c)). The mobilization of the granular layer is observed to be more uniform in the presence of vibrations. Another way to understand this is that the distinct triggering zones, which remain isolated for a considerable duration vibrations off, percolate much faster with vibrations on (see Section 4).

The decrease of avalanche angles due to vibration is observed in all simulations whatever the thickness and slope inclination, as shown in Fig. 7. Indeed, the orange curve, like the blue curve, is decreasing, but the important result here is that each of the avalanche angles in the orange curve (with vibrations) is lower than those in the blue curve by approximately 2 • . The only notable singularity is the avalanche angle for the thin thickness 𝐻 = 3𝑑, which is only 14.5 • , which is 5.5 • smaller than in the case without basal vibrations. Such a difference can be explained by the dominance of boundary conditions as discussed in Section 4.

Let us now investigate what happens at smaller angles, e.g. at 13 • but still for the thickness 𝐻 = 14.4𝑑, (Fig. 9(c,d)). Particularly visible in the snapshots (Fig. 9(c)), we observe that rearrangements are active, even at 𝑡 = 1.2 s (some areas are light blue on the surface after 1.2 s). However, these rearrangements remain local and are not sufficient to trigger a flow that would develop into a continuous flow. This is confirmed by the values of 𝜎 𝑙 , which both remain well below the criteria characterizing a uniform flow (Fig. 9(d)).

When the avalanche angle is exceeded, the effect of basal vibrations is to decrease the triggering time, as shown in Fig. 8. Indeed, for the same layer thickness, the triggering times represented by the orange curves are consistently lower than the triggering times without vibrations, i.e., the blue curves in Fig. 8. For example, for a thickness of 𝐻 = 6𝑑 and an angle 𝜃 = 19 • , the triggering time decreases by approximately 74 % (from 𝑡 𝑐𝑜𝑛 = 1.15 s without vibrations to only 0.3 s with vibrations). Similarly to the case without vibrations, the triggering times appear to converge to a limit but for angles much smaller (for 𝐻 = 11.4𝑑 and 14.4𝑑, the triggering times with vibrations are roughly the same at 18 • as those at 22 • without vibrations). Note that, once again, there is no clear influence of thickness on the triggering time.

Discussion

Force chains, ultrasound-induced deformation, and lubricated contacts

In the previous Section 3.2, we quantified how much basal ultrasonic vibrations reduces the avalanche angle and triggering time. We focus here on the interparticle mechanisms responsible for these effects, namely, the vibration-induced lubrication which is accounted for in our simulation by weakening of interparticle friction coefficients, as described by the Mindlin model(Section 2.3).

Figure 10 shows the grain configuration, at two consecutive times in the simulation 𝑡 = 0.051 s and 0.052 s (the numerical time step is Δ𝑡 𝑔 = 1 ms) for a slope angle 𝜃 = 15 • , a layer thickness 𝐻 = 14.4𝑑 and in a region between 𝑥 = 0.497 and 0.505 m. The black lines in Figs. 10(a,d) represent the force chains formed between the grains (their width correspond to normal force intensities 𝑓 𝑖 𝑗 𝑛 (see Section 2.1)). The arrows shown in Figs. 10(b,e) represent the deformation fields computed with the steady vibration model (the vector 𝒒 ∈ R 3𝑁 solution to Eqs. (2.10)-(2.11)). Note that the layer of grains in contact with the bottom has a unit displacement of 1, but for the sake of visibility of the deformation field, we did not use the same scale in Figs 10(b,e), 12(b), and 13(b, e, h, k). The force chains are again depicted in Figs. 10(c,f), but on a logarithmic scale (which explains the different thicknesses compared to the linear scale in Figs. 10(a,d). The colors represent the rate of change of interparticle friction coefficients, i.e., the quantity (𝜇 𝑝 -𝜇 ★ 𝑖 𝑗 )/𝜇 𝑝 (see Eq. (2.12)). Where the force chains are blue, the coefficients are slightly modified, while where they are red, the coefficients are reduced to zero (thus decreased by 100 %). Figure 10 illustrates two destabilization mechanisms of the granular layer revealed by our simulations.

The first mechanism, presented in Figs. 10(a)-(c), involves the transmission of deformation induced by basal vibrations in a preferential direction (in this case, horizontal around 𝑧 = 0.005 m and between 𝑥 = 0.496 and 0.501 m) through a relatively dominant (strong) force chain in the corresponding region (thicker force chain in this particular region (Fig. 10(a))). This transmission only slightly affects the interparticle friction coefficients along this force chain (the horizontal force chain appears as dark blue in Fig. 10(c)). On the contrary, the friction coefficients on contacts relatively perpendicular to this main force chain significantly changed (bright colors for vertical force chains surrounding the main force chain in Fig. 10(c)).

A second mechanism revealed by our simulations involves a kind of block destabilization. This mechanism is presented in Figs. 10(d)-(f). In this case, the destabilization does not follow a preferred direction but occurs in multiple directions (Fig. 10(e)), while remaining localized in a specific area of the granular layer (between 𝑧 = 0.002 and 0.006 m and 𝑥 = 0.496 and 0.501 m in Fig. 10(d)). A few strong force chains are present in this area (Fig. 10(d)), but the friction coefficients are primarily modified on the less dominant (weak) force chains within the disturbed zone (Fig. 10(f)). Note that the different response of weak and strong force chains when getting close to instability was also observed in [START_REF] Deboeuf | Memory of the unjamming transition during cyclic tiltings of a granular pile[END_REF].

In our simulations, we have noticed that the first mechanism (destabilization along a preferred direction) can trigger the second mechanism, as seen here since the figures are taken at consecutive times. However, this is not always the case. Sometimes, only the first or the second mechanism occurs. We have not observed that the second mechanism can trigger the first mechanism.

Nucleation time (delay) to reach uniform flows

When the basal plane is inclined, localized rearrangement zones appear both with and without vibrations (Figs. 6(a,c), and 9(a, c)). Without vibrations, these zones can remain isolated for a significant period if the slope is not too steep (Fig. 6(a)). However, beyond the avalanche angle, they can eventually merge (or "percolate") and result in a relatively uniform flow (Fig. 6(c)). In simulations with vibrations, in addition to this triggering mechanism through inclination, there are also interparticle vibratory mechanisms that trigger specific zones within the granular layer (Fig. 10). The difference with vibrations is that the triggering zones are much more numerous and percolate much faster (Figs. 9(a) and 8), leading to a quick homogenization of the flow (Fig. 9(b)) and resulting in lower avalanche angles compared to cases without vibrations (Fig. 7). In our simulations, this phenomenon of percolation transition is faster and more uniform with vibrations on, and it increases with the slope, as shown in Fig 11.

Basal vibration frequency and amplitude

The ultrasound basal vibration has a quite high frequency ( 𝑓 = 70 kHz) when compared to the system normal modes, which are mostly comprised between 5 kHz to 45 kHz (Fig. 3). In our simulations, the disturbance is more effective when the frequency is well-centered within the (2.12)). The basal vibration amplitude is the same as in Fig. 10, about 𝑈 0 /𝑑 = 10 -5 . medium's spectrum and the results presented are qualitatively similar at a lower frequency. For example, the simulation we conducted with 𝑓 = 30 kHz shows a more significant and less localized disturbance of friction coefficients (Fig. 12). Figure 12 highlights that the fundamental disturbance was transmitted more effectively from the bottom of the medium to the grains in the upper layers, and this was achieved more efficiently than at the higher frequency of 70 kHz (Fig. 10).

Regarding the choice of the amplitude of the basal vibration 𝑈 0 , we used the same amplitude for all simulations, 𝑈 0 /𝑑 = 10 -5 , which is equivalent to a few tens of nanometers. Figure 13 illustrates the effects of the ultrasonic vibrations on the coefficients of friction for the same frequency, 𝑓 = 70 kHz, but with amplitude 𝑈 0 ranging from 𝑈 0 /𝑑 = 10 -3 to 𝑈 0 /𝑑 = 10 -6 . It can be observed that, even though the solution of the steady problem remains similar in terms of the distribution of perturbations (Figs. 13(b,e,h,k)), simulations conducted with larger amplitudes significantly alter the friction coefficients, obviously more than those with smaller amplitudes (Figs. 13(c,f,i,l)). Furthermore, the choice of amplitude made in this article corresponds to an amplitude so small that the perturbations, although non-zero, remain small and generally localized (Fig. 13(i)).

By making these choices of frequency and amplitude of the ultrasound, we deliberately positioned the numerical simulations in an amplitude-frequency regime that makes disturbing the medium the most challenging, without inducing significant rearrangement of grain positions. The high-frequency domain is less explored than that of low-frequency vibrations, which are already well-known for their significant triggering effects, e.g., see [START_REF] Hanotin | Vibration-induced Liquefaction of Granular Suspensions[END_REF][START_REF] Lastakowski | Granular friction: Triggering large events with small vibrations[END_REF]. Decreasing (respectively increasing) the value of the basal vibration frequency (respectively amplitude) does not qualitatively alter the results presented in this article. The coupling of numerical models we present in this article is relevant in a frequency range from 5 kHz to 80 kHz (Fig. 3) and with amplitudes ranging from 10 -6 𝑑 to 10 -2 𝑑. We have not conducted studies beyond this range of parameters, but our model is likely to be less relevant in those cases since consideration of grain displacement at very low frequencies and/or large amplitudes may be necessary; see, e.g., [START_REF] Baumberger | Solid friction from stick-slip down to pinning and aging[END_REF][START_REF] Bureau | Jamming creep of a frictional interface[END_REF][START_REF] Hanotin | Vibration-induced Liquefaction of Granular Suspensions[END_REF]). In addition, in our simulations, the front systematically starts moving when the avalanche angle is exceeded (as illustrated in the inset Fig. 14(c), which shows the granular front, i.e., the right end of the mass in Fig. 5). Similarly, the experimental data also exhibit a configuration of a granular mass flowing near the front (Fig. 14(a)). In both experiments and simulations, we observe that vibration-triggered flows are more uniform (Figs. 14(b,d)) than the ones induced by gravity only (Figs. 14(a,c)).

Thus, we believe that the additional triggering mechanisms associated with basal vibrations, even if they are infinitesimal (ultrasounds in this article and in the experiments), are sufficient to uniformly destabilize the granular layer and consequently lead to a more uniform destabilization, resulting in a reduction of the avalanche angle of the mass.

Conclusion

In summary, we have developed a two-time scale numerical simulation of wave-induced friction weakening in 2D granular layers. Indeed ultrasound vibrations propagate with a time-scale of the order of 10 microseconds while grain motion occur at a time-scale of milliseconds. The triggering of granular flows is modeled by coupling the cocd model (Section 2.1) with the steady vibration model (Section 2.2.2), through the modification of interparticle friction coefficients 𝜇 𝑝 , which are altered by the vibrations using the Mindlin model (Section 2.3). This new two-time scale model has allowed us to investigate the nonlinear interaction between ultrasound and granular flows, in particular the vibrationinduced reduction of the interparticle friction coefficient through the acoustic lubrication, without the contact opening and rearrangements of grain positions.

As expected, we find that ultrasound vibration is predominantly supported by the strong-force chains, but the vibration-induced decrease of friction occurs mainly in the weak-force chains perpendicular to the strong contact forces, causing eventually STZ-like zones (shear transformation zones) on the mescoscopic scale. These local rearrangements nucleate to create a continuous flow through a percolation process with certain delay depending on the proximity to the failure, i.e. avalanche angles. The larger the vibration amplitude (or lower excitation frequency) is, the stronger the ultrasound induced destabilization is. The different in the behaviour of weak and strong force chains was also observed in the response of a granular packing (without vibration) during successive loading-unloading cycles close to the avalanche angle [START_REF] Deboeuf | Memory of the unjamming transition during cyclic tiltings of a granular pile[END_REF].

Compared to gravity-driven flow, ultrasound-induced flow appears more spatially homogeneous. This is consistent with the effective temperature role played by sound vibration. The qualitative agreement between these simulations and experimental observations of granular flows triggered by ultrasound supports our numerical modeling. Although further improvement on the vibration model is still needed, this work helps to highlight underlying mechanisms of landslide and earthquake triggering by seismicity; see, e.g., [START_REF] Durand | Repetitive small seismicity coupled with rainfall can trigger large slope instabilities on metastable volcanic edifices[END_REF].

where 𝜅 𝑖 𝑗 is a constant depending on grains' properties and 𝛿 𝑖 𝑗 is the overlap between the particles 𝑖 and 𝑗, characterizing grains deformation.

Note the difference between the distances D 𝑖 𝑗 , defined in Eq. (2.1) and 𝛿 𝑖 𝑗 : on the one hand, the distance D 𝑖 𝑗 measures how far the grains are from each other, on the other hand, the overlap 𝛿 𝑖 𝑗 models an overlap between the grains 𝑖 and 𝑗, which should be seen as quantifying a deformation of the bodies now considered as elastic. The two terms D 𝑖 𝑗 and 𝛿 𝑖 𝑗 evolve in opposite way. When the distance D 𝑖 𝑗 between 𝑖 and 𝑗 increases their overlap 𝛿 𝑖 𝑗 decreases. The scalars 𝑓 𝑖 𝑗 𝑛 can be provided by any model like cocd, as being the normal intensity of the contact force, while 𝐹 𝑛 is the normal force intensity in the Hertz theory framework.

The prefactor can be found from the Hertz theory of elastic contact (see the book of [START_REF] Andreotti | Granular Media: Between Fluid and Solid[END_REF] and the book of [START_REF] Johnson | Contact Mechanics[END_REF]). In the case of a contact between two spheres 𝑖 and 𝑗 of radii 𝑟 𝑖 , 𝑟 𝑗 , Young's moduli 𝐸 𝑖 , 𝐸 𝑗 , and Poisson's ratios 𝜈 𝑖 , 𝜈 𝑗 , we have 𝜅 𝑖 𝑗 = (4/3)𝐸 𝑖 𝑗 √ 𝑟 𝑖 𝑗 , where 𝐸 𝑖 𝑗 and 𝑟 𝑖 𝑗 are defined by 1/𝐸 𝑖 𝑗 = (1 -𝜈 2 𝑖 )/𝐸 𝑖 + (1 -𝜈 2 𝑗 )/𝐸 𝑗 and 1/𝑟 𝑖 𝑗 = 1/𝑟 𝑖 + 1/𝑟 𝑗 . However, to compare qualitatively the sound speed or elastic modulus of granular layers measured in laboratory experiments (see discussions in Section 2) with the present simplified Hertz (normal) contact model, we would need to upscale uniformly the coefficients 𝜅 𝑖 𝑗 by a factor of 100. This scaling does not qualitatively affect the results presented in this article since the basal vibration frequency is chosen with respect to the system normal modes; see Section 2.

A.2. First order Taylor expansion of Hertz contact force

Let 𝐽 𝑐 be the set of contacts defined by 𝐽 𝑐 = {(𝑖, 𝑗) | 1 ⩽ 𝑖 < 𝑗 ⩽ 𝑁 and ∃ 𝑓 0 𝑖 𝑗 𝑛 > 0}. We denote by 𝑁 𝑐 the cardinal of set 𝐽 𝑐 . Let us assume to be given a generalized normal force intensity vector 𝒇 0 ∈ R 𝑁 𝑐 , where 𝑓

A.3. Perturbation-induced position vector

We now consider given a generalized position vector of the mechanical system at the equilibrium that we denote 𝒄 0 ∈ R 3𝑁 ; see Section 2.1. In the previous section, the perturbation-induced force 𝑓 𝑖 𝑗 𝑛 can be expressed as a linear function of the infinitesimal perturbation 𝜀 𝑖 𝑗 of the overlap 𝛿 0 𝑖 𝑗 (see Eq. A 3). Similarly, we define the perturbationinduced position vector 𝒄 and 𝒆, the infinitesimal perturbation of the constant position vector 𝒄 0 , all belonging to R 3𝑁 , i.e., 𝒄 = 𝒄 0 + 𝒆, where, at the vibration time scale, 𝒄 and 𝒆 depend on time while 𝒄 0 remains constant. We can write 𝜀 𝑖 𝑗 , the infinitesimal perturbation of the overlap 𝛿 0 𝑖 𝑗 as the image of 𝒆 by the map N 𝑖 𝑗 , indeed, we have We define also the perturbation-induced generalized force intensity vector 𝒇 ∈ R 𝑁 𝑐 , resulting from the perturbation of the constant vector 𝒇 0 ∈ R 𝑁 𝑐 by the infinitesimal perturbation vector 𝒇 𝜀 ∈ R 𝑁 𝑐 , i.e., 𝒇 = 𝒇 0 + 𝒇 𝜀 .

Equations (A 2), (A 4), and (A 5) give 𝑓 𝜀 𝑖 𝑗 𝑛 = -3/2(𝜅 𝑖 𝑗 ) 2/3 𝑓 0 𝑖 𝑗 𝑛 1/3 N 𝑖 𝑗 𝒆, which can be written under its vector form as

𝒇 𝜀 = -KN 𝒆 ∈ R 𝑁 𝑐 ,
where the diagonal square matrix K contains the elastic properties of the system: ∈ R 𝑁 𝑐 × 𝑁 𝑐 ,

A.4. Wave equation

The assumption is made on the system described in Section 2.1 to be at rest, maintained by the gravity field. In this configuration, Newton's second law imposes that there exists a generalized reaction force vector that is necessary opposed to the global force vector 𝒘 0 = (𝒘 0 1 , . . . , 𝒘 0 𝑁 ) ∈ R 3𝑁 , applied on the system. Coupled with the non-overlapping condition D 𝑖 𝑗 ⩾ 0, it is equivalent to say that the inverse image of the vector -𝒘 0 by the map N 𝑇 (see Section 2.2.1), intersected with R 𝑁 𝑐 + (to have repulsive force only) is not empty, i.e., N 𝑇 -1 (-𝒘 0 ) ∩ R 𝑁 𝑐 + ≠ ∅. Let the generalized normal force intensity vector 𝒇 0 ∈ R 𝑁 𝑐 + , belonging to the inverse image of the vector -𝒘 0 by the map N 𝑇 , we have

N 𝑇 𝒇 0 = -𝒘 0 ∈ R 2𝑁 ,
whose local expression is exactly given by Eq. (2.7). The generalized force vector N 𝑇 𝒇 0 generates local normal repulsive forces between the spheres at the contact points. Furthermore, the generated pressures create deformation zones and the Hertz theory enables us to characterize these deformations by a set of scalars 𝛿 0 𝑖 𝑗 ∈ R (see Eq.(A 1)). Consequently, the 𝛿 0 𝑖 𝑗 can be given by an equation of type Eq. (A 2). As a result, we have shown how the forces computed at the grain-motion time scale by the grain-motion model cocd, which does not involve the Hertz theory can be linked to the latest at the vibration time scale, through the introduction of the overlaps 𝛿 0 𝑖 𝑗 . With the generalized mass matrix (masses only) defined by Eq. (2.8), and by considering the perturbation-induced generalized force vector N 𝑇 𝒇 ∈ R 3𝑁 , and applying Newton's second law, we finally obtain the unsteady wave equation resulting from the perturbation of 𝒄 0 by 𝒆. Indeed, we have At the end, a wave equation is defined for any 𝒆 ∈ R 3𝑁 M d 2 𝒆 d𝑡 2 + 𝚲𝒆 = 0, where the linear map defined by the matrix 𝚲 = N 𝑇 KN ∈ R 3𝑁 ×3𝑁 , can be seen as a kind of discrete Laplace operator.

HFigure 1 :

 1 Figure 1: (a) 2D schematic illustration of granular flows triggered by small-amplitude ultrasonic or seismic vibration (indicated by double arrows) where the granular layer of thickness 𝐻 is deposited on a slope at angle 𝜃 below the maximum angle of stability 𝜃 𝑚 .The inertial flow triggered by ultrasonic vibration is mostly uniform and continuous.(b) Sketch of the normalized stress-strain rate relation in a sheared granular medium: static state (V = 0), unstable flow (velocity-weakening) and stable flow (velocity strengthening). Under a shear 𝜇 between 𝜇 𝑠 and 𝜇 𝑑 , the metastable state can be switched abruptly to a flowing state by acoustic perturbation (lubrication).

Figure 2 :

 2 Figure 2: 2D schematic depiction of the two contact laws. (a) Representation of a three-disk situation: Disks 1 and 2 are not in contact, while disks 2 and 3 are. D 𝑖 𝑗 , (1 ⩽ 𝑖 < 𝑗 ⩽ 3) indicates the normal distance between disks 𝑖 and 𝑗, 𝑓 𝑖 𝑗 𝑛 is the normal force's intensity at the contact, 𝒏 2,3 (resp. 𝒕 2,3 ) is the unit normal (resp. tangential) vector at the contact between 2 and 3, 𝒇 𝑡 2,3 denotes the tangential force vector, and 𝒗 𝑡 2,3 is the tangential relative velocity vector between disks 2 and 3. (b) Graph representing the normal law. (c) Graph depicting the Coulomb friction law. 𝒙 • 𝒚 denotes the dot product of vectors 𝒙 and 𝒚. (d) Notations in 3D.

  𝑖 𝑗 𝑛 𝒏 𝑖 𝑗 + 𝒇 𝑖 𝑗 𝑡 on particle 𝑖. The reciprocity of this contact's action on both particles is evident from the definitions of 𝑓 𝑗𝑖 𝑛 and 𝒇 𝑗𝑖 𝑡 , derived from 𝑓 𝑖 𝑗 𝑛 and 𝒇 𝑖 𝑗 𝑡 . The normal force exerted on sphere 𝑖 due to this contact is 𝑓 𝑖 𝑗 𝑛 𝒏 𝑖 𝑗 , and 𝒇 𝑖 𝑗 𝑡 ∈ Π 𝑖 𝑗 represents the frictional (tangential) force, which lies in the plane orthogonal to 𝒏 𝑖 𝑗 .

Figure 3 :

 3 Figure 3: Histogram of the vibrational (i.e., eigen) modes of the quasi-static system considered at the vibration time scale, computed for a plane slope 𝜃 = 14 • , flow height 𝐻/𝑑 = 14.4 at 𝑡 = 1.2 s. There are about 32 000 modes presented in this figure.

  𝑖 𝑗 𝑛 ⩾ 0 for 1 ⩽ 𝑖 < 𝑗 ⩽ 𝑛. These normal force intensities 𝑓 𝑖 𝑗 𝑛 are provided by the grain-motion model resolving the normal, tangential, and collision laws (see Section 2.1 and Refs. Maury (2006); Martin et al. (2023a)) and they are not modified by the vibrations at the vibration time scale.

Figure 5 :

 5 Figure 5: Full domain simulated. The dashed rectangle represent the observation domain of Figs. 6 and 9.

Figure 6 :

 6 Figure 6: Granular assembly without basal vibrations ( 𝑓 = 0 kHz) at 𝜃 = 16 • (a, b) and 17 • (c, d) for 𝐻 = 14.4𝑑. (a, c) A sequence of 7 snapshots showing the evolution of a granular flow over time, ranging from 𝑡 = 0 to 1.2 s. The color scale represents the magnitude of grain velocity. The observation domain covers a distance of 0.4 m and corresponds to the dashed rectangle in Fig. 5. (b, d) Variation, over the same time period, of the normalized standard deviation 𝜎 𝑙 of grains' downslope velocity with respect to the average downslope velocity of the flow within two different layers decomposing the flow depth.

  m

Figure 7 :

 7 Figure 7: Phase space separating jammed and continuous regimes as a function of the avalanche angle 𝜃 𝑚 and the flow thickness ratio 𝐻/𝑑. Each point represents the avalanche angle for the five thickness ratios 𝐻/𝑑 = 3, 6, 8.7, 11.4, and 14.4. The blue curves represent these avalanche angles without basal vibrations ( 𝑓 = 0 kHz), contrary to the orange ones ( 𝑓 = 70 kHz).

Figure 8 :

 8 Figure 8: Delay time of avalanche triggering by gravity or (and) ultrasound vibration as a function of the slope angle for different thickness ratios 𝐻/𝑑 = 6, 11.4, and 14.4. The blue curves correspond to the measured times for flows without vibrations ( 𝑓 = 0 kHz), in contrast to the orange curves ( 𝑓 = 70 kHz).The stars indicate the maximum angles of stability corresponding to the avalanche angles measured for different thicknesses in Fig.7whereas the filled squares point to the inclined angles for which the flow has been initiated.

Figure 9 :

 9 Figure 9: Granular assembly with basal vibrations ( 𝑓 = 70 kHz) at 𝜃 = 16 • (a, b) and 13 • (c, d) for 𝐻 = 14.4𝑑. (a, c) A sequence of 7 snapshots showing the evolution of a granular flow over time, ranging from 𝑡 = 0 to 1.2 s. The color scale represents the magnitude of grain velocity. The observation domain covers a distance of 0.4 m and corresponds to the dashed rectangle in Fig. 5. (b, d) Variation, over the same time period, of the normalized standard deviation 𝜎 of grains' downslope velocity with respect to the average downslope velocity of the flow within two different layers decomposing the flow depth.

Figure 10 :

 10 Figure 10: (a, d) Snapshots of the contact force chains, (b, e) vibration-induced displacements, and (c, f) rate of vibration-induced perturbation of the interparticle friction, for 𝜃 = 15 • , 𝐻/𝑑 = 14.4, 𝑓 = 70 kHz. In (a, d), the black lines represent the force chains formed between the grains on a linear scale (the coefficients 𝑓 𝑖 𝑗 𝑛 in Section 2.2.2). In (b, e), the black arrows represent the computed vibration-induced displacement 𝒒 in Eqs. (2.10)-(2.11) (but not at the same scale). In (c, f), the lines between grains represent the force chains on a logarithmic scale, and the color map represent the rate of decrease in interparticle friction coefficients, i.e., the quantity (𝜇 𝑝 -𝜇 ★ 𝑖 𝑗 )/𝜇 𝑝 (see Eq. (2.12)). Snapshot times are 𝑡 = 0.051 s (a)-(c) and 𝑡 = 0.052 s (d)-(f).

Figure 11 :

 11 Figure 11: The percolation transition is accelerated with the presence of basal ultrasound vibrations when increasing the slope angle.

Figure 12 :

 12 Figure 12: (a) Snapshots of the contact force chains, (b) vibration-induced displacements, and (c) rate of vibration-induced perturbation of the interparticle friction, for 𝜃 = 15 • , 𝐻/𝑑 = 14.4, 𝑓 = 30 kHz, at 𝑡 = 0.072 s. (a) The black lines represent the force chains formed between the grains on a linear scale (the coefficients 𝑓 𝑖 𝑗 𝑛 in Section 2.2.2). (b) The black arrows represent the computed vibrational displacement 𝒒 in Eqs. (2.10)-(2.11). (c) The lines between grains represent the force chains on a logarithmic scale, and the color map represent the rate of decrease in interparticle friction coefficients, i.e., the quantity (𝜇 𝑝 -𝜇 ★ 𝑝 )/𝜇 𝑝 (see Eq.(2.12)). The basal vibration amplitude is the same as in Fig.10, about 𝑈 0 /𝑑 = 10 -5 .

U 0

 0 /d = 10 -3 U 0 /d = 10 -4 U 0 /d = 10 -5 U 0 /d = 10 -6

Figure 13 :

 13 Figure 13: (a, d, g, j) Snapshots of the contact force chains, (b, e, h, k) vibration-induced displacements (normalized by the basal amplitude), and (c, f, i, l) rate of vibration-induced perturbation of the interparticle friction, for 𝜃 = 15 • , 𝐻/𝑑 = 14.4, 𝑓 = 70 kHz, and amplitudes 𝑈 0 /𝑑 = 10 -𝑛 , for 𝑛 = 3 . . . 6, at 𝑡 = 0.018 s.

Figure 14 :

 14 Figure 14: Experimental velocity fields of granular flows (side view) in a dry packing during the first 0.15 s interval: (a) driven by gravity at the avalanche angle 𝜃 𝑚 = 24 • ± 1 • . (b) Triggered by ultrasounds at 𝜃 = 22 • . (c) Simulation snapshot at the front for 𝜃 𝑚 = 17 • , 𝐻/𝑑 = 14.4, 𝑓 = 0 kHz, at 𝑡 = 0.15 s. (d) Simulation triggered by ultrasounds ( 𝑓 = 70 kHz).

  for any (𝑖, 𝑗) ∈ 𝐽 𝑐 . We can define the overlap 𝛿 0 𝑗 remain constant at the vibration time scale and we study the effects of small perturbations around this equilibrium. We denote by 𝜀 𝑖 𝑗 ∈ R, the infinitesimal perturbation of 𝛿 0 𝑖 𝑗 and the perturbation-induced overlap by 𝛿 𝑖 𝑗 ∈ R, such that we have by definition𝛿 𝑖 𝑗 = 𝛿 0 𝑖 𝑗 + 𝜀 𝑖 𝑗 . (A 3)Only 𝛿 𝑖 𝑗 and 𝜀 𝑖 𝑗 depend on time while 𝛿 0 𝑖 𝑗 remains constant. Consequently, when considering the perturbation-induced overlap 𝛿 𝑖 𝑗 , a Taylor expansion at the first order of Eq. (A 2) provides the value of a perturbation-induced normal force intensity, denoted by 𝑓 is the infinitesimal perturbation of 𝑓 0 𝑖 𝑗 𝑛 and plays for the perturbation-induced force 𝑓 𝑖 𝑗 𝑛 the role played by the term 𝜀 𝑖 𝑗 for 𝛿 𝑖 𝑗 in Eq. (A 3).

𝜀

  𝑖 𝑗 = (𝛿 𝑖 𝑗 -𝛿 0 𝑖 𝑗 ) = -((𝒄 𝑖 -𝒄 𝑗 ) -(𝒄 0 𝑖 -𝒄 0 𝑗 )) • 𝒏 𝑖 𝑗 = -𝒆 𝑖 -𝒆 𝑗 • 𝒏 𝑖 𝑗 = -N 𝑖 𝑗 𝒆. (A5)

  , (𝜅 𝑖 𝑗 ) 2/3 𝑓 0 𝑖 𝑗 𝑛 1/3 , . . . . . . , (𝜅 𝑁 -1, 𝑁 ) 2/3 𝑓 0 𝑁 -1, 𝑁 𝑛 1/3

M d 2

 2 𝒄 d𝑡 2 = N 𝑇 𝒇 + 𝒘 0 ⇐⇒ M d 2 𝒆 d𝑡 2 = N 𝑇 ( 𝒇 0 + 𝒇 𝜀 ) + 𝒘 0 M d 2 𝒆 d𝑡 2 = -N 𝑇 KN 𝒆.

Table 1 :

 1 Statistics on computational time. 𝐻/𝑑 gives the normalized thickness of the granular layer, 𝑁 is the number of grains, 𝑡 0

	𝐻/𝑑	𝑁	𝑡 0 int	𝑡 70 int	(𝑡 70 int -𝑡 0 int )/𝑡 70 int
	3	8000 1.24 s 2.44 s	49 %
	6 13 319 2.32 s 4.03 s	42 %
	8.7 20 142 5.21 s 7.91 s	34 %
	11.4 24 000 10.29 s 14.46 s	29 %
	14.4 32 000 13.04 s 17.50 s	26 %

  ★ 𝑖 𝑗 . From the Mindlin model (see[START_REF] Léopoldès | Onset of Sliding in Amorphous Films Triggered by High-Frequency Oscillatory Shear[END_REF]), we consider that the decrease rate of the interparticle friction coefficient Δ𝜇 ★ 𝑝 /𝜇 𝑝 is approximately proportional to the ratio of the microscopic oscillating tangential force to the static normal force intensities times the coefficient of friction (this product being actually the bound of the tangential force 𝒇

				𝑖 𝑗
				𝑡 (see Eqs. (2.5)-(2.6))), i.e.,
	Δ𝜇 ★ 𝑝 /𝜇 𝑝 ∝ -𝛿f	𝑖 𝑗 𝑡 /(𝜇 𝑝 𝑓	𝑖 𝑗 𝑛 ), where 𝛿f	𝑖 𝑗

Acknowledgements. We thank MOSEK APS (2010) for the free academic license that made this article possible.

Funding. This article has been funded by the ERC Contract No. ERC-CG-2013-PE10-617472 SLIDE-QUAKES.

Declaration of interests.

The authors report no conflict of interest.

Author ORCIDs. H. A. Martin, https://orcid.org/0000-0002-7089-8958; A. Mangeney, https://orcid.org/0000-0002-3197-6087; X. Jia, https://orcid.org/0000-0003-0874-8235; B. Maury, https://orcid.org/0000-0001-5474-4134; A. Lefebvre-Lepot, https://orcid.org/0000-0002-9521-126X; Y. Maday, https://orcid.org/0000-0002-0443-6544; P. Dérand, https://orcid.org/0009-0003-0166-006X Funding acquisition, Methodology, Project administration, Supervision, Writing -original draft. P. Dérand: Conceptualization, Methodology, Experiment.

Appendix A. Wave equation

The purpose of this section is to present the derivation of the wave equation introduced by Eq. (2.9). We model the time-evolution of an infinitesimal perturbation of a given configuration 𝒄 0 ∈ R 3𝑁 (the generalized position vector), supposed to be at rest. We first introduce a few elements of the Hertz theory in Section A.1. In Section A.2, we describe the way the resulting perturbation of the overlaps 𝛿 𝑖 𝑗 is handled and its implication on the set of contact forces. Then, we show how this perturbation is related to the configuration itself in Section A.3 and derive the wave equation in Section A.4.

A.1. Embedding Hertz theory

In Section 2.2.2, the normal force between two grains 𝑖 and 𝑗 is characterized, at the local scale, by a positive scalar 𝑓 𝑖 𝑗 𝑛 > 0. The Hertz theory of contact provides a useful framework to model an elastic force at a contact between two particles. The general expression of the Hertz normal force between two grains is

𝑖 𝑗 , (A 1)


Subterráneas Hechas

Por El

Hombre En

Macizos Rocosos

L O Suárez-Burgoa

Palabras Clave

A SCHEME FOR CLASSIFYING THE SIZE OF MAN-MADE UNDERGROUND CHAMBERS AT ROCK MASS UN ESQUEMA PARA CLASIFICAR EL TAMAÑO DE CÁMARAS

Keywords: Size classification, chambers, rock mass Clasificación por tamaño, cámaras, macizo rocoso

No criteria for size classification of man-made underground chambers have been proposed until now. It could be pertinent to have a size classification scheme for characterizing underground spaces on rock mass, in order to avoid myths about the limits of these underground structures. In this paper, a simple volume-based size classification scheme for chambers is proposed by using a data base of 510 hydroelectrical projects located around the world that have their chambers inside the rock mass. The proposed classification has four categories.

INTRODUCTION

It is a common practice to apply the name chambers (also named rooms) to large manmade underground excavations, comparing them with other types of underground excavations (e.g., galleries, tunnels). Nevertheless, no criteria for size classification of these underground excavations have been proposed until now. Chambers can be differentiated with respect to other underground openings, because their three main dimensions (length, height and span, or width) are approximately the same. Galleries, shafts and tunnels are characterized by one of its dimensions being much larger than the other two. Another common characteristic of chambers is that they are ordinarily excavated by drilling and blasting, and this is typically done in phases because of their size. Other types of underground excavations, such as shafts, galleries and tunnels, can be achieved manually (for small galleries) or by mechanized methods (e.g. rise boring [for shafts] and Tunnel Boring Machine (TBM) [for tunnels]), without overlooking the usual method of drilling and blasting. It could be pertinent to have an undergroundsize classification scheme for characterizing the size of chambers at rock mass in a numerical unambiguous way; rather than only referring to particular, experienced but biased professional opinions, or particular positions of some group of persons following a particular interest. In countries that are in development, as is in South America, there was encountered many times common myths about the hazards that are involved when it is proposed some "huge" chamber; when in reality one is proposing a regular in-size chamber. The primary objective of this classification would be to ease the technical and scientific communication between people dealing with these kinds of underground openings. In this technical communication a simple volume-size classification scheme is proposed, to be applied to underground chambers, rooms or halls at rock mass. The proposed classification has just four categories, according to the volume of the opening.

CLASSIFICATION CRITERION

The classification criterion is based on the approximate volume of the opening that sustains the rock mass. This volume is calculated by multiplying the three main dimensions of the chamber: length (L), height (H) and span (S), without reference to details of the shape (e.g., a chamber dome).

3.

USED DATA The used data comes from the article of [START_REF] Hönish | The world's underground hydro power plants in[END_REF] and provided electronically by the same author (Hönish, 2007). This data base has 510 items, and fields recovered for each chamber were:

 chamber project name and country;  chamber length (L), height (H) and span (S);  real transversal area perpendicular to the length dimension (incomplete and available only in some chambers);  real volume (incomplete and available only in some chambers). The first four columns of Table 1 show ten of these data for randomly selected projects, as an example. The ten data were chosen randomly from the complete data base. First four columns were taken from [START_REF] Hönish | The world's underground hydro power plants in[END_REF]Hönish ( , 2007) ) and last four columns here calculated. Not available: na.

4.

DATA PROCESSING From the selected data fields it was calculated the following variables:

 rough area (A), determined by multiplying the height and the span A= H S;  the dimension ratio (R), which is the ratio of the length with the square root of the rough area (i.e. R= L A -0,5 );  the rough volume (Vr) that is obtained by multiplying the three defined lengths L, H and S;

 natural logarithm of Vr, i.e. log(Vr). The last four columns of Table 1 show these calculated values for the chosen samples. From all these data, univariate statistics were obtained, which are summarized in Table 2. Also, Figure 1 shows the histogram of these data for the variables Vr and ln(Vr); and Figure 2 the histograms of some ratios of interest as R, L/H, L/S and H/S. In all cases, the number of bins for the histograms was determined as a value close to the recommended by the Sturges (1926) equation, which gave a value of 10 bins for the 510 data items (an exception was for variable Vr where one wanted to have short ranges of bins).

In Figure 1b one can look that the logarithms of the volumes of the chambers have a normal tendency.

Vr in 10 3 m 3 from 0 to 1 000 in 50 bins of 10 units log(Vr) from 1,5 to 7 in 12 bins of 0,5 units

BENCHMARKS FOR CLASSIFICATION

As a first approach in this classification scheme, some reference benchmarks were chosen: an upper, an intermediate and a lower benchmark so as to define classification group limits from actual data presented at the database of [START_REF] Hönish | The world's underground hydro power plants in[END_REF]. As the upper benchmark, the biggest natural chamber of the world was chosen: the Sarawak natural chamber in Lubang Nasib Bagus in Malaysia. This chamber is approximately 700 m long and 396 m high and has a span of 230 m, with a total estimated volume of 64 × 10 6 m 3 . For an intermediate benchmark, some big man-made chamber constructed and reported to date was chosen: a chamber for the electrical station at the Robert-Bourassa Hydroelectrical Project in Quebec in Canada. This chamber is 483 m long and 45 m high and has a span of 22 m, with an estimated volume of 478×10 3 m 3 . For defining the lower benchmark of this classification, some difficulties became apparent because chambers can be as small as an exploratory or research room in the rock mass. The criterion for defining the lower benchmark of the classification scheme is based on the dimensions of a small exploration gallery. A common cross-section area 4 m high with a span of 4 m gives a cross-section area of 16 m 2 . In defining the third dimension of this lower benchmark, the shapes and sizes of some existing chambers were used to define a dimension ratio (R), which is the length (longest dimension) divided by the crosssectional equivalent dimension (i.e. R=L D -1 ), where D is the square root of the crosssection area perpendicular to L. From the data analyzed, one obtained a most frequent value of the dimension ratio (R) equal to 3,27 (i.e. the median of all the data). In this way, the third dimension of our lower reference volume is equal to 13 m (i.e., L= D R), and the approximate volume is equal to 209 m 3 . An intermediate lower benchmark criterion may be proposed from the dimensions of some actual big TBM (e.g. diameter of 15 m (D) and 177 m 2 of cross-section area). If the same obtained R value is applied, the length of L is 43 m. The lower reference volume of the chamber under this second criterion is be taken into account.

6.

CLASSIFICATION LIMITS To define the classification limits of the proposed benchmarks, some mathematical manipulations were performed. Small exploration gallery 0,2 R ratio from 1,5 to 18 in 9 bins of 2,0 units L/H ratio from 1 to 13,5 in 9 bins of 1,5 units L/S ratio from 1 to 24 in 10 bins of 2,5 units H/S ratio from 0,5 to 3,5 in 13 bins of 0,25 units These somewhat unusual numbers need to be manipulated to result in practical limits between size categories. Starting from the natural logarithm of the volumes of the assumed benchmarks (10 3 m 3 ), and computing the arithmetic media between these values, the limits may be given as actual volume values. The calculations and results are given in Table 3. The limits of the size categories are given in Table 4. 4 permitted to propose limits to the chamber-size categories, and the resulting size-classification scheme is shown in Table 5. In Figure 3 8.

CLASSIFICATION SCHEME The values of Table

CONCLUSIONS

The specific objective of the classification scheme presented here is to classify underground openings as a function of their volume. The classes are exhaustive and mutually exclusive, and are divided into four main classes: small, regular, large and huge. The criteria are based simply on the volume, approximately estimated. If the exact volume is determined, this would prevail over the approximate volume. The subdivision rules are based on a considerable number of manmade chambers. The classification is practical and robust and gives an economical classification without error and without possible declassification. It is applied only to chambers, rooms, halls or similar openings that form a volume-space in the rock mass in which the dimension ratio, defined as the ratio between the length and the square root of the product of the span and height, is less than 30 (i.e., R < 30). Even tough, this classification scheme was suggested after analysis a civil engineering chamber database, this classification may be used to classify also underground mining openings.

Figure 1 :

 1 Figure 1: Histograms of ratios variables of interest

Figure 2 :

 2 Figure 2: Histograms of ratios of interest

  these classification limits are shown with 16 chambers as examples.

Figure 3 :

 3 Figure 3: Man-made size chamber classification scheme a: LNG storage caverns Daejeon, Korea; b: Almenara Río Riachón Hydroelectric Project, Colombia; c: Transformer cavern Fljótsdalur hydropower project, Jökulsa í Flótsdal, Iceland; d: Powerhouse cavern Río Riachón Hydroelectric Project, Colombia; e: Valve house cavern Koyna hydroelectric project, Maharashtra, India; f: Transformer cavern Helms pumped storage project, Sierra Nevada mountain, California, USA; g: Transformer cavern Pumped Storage Plant Goldisthal. Thuringia, Germany; h: Transformer cavern San Carlos hydroelectric central, Antioquia, Colombia; i: Powerhouse cavern Hydroelectric system Miel 1, Colombia; j: Powerhouse cavern Porce III, Antioquia, Colombia; k: Powerhouse cavern San Carlos hydroelectric central, Antioquia, Colombia; l: Machine cavern Pumped Storage Plant Goldisthal. Thuringia, Germany; m: Powerhouse cavern Churchill Falls, Labrador, Canada; n: Desilting chambers Nathpa Jhaki hydroelectric power project; o: Powerhouse cavern Xiaolangdi hydro-project, Yellow River (Huanghe), China; p: Powerhouse cavern Robert-Bourassa generatin station, Québec, Canada; q: Cavern Xiluodu hydro-power station, Three Gorges Project, (Yangtze & Jinsha River) Sichuan and Yunnan, China.

Table 1 :

 1 Ten chambers of hydropower projects Project Name and Country L × H × S

		in m	Aa in m 2	Va in 10 3 m 3	A= H x S in m 2	R =L/ D Vr in 10 3 m 3	ln(Vr)
	Jajce 2, Bosnia	54 × 34 × 17	555	na	578	2,2	31	3,44
	Hinterm., Austria	59 × 24 × 16	290	20	384	3,0	23	3,12
	Bear Sw., U.S.A.	69 × 46 × 24	na	na	1104	2,1	76	4,33
	Acri 1, Italy	78 × 23 × 25	540	na	575	3,3	45	3,80
	Kilforsen, Sweden	97 × 33 × 18	485	48	594	4,0	58	4,05
	Waldeck 2, Germany	105 × 50 × 34 1390	106	1700	2,5	179	5,18
	Montezic, France	144 × 41 × 25 740	86	1025	4,5	148	4,99
	Wehr, Germany	219 × 33 × 20 570	135	660	8,5	145	4,97
	Bakun, Malaysia	300 × 52 × 32 1123	300	1664	7,4	499	6,21
	Vianden, Luxemb.	330 × 30 × 17 450	140	510	14,6	168	5,13
	Note:

Table 2 :

 2 Univariate statistics of the data

	Statistic Median Mean Std. Deviation Maximum Minimum	R 3,27 3,71 1,96 18,81 0,79	Vr in 10 3 m 3 55,49 84,31 96,42 985,60 5,18	Ln(Vr) 4,02 3,96 0,99 6,89 1,65

Table 3 :

 3 Benchmarks of the proposed classification

	Benchmark Upper Intermediate Robert-Bourassa Descriptor Sarawak natural chamber 64×10 3 Vr in 10 3 m 3 478 powerhouse chamber Lower 2 Big TBM section 7,7 Lower 1

Table 4 :

 4 Initial classification scheme limits

	Benchmark Upper Intermediate 478 Vr in 10 3 m 3 64×10 3 Lower 2 7,7 Lower 1 0,2	ln(Vr) 11,06 6,17 2,04 -1,56	Mean of ln(Vr) ---8,62 4,10 0,24	Vr in 10 3 m 3 ---5.5×10 3 61 1

Table 5 :

 5 Size limits of underground chamber classification

	Size classification Small Regular Large Huge	Volume limits (10 3 m 3 ) < 1 1 to 60 60 to 1 000 > 1 000

ACKNOWLEDGEMENTS

The author is most thankful to Dipl. Eng. K. Hönish for provide the author the electronic version of his data base.