Laurent Lacaze

Joris Bouteloup

Benjamin Fry

Edouard Izard

Immersed granular collapse: From viscous to free-fall unsteady granular flows

The collapse of a granular column in a liquid is investigated using numerical simulations. From previous experimental studies, it has been established that the dynamics of the collapse is mostly influenced by the Stokes number St, comparing grain inertia and viscous fluid dissipation, and the initial volume fraction of the granular column φ i . However the full characterisation of the collapse in the (St, φ i ) plane is still missing, restricting its modelling as a physical process for geophysical applications. Only numerical tools can allow providing a variation on the parameter space (St, φ i ) hardly reachable in experiments as well as a full description of the granular phase playing a major role on dense granular flows. For this purpose, a dedicated numerical model is used including a discrete element method to resolve the granular phase. The specific objectives of the paper are then two-folds: (1) the characterization of the dynamics of the collapse and its final deposit with respect to (St, φ i ) to complement available experimental data and (2) the description of the granular rheology according to these two dimensionless numbers including dilatancy effects. A simple predictive model stems from the obtained results, allowing to explain the evolution of the final deposit with (St, φ i ).

Introduction

Destructive natural phenomena such as snow avalanches, landslides, rock falls, and debris flows remain difficult to safeguard against. Given the complexity and the observed heterogeneity in dynamics of these gravity driven and unsteady multiphase-flows, it is natural to focus on canonical flows that can be controlled at the laboratory scale. For this purpose, several studies have been devoted to model configurations dealing with dry or wet granular flows.

Among other configurations, the collapse of an initial granular column in air and over a horizontal surface, referred to as dry granular collapse, reveals the behaviour of an unsteady granular flow starting from an initial unstable rest state and evolving towards a final deposit. This configuration thus shares many key features with natural situations, and it has naturally become the topic of several studies. It can be noted that a specific attention has been paid to the case of dry monodisperse-granular medium using laboratory experiments [START_REF] Lajeunesse | Spreading of a granular mass on a horizontal plane[END_REF][START_REF] Lube | Axisymmetric collapses of granular columns[END_REF][START_REF] Lajeunesse | Granular slumping on a horizontal surface[END_REF][START_REF] Lube | Collapses of twodimensional granular columns[END_REF][START_REF] Balmforth | Granular collapse in two dimensions[END_REF][START_REF] Lacaze | Planar collapse of a granular column: Experiments and discrete element simulations[END_REF] or numerical modelling at different scales [START_REF] Zenit | Computer simulations of the collapse of a granular column[END_REF][START_REF] Staron | Study of the collapse of granular columns using twodimensional discrete-grain simulation[END_REF][START_REF] Mangeney-Castelnau | On the use of saint venant equations to simulate the spreading of a granular mass[END_REF][START_REF] Staron | The spreading of a granular mass: role of grain properties and initial conditions[END_REF][START_REF] Lacaze | Axisymmetric granular collapse: a transient 3d flow test of viscoplasticity[END_REF][START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional navier-stokes model with a µ (i)-rheology[END_REF][START_REF] Girolami | A three-dimensional discretegrain model for the simulation of dam-break rectangular collapses: comparison between numerical results and experiments[END_REF][START_REF] Ionescu | Viscoplastic modeling of granular column collapse with pressure-dependent rheology[END_REF]. Even for such a canonical design, the dynamics of the collapse remains not fully understood, as many dry dense granular flows. Surprisingly, some simple features have nevertheless been observed in all these studies. In particular, the characterization of the dynamics and of the final deposit has shown that the collapse is mostly controlled by the aspect ratio a = H i /L i of the initial column with H i and L i its initial height and its initial horizontal length respectively.

Complexity has been increasingly added to the dry granular collapse to incorporate observable features in natural configurations, such as polydispersity or complex grainshape [START_REF] Phillips | Enhanced mobility of granular mixtures of fine and coarse particles[END_REF][START_REF] Degaetano | The influence of localised size reorganisation on short-duration bidispersed granular flows[END_REF][START_REF] Cabrera | Granular column collapse: Analysis of grain-size effects[END_REF], erodible bottom [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF][START_REF] Crosta | Numerical modeling of 2-d granular step collapse on erodible and nonerodible surface[END_REF] or the influence of the surrounding fluid [START_REF] Roche | On the run-out distance of geophysical gravitational flows: Insight from fluidized granular collapse experiments[END_REF][START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF][START_REF] Topin | Collapse dynamics and runout of dense granular materials in a fluid[END_REF][START_REF] Jing | Dynamics and scaling laws of underwater granular collapse with varying aspect ratios[END_REF][START_REF] Bougouin | Granular collapse in a fluid: Different flow regimes for an initially dense-packing[END_REF][START_REF] Jing | Flow regimes and dynamic similarity of immersed granular collapse: A cfd-dem investigation[END_REF]). In the latter situation, it has been highlighted that beyond a, the initial volume fraction φ i [START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF] and Stokes number St [START_REF] Bougouin | Granular collapse in a fluid: Different flow regimes for an initially dense-packing[END_REF]) can also play a significant role on the dynamics of the collapse. In other words, for given granular material and aspect ratio, the surrounding fluid can affect the granular collapse because of its influence on falling grain inertia through viscous dissipation, i.e. the St number, but also due to the initial compaction of the granular column, i.e. φ i .

The specific contribution of both St and φ i on the immersed granular collapse now requires to be characterized. For this purpose, a full picture of the influence of (St, φ i ) needs to be provided, as it has been for now only considered for a limited set of dimensionless parameters. In particular, the role of φ i observed at small St [START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF], and presumed to be negligible at large St according to the results obtained in dry configuration, suggests a combined influence of (St, φ i ) which would both originate from fluid viscosity. Thus, the influence of φ i has to be clarified as a function of St. Moreover, modelling granular collapse from physical analysis is required for larger scale situations. This means in particular to extract pertinent continuous models for the granular phase including its rheological behaviour. The granular collapse has already been shown to be an attractive test case for rheological models in the case of dry configuration [START_REF] Lacaze | Axisymmetric granular collapse: a transient 3d flow test of viscoplasticity[END_REF]. It has to be extended to the case of immersed situations.

In order to achieve these objectives, fluid-particle properties have to be continuously varied, and fluid-particle stresses have to be known. These objectives then suffer experimental limitations, and can, for now, only be achieved using numerical simulations. Here, a Volume-Average-Navier-Stokes/Discrete-Element-Method (VANS/DEM) coupling approach is used. This approach is referred to as the mesoscale approach in the following as it allows to solve the fluid phase at a scale slightly larger than the grain scale. This provides a reasonable scale of description to model laboratory scale configurations, keeping the Lagrangian description of individual grains. It has therefore been often used when dealing with immersed granular flows such as, among others, the immersed simple shear flow [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF]) and sediment transport [START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF][START_REF] Charru | Sediment transport and bedforms: a numerical study of two-phase viscous shear flow[END_REF][START_REF] Pähtz | Universal friction law at granular solid-gas transition explains scaling of sediment transport load with excess fluid shear stress[END_REF]. Even if most of the dynamics of the system is resolved with this approach, the properties and dynamics of the fluid flow in between grains remains modelled, and thus requires closure terms. This approach has nevertheless proved to be relevant for the above mentioned configurations and, particularly, to provide the relevant mechanisms, allowing to improve our understanding of the physics of these systems. This means that most of the required subscale physics of the fluid phase is captured by these closure models. Accordingly, this remains a pertinent scale approach to reach the main objectives depicted previously.

The paper is organized as follows. The numerical VANS/DEM method used in the paper, the collapse set-up and the dimensionless parameters are presented in section 2. In order to discuss the reliability of the mesoscale model used in this paper, an alternative resolved numerical approach for the fluid phase, solving part of the subscale physics, has also been tested. Results, comparisons and limitations are discussed in a supplementary material available at [link to be added]. The influence of (St, φ i ) on the dynamics of the collapse and its final rest state is then characterized in section 3, with a specific attention on the prevailing role of φ i when decreasing St. In section 4, the rheology of the granular material is extracted and characterised as a function of (St, φ i ), in the spirit of the µ -I model for dense granular flows. Prior to concluding, the link between the proposed rheological model and the obtained final deposit for immersed granular collapse is discussed in terms of a simplified predictive model in section 5.

A mesoscale approach for granular collapse modelling

The VANS/DEM numerical method used in the following has been explained in details in [START_REF] Charru | Sediment transport and bedforms: a numerical study of two-phase viscous shear flow[END_REF]. The method is thus only briefly recalled for record. Nevertheless, we pay attention here to highlight the terms that require closure models and the strategy adopted accordingly (see also the supplementary material available at [link to be added] for a discussion on these models). The physical setup and the associated dimensionless parameters used are then given.

Granular phase: DEM

The dynamics of the granular phase is solved using a classical Discrete Element Method (DEM). The motion of each solid particle j, with j ∈ [1, N p ] (N p being the number of particles), submitted to gravity acceleration, solid contact force with other particles and hydrodynamics force induced by a surrounding fluid is obtained by integrating Newton's equations for linear and angular momentum of a solid sphere of mass m j ,

m j dv j dt = m j g + F h j + k =j F c kj ; m j d j 2 10 dω j dt = Γ h j + k =j Γ c kj , (2.1)
where v j and ω j correspond to the linear velocity and the angular velocity respectively. F h j and Γ h j are the hydrodynamic force and torque exerted on each particle j respectively and F c kj and Γ c kj are the solid contact force and torque, respectively, exerted by a particle k on j if they are in contact. Note that with the previous formulation, F h j then includes the buoyancy contribution. The model of the hydrodynamics force on each particle will be specified in the next section.

Solid contacts between particles are modelled using a soft sphere approach, i.e. by allowing a small overlap between particles to mimic the deformation of real grains. This overlap is then used to calculate the contact force between grains, using a linear springdashpot model. The tangential force is limited by a Coulomb threshold allowing sliding between grains. Details of the model can be found in [START_REF] Izard | Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid[END_REF]. We only recall that the solid contact is then parametrized by the coefficient of restitution e and the coefficient of friction µ p between the two particles in contact and the stiffness of the considered material k n , or equivalently the contact time t c . In the case of granular material, we impose t c d/g to ensure rigidity of the material for dry situation [START_REF] Baran | Velocity correlations in dense gravity-driven granular chute flow[END_REF]. Here, t c = 2. 10 -3 d/g, or equivalently the stiffness of the particle is k n = 2. 10 5 mg/d. Note that in such limit, the stiffness of the material no longer influences the dynamics of the granular material [START_REF] Baran | Velocity correlations in dense gravity-driven granular chute flow[END_REF]. Actually, we also have t c ρd 2 /η, with ρ and η are the density and the viscosity of the surrounding fluid, respectively. The latter constraint ensures that the contact time is smaller than the diffusive time scale in the vicinity of the moving particle.

Fluid phase: VANS model

The fluid phase is solved at a scale larger than the grain size, as sketched in figure 1. For this purpose, the fluid-phase equations to be resolved are derived from the Navier-Stokes equations spatially averaged over a spatial scale larger than the grain size [START_REF] Jackson | The dynamics of fluidized particles[END_REF]. In other words, the NS equations are averaged over a volume V f of fluid contained in a volume of reference V larger than a solid particle and similar to the mesh cell volume (see figure 1 for a sketch). The Volume Averaged Navier-Stokes (VANS) mass and momentum equations read [START_REF] Jackson | The dynamics of fluidized particles[END_REF])

∂ε ∂t + ∇. (ε u f ) = 0, (2.2) ε D u f Dt = εg + 1 ρ ∇ • S - 1 ρ n f p/f p , (2.3)
where ε = 1 -φ is the local fluid volume fraction (φ is the particle volume fraction), • f and • p denote the average operator over the fluid phase and particle phase within the volume V, respectively. D/Dt is a fluid material derivative and is defined accordingly with respect to the fluid velocity u f as D/Dt = ∂/∂t + u f • ∇. The fluid-particle interaction force averaged over the particles within V is denoted n f p/f p , with n the number of particles per unit volume. To finish with, S is an effective stress tensor for the fluid which has to be specified. Note that using this volume averaged formulation, different contributions emerge in the stress S including the average fluid stress tensor over V f as well as what is referred to as the traction term in [START_REF] Jackson | The dynamics of fluidized particles[END_REF], defined at the interface between the particles and the fluid. The fluid-particle interaction force n f p/f p is simply related to the hydrodynamics force on each grain F h j in the lagrangian formulation (2.1) as

n f p/f p = 1 V j∈V F h j .
(2.4)

In the VANS/DEM model, particle-fluid interaction is not resolved at the grain scale and F h j , and thus n f p/f p , have therefore to be modelled. According to [START_REF] Jackson | The dynamics of fluidized particles[END_REF], the fluid-particle interaction force can be split into a buoyancy contribution and a local interaction force. For sake of simplicity, we assume that the remaining contribution is only a drag force. We thus write (2.5) where F M D corresponds to a drag model force at the scale of the grain, superscript M stands for Modeled. We choose in the following to model this drag contribution as (see [START_REF] Richardson | The sedimentation of a suspension of uniform spheres under conditions of viscous flow[END_REF][START_REF] Maurin | A minimal coupled fluid-discrete element model for bedload transport[END_REF], for instance)

n f p/f p = φ∇ • S + nF M D ,
nF M D = 18η d 2 φ (1 -φ) -ξ 1 + 1 60 Re p ( v p -u f ) , (2.6)
with Re p = ρd| v p -u f |/η, and ξ is a constant whose value lies in the interval [1,3]. Note that such force model does not provide any torque on the particle that would be induced by the fluid, then Γ h j = 0 in (2.1) for this specific method. The only torque applied to each grain therefore comes from solid contact. This is quite a crude approximation, but it provides the most simple model leading to the expected dynamics of the collapse (see the supplementary material available at [link to be added] ).

Moreover, we assume that the deviatoric part of the stress can be simply written as a generalized viscous stress, leading to a total stress S of the form

S = -p f I + η ef f ∇ u + T ∇ u ≡ -p f I + η M ef f ∇ u + T ∇ u , (2.7)
where • stands for an average over the mixture and

u = ε u f + φ v p is its average velocity, with φ v p = 1 V j∈V V j v j .
The viscosity η ef f is an effective viscosity. This viscosity has to be modeled and will thus be denoted η ef f ≡ η M ef f . Note that in (2.7), mixture velocity has been chosen instead of the fluid velocity to model the viscousdeviatoric stress component. In the literature both choices can be found as long as they provide an actual deviatoric tensor for the viscous stress (see [START_REF] Jackson | The dynamics of fluidized particles[END_REF], for a review). Then, if the fluid velocity is chosen, its trace contribution has to be subtracted, as the fluid phase is not divergence free in the case of the fluid [START_REF] Baumgarten | A general fluid-sediment mixture model and constitutive theory validated in many flow regimes[END_REF]. Using the mixture velocity allows to satisfy straightforwardly the previous constraint as it is divergence free, and moreover it has been shown to appear quite naturally in the case of dilute Stokes flow [START_REF] Jackson | The dynamics of fluidized particles[END_REF][START_REF] Zhang | Momentum and energy equations for disperse two-phase flows and their closure for dilute suspensions[END_REF]. The latter approach has also been used to describe dense situations [START_REF] Guazzelli | Rheology of dense granular suspensions[END_REF]. This definition has thus been chosen here regarding the state of our knowledge. Moreover, the following model will be used for the effective viscosity,

η M ef f /η = 1 + 5 2 φ + 7.6φ 2 + ζφ 3 , (2.8)
where we recognize the Einstein viscosity at O(φ), the Batchelor viscosity for hard sphere at O(φ 2 ), and an extra O(φ 3 ) term to account for higher order correction in such dense configuration.

In most of the simulations discussed in this paper, ξ = 1 in (2.6) and η M ef f /η = 1 + 5 2 φ in (2.8). However, a discussion on the influence of these models is given in section 3.4 including ξ = 2, ξ = 3 and higher order viscosity terms. Moreover, a resolved numerical approach at the grain scale has been used to show the relevance of these different models and to obtain an estimation of ζ ≈ 16 in (2.8) for the flows studied in this work (see the supplementary material available at [Link to be added] ).

The VANS equations (2.2) and (2.3) are solved numerically to obtain u f and p f on a regular meshgrid ∆x V = ∆y V = ∆z V = 2d, where the cell volume matches the elementary volume V (see figure 1). Note that we take advantage of the incompressibility of the mixture (fluid+grains) to solve a divergent-free equation for the mixture phase instead of equation (2.2), allowing the use of standard numerical algorithm developed for incompressible flows. ε and v p are obtained by averaging DEM results over the fluid cell. For more details on the numerical algorithm, the reader can refer to [START_REF] Charru | Sediment transport and bedforms: a numerical study of two-phase viscous shear flow[END_REF].

Setup and dimensionless numbers

A typical sketch of the configuration considered in this study is shown in figure 2. The computational domain consists of a rectangular box (L x , L y , L z ) in (x, y, z) with (x, y) the main propagation plane, y being opposed to gravity (see figure 2), and z is the third direction out of plane. Boundary conditions for both the fluid phase and the granular phase are periodic in the z direction. For the fluid phase, a no slip condition is imposed at the walls located at x = 0 and y = 0, while a slip condition is imposed at the walls located at x = L x and y = L y , ∂ u f /∂n = 0, n being the normal to the wall. Grains of diameter d are glued on the bottom plane -on a square grid centred at y = 0-to prevent the granular material from rolling on the bottom.

At the left side of the domain, x = 0, a rectangular column of base (L i , L z ) is filled up to a height H i with N p spherical grains of mean diameter d and same density ρ p (see figure 2(a)). Only a small dispersion -±5% in diameter of uniform distribution-is imposed to the grain diameter to avoid crystal-like pattern in the medium while keeping a monodisperse behaviour, i.e. no segregation is observed for such small dispersion in grain diameter. The geometry of the system is unchanged for all simulations. In particular, the aspect ratio of the initial column is set to a = H i /L i = 0.5, i.e. H i = L i /2, as its influence on the collapse has already been reported in several studies. The dimensionless base length of the column is L i /d = 64. The dimensionless size of the computational domain is such as L x /d ≈ 192, L y /d ≈ 51, and L z /d = 8.

At t = 0, the column is released in a liquid of density ρ and viscosity η. For t > 0, the granular medium first collapses during a so-called slumping phase, and eventually stops on a finite time t = t f . The final deposit can be characterized by two lengths, say a final spreading length L f and a final maximum height H f (see figure 2(b)), or equivalently, L f and a deposit slope tan α which will be defined later on.

At the grain scale, the dynamics is controlled by several dimensionless parameters, including solid-solid interaction and fluid-solid interaction. For the solid-solid interaction, first, the dimensionless parameter characterising the rigidity is defined as κ = (τ i /t c ) 2 with τ i = d ρ p /∆ρgH i a characteristic time of rearrangement imposed by the granular pressure (see for instance MiDi 2004, here the characteristic pressure is the granulostatic one at the bottom of the column using the apparent weight of the granular material ∆ρ = ρ p -ρ f ). The value of this parameter in the present simulations is κ = 5. 10 5 ensuring the rigidity of the material as mentioned previously [START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF]. In this rigid limit, the significant parameters characterising the solid-solid interaction are then, the coefficient of restitution e and the coefficient of friction µ p . In the following, they are set to e = 0.87 and µ p = 0.25. For the fluid-solid interaction in such gravity driven flow, two dimensionless parameters can be built upon the fluid and particle properties, say a Stokes number St and a density ratio r. The density ratio r is usually defined as r = (ρ p /ρ) 1/2 and is kept constant in the following, r = 1.6 corresponding to glass into water as mostly used in laboratory experiments. Its influence on the collapse has been reported in [START_REF] Bougouin | Granular collapse in a fluid: Different flow regimes for an initially dense-packing[END_REF], and is beyond the scope of the present paper. St can be defined in different ways, and we choose here to follow [START_REF] Bougouin | Granular collapse in a fluid: Different flow regimes for an initially dense-packing[END_REF] as

St = 1 18 √ 2 (ρ p ∆ρgd 3 ) 1/2 η .
(2.9)

In order to vary St, the viscosity of the fluid η is varied over 5 decades to range in St ∈ [6. 10 -3 , 60]. This range of St allows to cover both the viscous regime and the freefall regime as defined in Courrech du [START_REF] Courrech Du Pont | Granular avalanches in fluids[END_REF] and [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF]. Note that the validity of the model (2.8) for η M ef f when increasing St in this range is probably questionable, as fluid inertia at the scale of the grains could become not negligible. However, the small-scale fluid-inertia contribution should be limited in a dense granular configuration, and mostly dominant close to the upper surface of the granular medium. This is therefore expected to decrease quickly in the granular medium and not to be a dominant effect in the present configuration. Fluid inertia is therefore only accounted for through the drag model (2.6).

At the scale of the initial column, beyond the initial aspect ratio a set constant in the present study, recall a = 0.5, the column is characterized by its initial volume fraction defined as

φ i = j=1:Np V j L i H i L z , (2.10)
with V j the volume of grain j. The volume fraction of the initial column φ i is varied by modifying the filling procedure. We cover the range φ i ∈ [0.57, 0.63]. Note that the total number of grains is N p ≈ 20 000 whose exact value depends on φ i . The set of dimensionless parameters used is recalled in Table 1.

3. On the influence of St and φ i on immersed granular collapses

Preliminary considerations and typical observations

Experimental observations of dry and immersed granular collapses are noticeably different, particularly, concerning the influence of φ i reported in both cases. Even if differences on the transient and the final profile for varying φ i was observed for the dry case in an almost similar configuration [START_REF] Daerr | Sensitivity of granular surface flows to preparation[END_REF], they remain relatively small. Then the influence of φ i was not discussed any longer in the literature for the case of dry granular collapse. On the other hand, when dealing with immersed granular flows, [START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF] reported the influence of φ i as one of the dominant effects on the granular pile evolution. Their experiments were performed for St = {0.035, 0.065} according to the definition (2.9), i.e. in a viscous regime at small St. These experimental observations can be recovered using the VANS/DEM approach as reported in figure 3. In particular, one shows the temporal evolution of the dimensionless height profiles of the granular material h(x/d)/d for different values of (St, φ i ). We focus here on St = {6. 10 -3 , 60} and φ i = {0.57, 0.63}. In a viscous dominated situation, St = 6. 10 -3 (figures 3(a)-(b)), φ i clearly influences both the dynamics and the final rest state of the collapse. In particular, the initiation of the collapse is clearly delayed for the dense initial situation φ i = 0.63 with an initiation of the collapse at the right upper corner of the initial column (figure 3(b)) while no delay is observed for the initial loose configuration for which the collapse is initiated at the right bottom corner (figure 3(a)). Moreover, the spreading length is significantly more important for the loose situation. These observations are in qualitative agreement with experimental observations at small St [START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF]. On the other hand, for a particle-inertia (free-fall) configuration, St = 60 (figures 3(c)-(d)), the influence of φ i is less obvious. For sure, the spreading length is not affected which explained that no influence of φ i was reported in experimental studies dealing with dry collapse. However, a small difference can be observed close to the summit of the deposit. This observation reflects results reported by [START_REF] Daerr | Sensitivity of granular surface flows to preparation[END_REF]. Numerical results reported in figure 3 highlight the relevance of the mesoscale approach to provide the main behaviours obtained in previous experiments.

Granular morphology during collapse: a simple geometrical model

In order to highlight similitudes and differences of collapses in the range of (St, φ i ) considered here, the main oder morphological properties of the granular material during collapse is first considered. For this purpose, the trajectory of the centre of mass (2Y g -L i )/L i ,(2X g -L i )/L i during the collapse, i.e. for t ∈ [0, t f ], is shown in figure 4 for different values of (St, φ i ). Surprisingly, one observes that all trajectories remarkably collapse onto a single curve prior reaching the final deposit. The only difference between all the cases considered here is the position at which the trajectory of a given (St, φ i ) stops onto this universal curve (see in particular the inset of figure 4). Note that this observation confirms experimental results reported by [START_REF] Bougouin | Granular collapse in a fluid: Different flow regimes for an initially dense-packing[END_REF] for constant φ i ≈ 0.64. Obviously, a closer investigations indicates some small deviations, mostly at early stages of the collapse particularly for small φ i and small St (light grey squares). Assuming these deviations to be of second order, we focus here on the main curve holding all these trajectories.

The trend of the trajectory of the centre of mass can be predicted by deriving trajectories of simple geometric models. In particular, according to the shapes observed for small spreading length situations and larger ones, one considers that the shape of the collapse remains either trapezoidal or triangular during the entire collapse. We obtain for the trapezoidal shape and the triangular shape, respectively,

2Y g -H i H i ∼ - 2X g -L i L i 1/2 and 2Y g H i ∼ 2X g L i -1 . (3.1)
These two solutions are plotted in figure 4 with dashed line and solid line, respectively. Such simple solutions are shown to predict quite well the numerical data far from the initial state, i.e. for (2X g -L i /L i , 2Y g -H i /H i ) sufficiently far from (0, 0). In the latter case, trapezoidal and triangular predictive models give an identical trend. However, closer to the initial state, i.e. (2X g -L i /L i , 2Y g -H i /H i ) ≈ (0, 0), the trapezoidal model gives a better estimation of the collapse evolution in the range of parameters considered here (dash line). This is in accordance with the small a = 0.5 aspect ratio considered here, for which the granular medium remains in a trapezoidal shape during most of the collapse and particularly at early stages, even if the final rest state can be closer to the triangular shape for some cases (small φ i and small St for instance). Note that it has been shown that larger a situations are better estimated by the triangular model, as shown in [START_REF] Bougouin | Granular collapse in a fluid: Different flow regimes for an initially dense-packing[END_REF] for laboratory experiments.

Even if these simple models give a first insight of the dynamics of the collapse, the final rest state remains unpredicted. In order to close model (3.1), i.e. to predict the final state of the granular collapse, it would require to prescribe the time scale of the collapse and one of the final state morphological property, as the final spreading length or the final height or even deposit slope. This is discussed in the next section.

Final state morphology and time scale

In order to quantify the combined influence of (St, φ i ) on both the final state and the time scale of the collapse, figure 5 shows the final runout r = (L f -L i )/L i (a), the mean deposit slope tan α (b) and the dimensionless collapse time scale T 95 /T c (c), as a function of St and for several φ i (φ i = 0.57 green symbols, φ i = 0.59 red symbols, φ i = 0.63 blue symbols). In figure 5(b), the mean deposit slope tan α is defined as tan α = L f /H f for triangular deposit and tan α = (L f -L b )/H f for trapezoidal deposit (see inset of figure 5(b) for sketches and corresponding deposit lengths). Also shown in figure 5(b) is the maximum deposit slope tan α t obtained at the top of the deposit (dot symbols). In figure 5(c), T 95 corresponds to the time at which 95% of the final spreading length is reached by the front of the granular avalanche. Two different characteristic times T c have been considered here (full symbols and open symbols in the figure 5(c)) and defined as in [START_REF] Bougouin | Granular collapse in a fluid: Different flow regimes for an initially dense-packing[END_REF]. A so-called free-fall characteristic time

T c = T i = 2ρpHi
(ρp-ρ)g corresponds to the time that a particle needs to fall from a height H i with a constant acceleration induced by gravity g, assuming no interaction with the surrounding fluid (full symbols in figure 5(c)). A second characteristic time corresponds to a viscous time scale that a grain needs to fall from a height H i when its velocity remains constant as an equilibrium between weight and viscous drag. In this case, one can write

T c = T v = 18ηHi
(ρp-ρ)gd 2 and the corresponding results are shown by empty symbols in figure 5(c).

When St is large enough, St 10, figure 5(a) shows that (L f -L i )/L i is independent of φ i and reaches the expected dry situation (grey horizontal line from experiments of [START_REF] Lajeunesse | Granular slumping on a horizontal surface[END_REF] and black horizontal lines from dry DEM simulations; dry DEM simulations are preformed by removing the fluid solver on the same granular configuration). The spreading length then clearly decreases with decreasing St when the initial packing is dense enough, φ i = 0.63 and to a lesser extent φ i = 0.59, as could be expected from the role of viscous dissipation. Yet, an opposite trend is observed for the initial loose packing, φ i = 0.57. In particular, (L f -L i )/L i increases with decreasing St, reaching a maximum around St = 0.12. For St < 0.12, the spreading length slightly decreases with St as for the dense situations, but keeping a value larger than the large St limit one. For all φ i reported here, (L f -L i )/L i reaches a plateau when St 1, whose value increases for decreasing φ i . Accordingly, figure 5(b) shows a similar trend for the deposit slope tan α (squares). A noticeable difference is however observed at large St where a small difference on the deposit slope tan α (as it is defined here) is obtained. This remains in the range of values reported from experiments (horizontal grey line from [START_REF] Lajeunesse | Granular slumping on a horizontal surface[END_REF]) and dry DEM simulations (horizontal black lines). This difference at large St quantifies observations reported in figure 3(c-d), and is probably a signature of previous experimental observation on a similar configuration [START_REF] Daerr | Sensitivity of granular surface flows to preparation[END_REF]. Finally, the maximum slope close to the top of the deposit tan α t support previous observations (dots in figure 5(b)). Small differences between tan α t and tan α somehow measures an apparent uncertainty between macroscopic behaviour induced by (St, φ i ) on the shape of the final deposit and a local manifestation of different spatio-temporal dynamics, close to the front or close to the upper part of the deposit. Then, accordingly, the influence of (St, φ i ) is considered as remaining small for St 1.5, above which tan α t no longer depends on St or φ i .

At large St, the time scale of the collapse scales with a free-fall situation for which the fluid is disregarded (see full symbols in figure 5(c)). This time scale then strongly increases when St decreases (full symbols in figure 5 ef f /η = 1 + 5 2 φ) and St ≈ 6. 10 -3 (dark squares; blue squares correspond to the same simulations but shifting φi of ≈ 0.012, value based on the difference of rheological parameters between experiments and simulations, as explained in section 4.1), experimental data from [START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF] for St ≈ 5. 10 -2 (small dots) and for St ≈ 5. 10 -2 and a ≈ 0.5 or a ≈ 0.65 (big dots).

collapse, viscous vs. free-fall, changes for 1 < St < 10, actually close to St ≈ 1.5, as the transition observed for the deposit shape mentioned previously. This confirms the influence of the viscous dissipation on the final deposit. Yet, fluid viscosity can act on a very different manner on the dynamics and the deposit when St → 0, depending on the value of φ i . A non intuitive consequence is a possible enhancement of the spreading length due to the fluid viscosity. This means that the influence of the fluid viscosity is not only to slow the collapse down as could be expected through St whatever φ i , but plays another role through φ i at given St. The only other source of dissipation is obviously the granular friction, which has therefore to be strongly affected by φ i . The rheological properties of the granular material then plays a major role on these observations. This will be discussed in section 4.

An effective

St definition: unifying closure models and laboratory experiments. Before discussing the rheological properties, we focus on the generalization of the above mentioned results regarding the fluid phase closure models and with respect to experimental data available in the literature. So far, we have provided results for a given set of parameters of the closure models (2.6) and (2.8). According to the results discussed in the supplementary material available at [link to be added], this specific set of parameters should contain all the required physics of the fluid phase at the microscale, smaller than d, to provide the expected behaviour of the collapse, at least qualitatively. This would mean that the influence of a specific choice of the closure model, in the range of the one proposed in section 2.2 could only affect quantitatively the dynamics. We will show in this section that this quantitative influence of the closure models can actually be simply accounted for by defining an adequate effective Stokes number. Moreover, this will be then discussed in light of experimental results to highlight their predictability from the present simulations.

We have shown previously that the collapse is strongly affected by the transition from the viscous regime towards the free-fall regime. This transition has been shown to be controlled by the St number according to many configuration involving fluidparticle interactions. One can thus anticipate that closure models, affecting viscous dissipation at the microscale, should modify the critical range of St characterizing the transition from the viscous regime towards the free-fall regime. However, the qualitative trend obtained previously by varying (St, φ i ) should be maintained. To highlight that assumption, different models have been considered for φ i = 0.63 and varying St. Results are reported in figure 6(a). The different symbols correspond to different closure models as (ξ = 1, η M ef f /η = 1+ 5 2 φ) (squares; circles correspond to the same closure models but with an extra lubrication force added into the DEM contact model as in [START_REF] Izard | Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid[END_REF]

, (ξ = 2, η M ef f /η = 1 + 5 2 φ + 7.6φ 2 + 16φ 3 ) (up triangle), (ξ = 3, η M ef f /η = 1 + 5 2 φ + 7.6φ 2 + 16φ 3 ) (down
* = Stη/η M ef f uses η M ef f /η = 1 + 5 2 φ + 7.6φ 2 + 16φ 3
, even if the first order Einstein viscosity was used here, as lubrication adds dissipation which is not incorporated in the fluid phase solver. Then, one shows that models used for drag and/or viscosity, and even lubrication, only shift transition from viscous to free-fall regimes that could be easily understood using an appropriate definition of the St number. This is a remarkable results as it provides the relevance of the mesoscale approach for modelling immersed granular flows, even if the accurate closure models to be used are still debated.

A comparison of the results with available experimental data is given in figure 6. More particularly, figure 6(a) compares the influence of St in a dense configuration φ i ≈ 0.64 on the spreading length as experimentally studied by [START_REF] Bougouin | Granular collapse in a fluid: Different flow regimes for an initially dense-packing[END_REF], and figure 6(b) shows the comparison of the deposit angle as a function of φ i at small St as reported by [START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF]. In figure 6(a), the final spreading length is shown relatively to the dry one. This is done to ensure comparability between the numerical simulations and the experiments for which the bottom surface is different. In figure 6(b), the comparison is done on the deposit slope tan α. Moreover, note that, here, the values of St between numerical simulations and experiments are slightly different. However, their small values ensure a viscous regime, for which the final state remains roughly unchanged up to St ≈ 10 -1 . In both cases (figure 6(a) and figure 6(b)), a good agreement is obtained between the numerical simulations and the experimental data. In particular, the evolution of the final deposit as a function of the (St, φ i ) is clearly captured and, moreover, the quantitative ranges of evolution are in reasonable agreements. Note that, as the value of φ i is very sensitive, we also provide in figure 6(b) results for which φ i is slightly shifted for the numerical data (blue square; the shift value is not arbitrary and is based on the difference of equilibrium rheological states obtained in experiments and simulations as explained in section 4).

It can be noted that adding φ i in St * only accounts for the rate of fluid dissipation induced by compaction in the granular pores. This allows to capture the relevant range of variation of the final state as a function of fluid dissipation (figure 6(a)) but not to provide, solely, an explanation for the influence of φ i at given St reported in the previous section. Actually, the latter highlights a too significant influence on the final state to be only attributed to the definition of St * . This shows again that the obtained results are to be linked to a significant influence of the rheological properties of the granular phase.

Rheological model of the granular phase

Following [START_REF] Lacaze | Axisymmetric granular collapse: a transient 3d flow test of viscoplasticity[END_REF], we use here a coarse-graining approach to extract the equivalent stress tensor of the granular medium in this unsteady configuration. We do not recall here this methodology which has become quite standard and refer the reader to [START_REF] Goldhirsch | On the microscopic foundations of elasticity[END_REF] for details. It should be noted that the granular collapse is an interesting configuration to test and validate rheological models -as equilibrium models obtained from steady and 1D shear flows-in the case of unsteady and multidirectional shear configurations. Moreover, it allows to extract out-of-equilibrium local behaviours, that could be significant for unsteady configurations and required to improve models. In return, the averaging procedure of the coarse graining method has to be localised in space and time to extract the local component of the shear and stress tensors, which are space and time dependent [START_REF] Lacaze | Axisymmetric granular collapse: a transient 3d flow test of viscoplasticity[END_REF]. This usually leads to relatively dispersed results, as it will be shown in the following. Note, however, that the granular collapse allows to cover a wide range of dynamical properties such as shear rate and stress contribution, allowing to extract most of the rheological law, during a single event.

Here, coarse graining is performed on a regular grid in the (x, y) plane of resolution 2d at different times of the collapse. The averaging procedure is performed over volumes of gaussian shape in the (x, y) of standard deviation d and invariant in the z direction. From the coarse-grained results, the rheological law will be characterized by the local volume fraction φ, which can also be refedred to the local state of the granular material, and the coefficient of local effective friction µ defined as µ = τ p / p p , (4.1)

with both the deviatoric contribution of the granular stresse tensor τ p and the granular pressure p p being obtained from the coarse graining method. Once again, in such configuration all these quantities depend on (x, y) and t. In the case of a 3D flow, . refers to the second invariant of a tensor. The µ(I)-rheology as defined in [START_REF] Jop | A constitutive law for dense granular flows[END_REF] is used as the relevant rheological model, but including the extension of [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF] defined when transition from viscous to free-fall granular flow occurs, as observed for instance when increasing St (see previous section). Note that in gravity driven configuration, a distinction is sometimes made between inertial and free-fall regimes, for which inertial referred to as the state of the fluid drag at large Re p . However, [START_REF] Bougouin | Granular collapse in a fluid: Different flow regimes for an initially dense-packing[END_REF] shows that a fluidinertial regime that would differ from the free-fall regime is hardly observable in the range of properties of the immersed granular collapse covered here. Then, we do not clearly distinguish the inertial regime and free-fall regime as they both share quite common features for the granular phase. In [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF], the high St limit is referred to as an inertial regime in their neutrally buoyant configuration. However, inertial refers to the inertia of the grain and the rheology then shares the same features as for dry situation. This is basically the same as the free-fall case considered here, in which gravity now drives particle inertia through their apparent weight, i.e. weight and buoyancy.

The rheological model defines the effective friction µ of the granular material as function of a single dimensionless number, comparing a time scale of macroscopic deformation of the granular media induced by a shear γ p and a microscopic time scale of rearrangement of the media due to a confining granular pressure p p . Depending on the flow regime, [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF] has proposed two different relevant definitions for this dimensionless number, say J in the viscous regime and I in the free-fall regime. [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF] then proposed a combination of these different definitions I and J, labelled K, unifying a large range of St, from the viscous regime to the free-fall regime.

One therefore uses the dimensionless number K which can be written as

K = J + βI 2 , with J = η γ p p p I = γ p d p p /ρ p , (4.2)
where β is a constant, found to be ≈ 0.65 in [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF]. It should be noted that this value of β was found for 2D simulations. In the present case of 3D simulations, a different value could be obtained. However, in the range of parameters considered here, small variations of β around this reference value does not show any significant improvement for universal collapse of the rheological data. It has been chosen to be set to the value already proposed by [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF]. Note that dedicated 3D simulations should help in adjusting more accurately this parameter for 3D configurations.

4.1. Influence of St in the loose packing configuration In this section, we focus on the initial loose packing configuration, i.e. φ i = 0.57. This initial case is first considered as most of the initial column contributes to the granular motion, while for large φ i , most of the bottom-left corner remains static during the collapse. It will be shown that this configuration more clearly highlights an equilibrium state according to the rheological and state parameters (φ, µ), as obtained in steady flows.

The discussion on the rheological results obtained from the coarsening method is first discussed over the entire collapse without making distinction of specific regions, i.e. over the domain A ∈ {x > 2d, y > 2d, y < h(x, t)}, with h the time-dependent height profile of the granular medium, and over a linear time stepping scale, as t ∈ [0 T 95 ] with a time step T 95 /10. Note that T 95 is a function of St, then this time stepping account for the different dynamics of the collapse to cover regularly the full range of evolution. Moreover, boundaries are excluded from the analysis as they would need a specific attention such as an adapted coarse-graining approach [START_REF] Weinhart | Closure relations for shallow granular flows from particle simulations[END_REF].

Only the influence of St is thus considered here. The results obtained from the coarse graining approach are shown for (St, φ i ) = (6. 10 -3 , 0.57) (light green symbols) and (St, φ i ) = (6, 0.57) (dark green symbols) in figure 7. Note that larger values of St have been performed up to St = 60, for which no significant difference is observed with St = 6. St = 6 is thus discussed here without altering the analysis. The volume fraction φ and the effective friction µ are plotted as a function of K in figure 7(a) and figure 7(b) respectively. We can first note that, beyond an important dispersion of the results as anticipated, φ and µ show a clear evolution with respect to K, thus highlighting the relevance of the parameter K as the controlling dimensionless parameter to describe the rheology of immersed granular material. Moreover, the two different cases (St, φ i ) = (6. 10 -3 , 0.57) and (St, φ i ) = (6, 0.57) nearly collapse on a similar trend curve, even if a slight difference is actually observes as will be discussed in the following. These observations reinforce the role of K whose purpose is to unify the different regimes from viscous to free-fall.

The results obtained here show that at small φ i , the rheology is not significantly affected by St, and that the obtained critical volume fraction φ c = φ(K → 0) is larger than the initial state one. These observations suggest that the rheological state is close to a somehow universal one, independent of the initial state, which should therefore be expected to be the one obtained for steady state systems. We thus compare these results to the model proposed by [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF] obtained for steady state and considered here as the equilibrium state rheology. This model can be written as

φ eq = φ c -a √ K, µ eq = µ c + ∆µ K 0 /K + 1 , (4.3)
where subscript eq stands for equilibrium. In figure 7, black solid lines correspond to the model (4.3) with φ c = 0.592, a = 0.18, µ c = 0.24, ∆µ = 0.6 and √ K 0 = 0.2. These values correspond to the best fit that can be obtained from the data shown in figure 7 for St = 6. 10 -3 which shows less dispersion. Moreover, the values for ∆µ and √ K 0 are consistant with results obtained for the dry collapse [START_REF] Lacaze | Axisymmetric granular collapse: a transient 3d flow test of viscoplasticity[END_REF]. This leads to values of the fitting data which are in reasonable agreement with what is obtained in the literature. It can also be noted at this point that usual value of φ c found in experiments are around 0.58 (see for instance [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF]. The ≈ 0.012 difference of φ c between experiments and simulations has thus been used as the shift value for φ i in figure 6(b) (blue symbols).

The agreement between numerical results and the equilibrium model in figure 7 confirms the reliability of this model in the case of an unsteady and 3D configuration for a large range of St, at least for φ i = 0.57. However, a closer inspection of the results suggest that the influence of St is not inexistent, and that the spatio-temporal domain of extraction of the rheological characteristics requires a specific attention as it highlights different states during collapse. This point is discussed in the following, and will also be considered attentively in the next section.

As mentioned above, the general rheological trend given in figure 7 is obtained at 10 regularly spaced time steps for t ∈ [0 T 95 ]. At large St, this includes both accelerating and decelerating stages of the collapse, while at small St it only includes decelerating stages. This is highlighted in figure 8(a) in which the average value < K > is plotted as a function of t/T 95 for St = 6. 10 -3 and St = 6, and with < K > being the spatial average of K over the non static area A D ∈ {x > 2d, y > h s (x), y < h(x, t)}, where h s is the height profile delimiting the granular region which remains static for t ∈ [0 T 95 ]. We now focus on the decelerating stage, referred to as the resting stage, which is characterized as the time interval of decreasing < K >. Note that according to results shown in figure 8(a), the case St = 6. 10 -3 is not significantly affected by this new procedure as previous results were already obtained in the resting stage. However, this is somehow different for St = 6, for which part of the coarse-grained results are removed with this new procedure. It explains the more important dispersion of results for St = 6 in figure 7. Using this new procedure, each St is investigated independently. Insets of figure 7 show the average value of φ and µ, respectively, obtained on A D during the resting stage, for the two values of St previously. The average is obtained here by binning the spatio-temporal rough data in log-scale compartments of K. This first confirms the K-trend of both φ and µ. Moreover, µ and φ merges at large K for the different St, typically K > 10 -3 here, i.e. at the beginning of the resting stage. This would naturally support the fact that the granular material has reached a universal equilibrium state. However, according to the definition of φ c and µ c in (4.3) and obtained as K → 0 during the resting stage, this refined investigation shows that they are both functions of St. The dependence of φ c and µ c with St is shown in figure 8(b).

The evolution of µ c (St) and φ c (St) shown in figure 8(b) highlights a rapid transition around St ≈ 1 delimiting two asymptotic behaviours at small St and large St, at which plateaus are obtained. In particular, for both µ c and φ c , values on these plateaus are found to be slightly larger for St 1 than St 1. Note that the value found at St 1 is very close to the one reported for dry collapses in [START_REF] Lacaze | Axisymmetric granular collapse: a transient 3d flow test of viscoplasticity[END_REF]. However, the physical reason for the sudden decrease of µ c at small St remains unclear, even if viscous/free-fall transitions are usual observations for several properties in different configurations (see for instance [START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF], for case of particle bouncing). It can be at least noted that this is associated with the state of the granular medium for K → 0 as φ c follows a similar trend, such as µ c decreases with φ c . This indicates an influence of the history of the collapse on the final state, through St, which leads to a small deviation from a universal equilibrium law.

For an unsteady configuration such as the collapse, the state of the granular system and its associated rheological law, described here by φ and µ during the resting stage, i.e. for decreasing K, are close to the equilibrium law obtained for steady flows. However, Notes: the dilatation effect for φi = 0.57 is hardly observable, its correction tan ψ has therefore been neglected for this case (green dots in (c)); the dilatation angle correction tan ψ is only observable when grain motion is significant enough, data below J ≈ 10 -4 , i.e. in the plastic region, have therefore been disregarded for φi = 0.63 (blue dots (c)). Solid lines in (a) and (c) correspond to model (4.5) with φeq = φc. Solid line in (b) corresponds to model (4.3)

the unique function of K with constant parameters in (4.3) as obtained in steady configurations is not clear in the state. Here, different states can exist for the same value of K. In the loose configuration considered here, it has been shown to be possibly modeled by including a St dependency in µ c and φ c , keeping the trend of the equilibrium model.

To finish with, note that the obtained evolution of µ c (St) is in line with the conclusions drawn previously from the final state shape of the initial loose packing case, i.e. the spreading length increases while µ c decreases for decreasing St (see figure 5(b)).

Influence of φ i at small St

The same procedure can be followed to identify the influence of the initial volume fraction φ i . As the influence of φ i is clearer at small St, we focus here on St = 6. 10 -3 . In this case, the granular flow remains in the viscous regime and then K = J. Figure 9 reports the coarse-graining results obtained over the temporal interval t ∈ [0 T 95 ] with a time step T 95 /10 and over the spatial domain A ∈ {x > 2d, y > 2d, y < h(x, t)}, as done in the previous section. More particularly, the volume fraction φ and the effective coefficient of friction µ are plotted as a function of J for (St, φ i ) = (6. 10 -3 , 0.57) (green symbols) and (St, φ i ) = (6. 10 -3 , 0.63) (blue symbols) in figure 9(a) and figure 9(b), respectively. Note again that the equilibrium law is recovered for φ i = 0.57 on this spatio-temporal domain as discussed in the previous section. Yet, as observed in figure 9, the case φ i = 0.63 does not show the same trend as for φ i = 0.57; the influence of the initial decompaction of the granular medium prior collapsing cannot be disregarded from the rheological point of view. This behaviour has been reported in [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF], specifying a dilatancy angle ψ having to be accounted for when defining the effective friction angle. Following their analysis, we define the dilatancy angle ψ as

tan ψ = ∇ • v p .
(4.4)

The rheological model (4.3) can then be modified to include the influence of the dilatancy angle ψ as discussed in [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF] and [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] tan ψ = b(φ oeq -φ eq ),

µ oeq -tan ψ = µ c + ∆µ K 0 /J + 1 , (4.5)
where subscript oeq stands for out of equilibrium. According to this model, µ -tan ψ is plotted as a function of J in figure 9(c), with the same symbols as previously. For φ = 0.63, only values above J ≈ 10 -4 are reported here. As will be shown later, this is actually similar to considering only the moving region of the granular medium after the initial stages of the collapse corresponding to the initial expansion of the granular material. Below J ≈ 10 -4 , the µ -tan ψ does not collapse on the same curve (not shown here). The reason is probably to be attributed to the connection with the plastic region which is not well captured by the µ(J) rheology (probably to be related to the behaviour of the µ(I) rheology at small I as discussed in [START_REF] Barker | Partial regularisation of the incompressible µ(I)rheology for granular flow[END_REF]. This very specific behaviour should be addressed in future works but is beyond the scope of the present paper, as it does not affect significantly the physics and the dynamics of the collapse. However, for J > 10 -4 , one obtains a collapse of µ -tan ψ on a single function of J for φ i = 0.57 and φ = 0.63. In order to predict simply the rheological behaviour shown in figure 9, we can first assume that φ eq ≈ φ c in (4.5) as it is shown to be roughly constant in the inset of figure 9(a), where tan ψ is shown to be a single linear function of φ whatever J. This is moreover shown here to be equal to φ c = 0.592 (intersection of the solid line and tan ψ = 0 in the inset of figure 13(a)), similar to the value obtained with model (4.3) for (St, φ i ) = (6. 10 -3 , 0.57). Solid lines in figure 9(a,c) correspond to the model (4.5) with φ eq = φ c keeping the values of µ c , ∆µ and K 0 obtained previously and b = 5.5. This model prescribes relatively well the trends of rheological parameters obtained from the simulations, particularly for µ -tan ψ as a function of J and tan ψ as a function of φ, highlighting the influence of the deviation of the volume fraction φ with respect to a reference value φ c . Nevertheless, this model does not allow to provide the trend of φ oeq as a function of J, i.e. model (4.5) is not closed as we basically end up with one equation for two unknowns (φ oeq , µ oeq ).

A slightly more advanced approach to close the rheological model (4.5) in our configuration would be to prescribe φ oeq . This is not necessarily an obvious task without solving an appropriate time dependent equation. For sake of simplicity, we propose in the following a prediction based on the expected behaviour of φ oeq , or more likely on tan ψ, during the collapse. For that purpose, we rewrite the dilatancy angle ψ such as (4.6) where t = t/t m with t m = η/ p p is a microscopic time scale of rearrangement of the granular material in the viscous regime [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF]. As dilatation/compaction are precisely a rearrangement process, then we can assume that 1 φ Dφ D t = O(1) and then tan ψ ∼ J -1 . That said, we seek for a solution which is also regular when J → 0. φ(J → 0) can take several values as shown in figure 9(a). However, these values can be anticipated to be either φ i or φ c , depending on the situation. These situations can be distinguished as J → 0 at two different stages of the collapse, say at the very beginning of the collapse and when reaching its final rest state. In particular, the former case φ(J → 0) = φ i is more likely to be expected at the initial stage of the collapse, when the granular medium starts from an imposed volume fraction φ i . This will be referred to as the initiation stage. The latter one, φ(J → 0) = φ c , on the other hand, is more likely to be reached at the end of the collapse when the granular material returns to rest. As already discussed in the previous section, this resting stage could be considered as mainly described by its equilibrium state. This can be summarized as φ(J → 0, t → {0; t f }) = {φ i ; φ c }. A solution that would be consistent with (4.5) and (4.3) when J → 0 should then have the form tan ψ → b(φ i -φ c ) for an initiation state and tan ψ → b(φ c -φ c ) = 0 for a resting state. We thus propose the following simple model that would account for out-of-equilibrium state during the granular collapse

tan ψ = - 1 φ Dφ Dt γ-1 = - 1 φ Dφ D t J -1 ,
φ oeq = φ eq + ({φ i ; φ c } -φ c ) J 1 J + J 1 , µ oeq = µ eq + b ({φ i ; φ c } -φ c ) J 1 J + J 1 , (4.7)
where {φ i ; φ c } takes one of these values depending on the considered stage, i.e. initiation versus resting. Note again that with this description, the resting stage is nothing but an equilibrium state. J 1 in (4.7) is a regularization term allowing to reach a finite nonzero value when J → 0. Its value is found here to be J 1 ≈ 10 -3 to fit the numerical data. This can be seen as the transition from a dominance of dilatation/compaction effect towards a rheology controlled by the equilibrium state, this transition being found here to be controlled by the number J.

In order to show the relevance of the above mentioned model, one proposes in the following to extract and to separate results from simulations during the so-called initiation stage and resting stage, labelled in the following (.i) and (.ii) respectively. For that purpose, the distinction between these two stages has to be prescribed. We propose here to use the temporal evolution of averaged quantities over the granular medium. We thus plot the evolution of J and φ as a function of t/T v in figure 10(a) for both φ i = 0.57 (left column, labelled L. for loose) and φ i = 0.63 (right column, labelled D. for dense). Recall that T v is a viscous time scale. Here, the procedure to obtain • is the same as the one explained in section 4.1 but over the domain A Dt ∈ {x > 2d, y > h s (x, t), y < h(x, t)}, with h s delimiting flowing region and plastic region at each time step. Note that domain A Dt is slightly more restrictive than A D as used previously in section 4.1.

For the loose packing configuration L., one observes that the dynamics of the collapse, measured as J , quickly reaches its maximum inertial state, and then J → 0. At the same time φ increases during all the collapse (see figure 10 also corresponds to a change in slope of J (t/T v ). This time is chosen as a delimitation of region (L.i) and region (L.ii). For the dense packing configuration, the evolution of J and φ is slightly different (see figure 10(a) D.). In particular, a first initial stage (D.0) is observed during which J ≈ 0. This corresponds to the initial dilatancy during which both J and φ evolve very slowly. Then J increases suddenly. The rest of the dynamics can also be separated into two stages (D.i) and (D.ii) roughly similarly to the loose configuration, therefore corresponding to the initiation stage and resting stage.

Here, however, the delimitation is characterized by an inversion of evolution of φ (t/T v ) which decreases during (D.i) and increases during (D.ii). Note that in both cases (loose and dense), the chosen delimitation between stage (.i) and (.ii) more clearly corresponds to a change in the evolution of the state of the granular medium through the evolution of φ . However, it is clear from figure 10(a) that this separation also delineates dynamical stages of the collapse, as for instance J is maximum during the initiation stage (.i), while it slowly goes to zero during the resting stage (.ii). On the other hand, stages (.i) and (.ii) actually correspond to significant difference in the mass evolution, and thus the runout, between loose and dense configurations (as shown by the snapshot of the grains position in figure 10(a); black dots correspond to moving grains). In particular, T 95 is obtained at the end of stage (L.ii) for the loose configuration, while it is at the end of (D.i) or only early (D.ii) for the dense one. Note that the latter observation explains coarse-grained results shown in figure 9, as will be discussed later on.

Using the above-mentioned delimitation in time, the rheological and state variables µ and φ are plotted as a function of J during stage (.i) (red dots) and stage (.ii) (green dots) in figure 10(b). Here, coarse-graining is performed over the spatial domain A Dt . Even though coarse-grained results are a bit sparse, trends can be observed. These coarsegrained results are compared to the model (4.7) for φ(J → 0, t → 0) = φ i (solid line) and φ(J → 0, t → t f ) = φ c (dash lines). The qualitative agreement between the model and the simulations confirms the assumptions used to obtain model (4.7).

To finish with, we come back to results obtained in figure 9. Recall that there, coarsegrained results were obtained in the interval t ∈ [0, T 95 ] over a regular time grid in a lin-scale, and also includes the static region. According to results reported in figure 10 and the scale of T 95 for the loose and dense configuration, the spreading phase t ∈ [0, T 95 ] mostly lasts at the end of (L.i) and during (L.ii) for loose configuration while it mostly lasts during (D.i) for the dense one. In figure 11, one compares model (4.7) with this previously obtained coarse-grained results. Solutions of (4.7) are shown here for {φ i ; φ c } = φ i = [0.57 : 0.01 : 0.63] (thin lines), with highlights on {φ i ; φ c } = φ c (dash line) and {φ i ; φ c } = φ i = 0.63 (solid line). We can first conclude that the model proposed is in good agreement with numerical data showing its relevance for out-of-equilibrium situations. Moreover, the dense packing configuration highlights more clearly an out-ofequilibrium law model with {φ i ; φ c } = φ i while the loose packing situation shows an equilibrium state, i.e. {φ i ; φ c } = φ c , during most of the spreading phase t ∈ [0, T 95 ]. This latter remark does not mean that the loose configuration is not affected by the initiation stage, but that it happens on a short time scale compared to the entire collapse one, unlike the dense configuration.

Discussion

A simple phenomenological model is proposed in an attempt to provide a link between the rheological models obtained in the previous section and the morphology of the final state of the collapse. If the inertial acceleration of the granular slumping is not a dominant contribution of the collapse, which seems reasonable for immersed configurations and more particularly for small a, the final state is controlled by a balance between the macro pressure gradient, linked to the height gradient, and the friction term close to threshold. In other words, we suppose a quasi-static evolution of the collapse preventing the granular medium to spread further than its state at the threshold of motion. Obviously more complex situations could be imagined particularly when increasing a, but it is shown here that this assumption is sufficient to explain the influence of (St, φ i ) discussed so far. Assuming the final deposit to have a trapezoidal shape, one simply obtains

tan α = µ and L f -L i L i = a 2µ , (5.1)
with µ some effective friction coefficient at the macroscopic scale to be determined. Then the deposit only depends on the model prescribed to µ for a given a. According to the results obtained in the previous section, the friction parameter, even at threshold can vary for varying St and φ i . This leads to a finite interval of possible deposit slope, even for a quasi-static situation. Based on (4.6) and (4.7), we propose here to write this friction parameter as

µ = µ c (St) + b (φ i -φ c ) Γ -1 t φ , (5.2)
with t φ a time scale associated with the evolution of φ, from φ i to φ c during the initiation stage, and Γ -1 is a macroscopic time scale of deformation of the granular medium. Defining these two time scales is not straightforward for such a predictive model as different stages of the collapse have been shown to highlight different behaviours. A key point is thus to anticipate the stage controlling the final deposit. We thus base their definition on the observations extracted from figure 5 and figure 10. In particular, we have shown that φ i influences the final deposit mostly at small St. Moreover, φ i has been shown to affect the rheological behaviour during the initiation stage. We thus anticipate that the initiation stage is of main importance at small St while it probably does not affect much the final state at large St. Accordingly, Γ -1 /t φ should be order one at small St and decreases with increasing St.

According to the definition of t φ as given in the previous section, i.e. a microscopic time scale of rearrangement, we write t φ = η/∆ρgH i + βd ρ p /∆ρgH i , obtained for a granulostatic pressure over the column height, and accounting for a viscous to free-fall transition. Note that β is used here as in the definition of K for the sake of simplicity, as the transition parameter from viscous to free-fall regime. t φ then accounts for the difference in time scale for rearrangement depending on the state, viscous vs free-fall. According to the previous discussion, this time scale should be compared to a viscous time scale of macroscopic deformation of the granular medium, then Γ -1 /t φ would actually be order one at small St. This therefore implies that the relevant deformation accounting for dilatancy/compaction of the initial granular column, is associated with a viscous scale. It is thus written as Γ -1 = η/∆ρgH i , imposed by the balance between pressure gradient and viscous shear. This leads to

µ = µ c (St) + b (φ i -φ c ) 1 1 + βSt H i /d .
(5.3)

Figure 12 shows a comparison of this model for both the spreading length and final deposit angle with numerical results. Here, each line corresponds to the St dependency of the final deposit through the dilatancy/compaction term in (5.3), µ c being either the value obtained in figure 8(b) at large St (dashed lines) or at small St (solid lines). Even if the quantitative agreement seems poor, it clearly captures the trend of the evolution of the final state depending on (St, φ). Then the ingredient to understand and predict the final deposit seems to be captured. For sure, the specific evolution of the collapse and the different stages observed depending on (St, φ) has not been explicitly considered here, as for instance loose and dense configurations do not have a similar initiation stage and resting stage. A better description and understanding of the physics of the immersed granular collapse would thus require an intermediate description between the meso-scale approach and such a simple model. In particular, Euler-Euler simulations would be necessary to show the relevance of the local rheology obtained in the previous sections in the concept of continuum modelling. A full shallow-layer prediction on the unsteady flow could then be used to confirm the pertinent time scales required for a simple predictive model as discussed here.

Conclusion

Numerical simulations of an immersed granular collapse have been reported. The objective of the paper was to provide a characterization of the immersed granular collapse with respect to the Stokes number St and the volume fraction of the initial column φ i , which have been reported as the two parameters controlling the dynamics and the final deposit in experimental studies. The relative influence of these two parameters have been studied here using numerical simulations at a scale referred to as the mesoscale according to the scale of the fluid resolution, the granular phase being solved using a Lagrangian DEM approach. This method has been referred to as the VANS/DEM approach.

The influence of the two control parameters (St, φ i ) which characterize the immersed granular column in this case can thus be investigated. First, their influences on the dynamics and the final deposit have been considered. In particular, the influence of φ i is quite significant on the dynamics of the flow at small St, i.e. in the viscous regime as reported by [START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF], but tends to disappear at large St, probably explaining why its influence has been mostly not reported in dry granular flows. However, when St → 0, varying φ i can lead to very distinct behaviours, as it can enhance the mobility of the granular material compared to the dry case for small φ i , while the spreading length decreases with decreasing St for large φ i .

These simulations have been then used to provide visco-plastic rheological models proposed, based on models obtained for steady configurations in the literature. In particular, the well known µ(I) rheology has been considered as the rheological base model, but accounting for both an extension to the inertial number K as defined in [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF] unifying free-fall and viscous granular flows, and a dilatancy ψmodel. In particular, the rheology is characterised here by φ(K) and µ(K) obtained using a coarse-graining method. The rheological behaviour has been studied by separating two stages of the collapse: an initiation stage dominated by a dilatancy/compaction process and a resting stage assumed as being characterized by an equilibrium rheological law as in steady configuration. The initiation stage has been shown to be strongly influenced by φ i . Accordingly, an extension of the equilibrium rheological model accounting for dilatancy/compaction effect has been proposed, then referring to out-of-equilibrium rheological states. This shows to be pertinent for the configuration studied in this paper. However it would have to be confirmed in other situations, as the assumptions made here to obtain the model could be only relevant for the considered flow. The resting stage has been shown to slightly depends on St. However, this latter observation remains unclear as it has no clear physical support. This would deserve a specific attention, particularly on the process of compaction during the resting stage, probably influenced by St.

To finish with, a link between this rheology and the shape of the final deposit has been highlighted assuming a quasi-static evolution state towards an equilibrium between the pressure gradient and a friction term at threshold. A comparison with the numerical results shows that the influence of (St, φ i ) on the final state can be captured by this simple model. A more refined model would be required to improve the quantitative evolution of the final state in the (St, φ i ) parameter space.

Figure 1 .

 1 Figure 1. Sketch of the spatial scales used for the VANS/DEM approach in a 2D (x, y) cross-section. (∆xV , ∆yV ) define the grid size for the fluid phase resolution in the 2D plane. Grey disks are the 2D slice through the 3D spherical grain of diameter d in this 2D cross-section.

Figure 2 .

 2 Figure 2. Sketch of the 3D setup in a 2D streamwise (x, y) plane at a fixed z coordinate. Unit vector z is in the out of plane direction. (a) initial configuration at t = 0 and (b) final deposit for t t f .

Figure 3 .

 3 Figure 3. Dimensionless height profiles h(x/d)/d of the granular medium at different time: (a) (St, φi) = (6. 10 -3 , 0.57) and time step ∆t g/d = 150, (b) (St, φi) = (6. 10 -3 , 0.63) and time step ∆t g/d = 150, (c) (St, φi) = (60, 0.57) and time step ∆t g/d = 15 and (d) (St, φi) = (60, 0.63) and time step ∆t g/d = 15.

Figure 4 .

 4 Figure4. Trajectories of the dimensionless position of the centre of mass ((2Yg -Li)/Li,(2Xg -Li)/Li) (inset: log-log representation). Light grey symbols, dark grey symbols and black symbols correspond to φi = 0.57, φi = 0.59 and φi = 0.63, respectively (squares St = 6. 10 -2 and circles St = 6. 10 -1 ). Solid line (triangular shape) and dash line (trapezoidal shape) correspond to model (3.1).

Figure 5 .

 5 Figure 5. (a) dimensionless runout r = (L f -Li)/Li, (b) mean deposit slope tan α as defined in the text (squares) and maximum deposit slope tan αt (dots) as sketched in the inset and (c) dimensionless time scale T95/Tc, all as a function of St (green symbols φi = 0.57, red symbols φi = 0.59 and blue symbols φi = 0.63). Full symbols in (c) corresponds to Tc = Ti while open symbols are for Tc = Tv (see text for definitions of Ti and Tv). In (a) and (b), horizontal small black lines at large St correspond to dry simulations of the present configuration for φi = 0.57 and φi = 0.63, respectively, while dark grey lines are extracted from the experimental results of Lajeunesse et al. (2005).

Figure 6 .

 6 Figure 6. (a) Ratio of the spreading length L f -Li to the dry one L d f -Li as a function of the effective Stokes numbers St * = St × η/η M ef f (open symbols) or St * = St/ φi(1 -φi) -ξ (full symbols) for the initially dense configuration φi = 0.63 (see text for details): simulations for(ξ = 1, η M ef f /η = 1 + 5 2 φ) (squares)and an extra lubrication model (circles), (ξ = 2, η M ef f /η = 1+ 5 2 φ+7.6φ 2 +16φ 3 ) (up triangle) and (ξ = 3, η M ef f /η = 1+ 5 2 φ+7.6φ 2 +16φ 3 ) (down triangle); and experimental data from Bougouin & Lacaze (2018) (crosses). (b) mean deposit slope tan α as a function of φi: numerical simulations for (ξ = 1, η M ef f /η = 1 + 5 2 φ) and St ≈ 6. 10 -3 (dark squares; blue squares correspond to the same simulations but shifting φi of ≈ 0.012, value based on the difference of rheological parameters between experiments and simulations, as explained in section 4.1), experimental data from[START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF] for St ≈ 5. 10 -2 (small dots) and for St ≈ 5. 10 -2 and a ≈ 0.5 or a ≈ 0.65 (big dots).

  triangle). One shows in this figure that providing a relevant definition of St, say St * , allows almost to collapse the different cases on a master curve. In particular, this modification of St includes the influence of φ i through the effective models as St * = St×η/η M ef f (open symbols in figure 6(a)) or St * = St/ φ i (1 -φ i ) -ξ (full symbols in figure 6(a)), i.e. accounting for the φ i contribution included in either the effective viscosity closure model or the effective drag closure model, respectively. Note that for the specific situation of an additional lubrication force in the contact model (circles), the definition of St * is chosen as St

Figure 7 .

 7 Figure 7. Volume fraction φ (a) and effective coefficient of friction µ (b) (inset: average data on a regular K-scale from rough local data) as a function of K for (St, φi) = (6. 10 -3 , 0.57) (light green dots) and (St, φi) = (6, 0.57) (dark green dots). Solid lines correspond to model (4.3) with φc = 0.592, a = 0.18, µc = 0.24, ∆µ = 0.6 and √ K0 = 0.2. Dash Lines correspond to the model used at large St and extracted from figure 8(b), i.e. φc = 0.6 and µc = 0.29, keeping the other parameters equal to the former case.

Figure 8 .

 8 Figure 8. (a) Evolution of < K > as a function of time t/T95 over the domain AD (see text for details) for (St, φi) = (6. 10 -3 , 0.57) (black squares) and (St, φi) = (6, 0.57) (grey squares). (b) µc (black squares) and φc (grey dots) as a function of St. Vertical lines in (b) indicate the typical uncertainty in evaluating µc and φc.

Figure 9 .

 9 Figure9. Volume fraction φ (a) and effective coefficient of friction µ (b) as a function of J for (St, φi) = (6. 10 -3 , 0.57) (green dots) and (St, φi) = (6. 10 -3 , 0.63) (blue dots). Angle of dilatation tan ψ (4.4) as a function of J for (St, φi) = (6. 10 -3 , 0.63) (Inset of (a)). Corrected coefficient of friction µ -tan ψ as a function of J (c) (symbols are similar to figures (a) and (b)). Notes: the dilatation effect for φi = 0.57 is hardly observable, its correction tan ψ has therefore been neglected for this case (green dots in (c)); the dilatation angle correction tan ψ is only observable when grain motion is significant enough, data below J ≈ 10 -4 , i.e. in the plastic region, have therefore been disregarded for φi = 0.63 (blue dots (c)). Solid lines in (a) and (c) correspond to model (4.5) with φeq = φc. Solid line in (b) corresponds to model (4.3)

Figure 10 .

 10 Figure10. Time dependent rheology of the viscous collapse, St = 6. 10 -3 , for φi = 0.57 (loose L.: left column) and φi = 0.63 (dense D.: right column). (a) J (black symbols) and φ (grey symbols), averaged over domain ADt (see text for details), as a function of time (upper plots show snapshots of the granular medium within the different phases with static grains in grey and moving grains in black). (b) state and rheological parameters φ and µ as a function of J for the different phases as labelled in each figure. (.0), (.i) and (.ii) correspond to the different stages of the granular flow according to its dynamical behaviour and rheological state: (.0) is an initial decompaction stage prior collapsing only observation for dense configuration (D.), (.i) correspond to the initiation stage as an accelerating phase of the collapse and (.ii) is the resting stage during which the granular material goes to rest.

Figure 11 .

 11 Figure 11. Comparison of the coarse-graining results reported in figure 9 for φi = 0.63 (blue dots) and φi = 0.57 (green dots) with the model (4.7) (lines). (a) φ as a function of J and (b) µ as a function of J. Dash lines correspond to the equilibrium rheological model, i.e. φi = φc, while solid lines show the out-of-equlibrium model for φi = [0.57 : 0.01 : 0.63].

Figure 12 .

 12 Figure 12. Geometrical characteristics of the final deposit, as the deposit angle (a) and the dimensionless spreading length r = (L f -Li)/Li (b) as a function of St and different φi (green: φi = 0.57; red: φi = 0.59, blue: φi = 0.63). Symbols correspond to mesoscale VANS/DEM simulations (reported from 5(a-b), see the corresponding caption for details) and lines are model (5.3) with µc(St → 0) = 0.24 (solid lines) and µc(St 1) = 0.29 (dash lines).

Table 1 .

 1 Range of dimensionless parameters covered.

	a Li/d Lz/d	φi	r	St
	0.5 64	8	[0.57, 0.59, 0.63] 1.6	[6. 10 -3 : 60]
					(13 values evenly
					distributed on a log-scale)

† Email address for correspondence: laurent.lacaze@imft.fr

Acknowledgments

The authors thank Annaig Pedrono and Thomas Bonometti for their support in the development of the IBM/DEM and the VANS/DEM. We acknowledge the anonymous reviewers, whose thorough reviews helped to significantly improve the manuscript. Some of the computational time was provided by the Scientific Groupment CALMIP (project P1027), the contributions of which is greatly appreciated. This study has been supported by the 'Agence Nationale de la Recherche' in the frame of the project ANR MODSED No. ANR-12-JS09-0012.

Supplementary material

A supplementary study highlighting the relevance and the limitation of the models used in the present paper using a smaller scale of resolution for the fluid phase is available at ......