# Transient Rheology of Immersed Granular Materials

Zhuan Ge  Zhejiang University, 866 Yuhangtang Road, Hangzhou 310058, Zhejiang,
China
Key laboratory of Coastal Environment and Resources of Zhejiang Province
(KLaCER), School of Engineering, Westlake University, 18 Shilongshan Street,
Hangzhou, Zhejiang 310024, China.
Teng Man<EMAIL_ADDRESS>Key laboratory of Coastal Environment and
Resources of Zhejiang Province (KLaCER), School of Engineering, Westlake
University, 18 Shilongshan Street, Hangzhou, Zhejiang 310024, China.
Herbert E. Huppert  Institute of Theoretical Geophysics, King’s College,
University of Cambridge, King’s Parade, Cambridge CB2 1ST, United Kingdom
Sergio Andres Galindo-Torres<EMAIL_ADDRESS>Key laboratory of
Coastal Environment and Resources of Zhejiang Province (KLaCER), School of
Engineering, Westlake University, 18 Shilongshan Street, Hangzhou, Zhejiang
310024, China.

###### Abstract

In this letter, we investigate the transient rheological behavior of immersed
granular flows using both experiments of submerged granular column collapses
and corresponding numerical simulations. The simulations are performed with
the lattice-Boltzmann method (LBM) coupled with the discrete element method
(DEM) and provide a significant amount of data of the stress and deformation
conditions at different positions and times during the granular collapse. We
derive a new dimensionless number $\mathcal{G}$ that can unify the rheology of
transient granular flows in different regimes for all the simulation data
points. $\mathcal{G}$ smoothly transforms from an inertial number into a
viscous number, unifying both extremes of the rheology law. We also show the
need to introduce the kinetic stresses to achieve a universal relation. The
findings establish a transient constitutive framework for visco-inertial
granular flows, and are important for a better understanding of granular-fluid
mixtures in both natural and engineering situations.

††preprint: APS/123-QED

$Introduction-$Granular flows are ubiquitous in natural phenomena, such as
landslides, debris flows, and rock falls [1, 2, 3], and they can exhibit
different flow behaviors akin to solids, fluids, or gases [jaeger1996granular,
21, 22, 7]. Complex environmental conditions, transient fluctuations, and
highly dissipative interactions make it difficult to obtain a unified
constitutive law for their flow characteristics. Following pioneering works
[4, 5, 6] on dry granular flows in steady state conditions, it has been
determined that the apparent frictional coefficient $\mu$ can be considered a
sole function of the inertial number $I$. This inertial number
$I=\dot{\gamma}d/\sqrt{P/\rho_{s}}$ is defined as the ratio of a microscopic
time scale ($\sqrt{d^{2}\rho_{s}/P}$) to a macroscopic deformation time scale
($1/\dot{\gamma}$) [7, 5], where $\dot{\gamma}$ is the shear strain rate, $d$
is the averaged particle diameter, $P$ is the pressure applied to the granular
sample, and $\rho_{s}$ is the particle density. Lacaze et al.[8] verified this
theory through transient granular column collapse experiments, showing a
successful application of the $\mu(I)$ theory to granular flows.

At the other end of the spectrum, in the case of fully submerged granular
flows, Boyer et al.[9] showed the apparent friction $\mu(I_{v})$ of dense
suspensions is a function of the viscous number
$I_{v}=\eta_{f}\dot{\gamma}/P$, where $\eta_{f}$ is the fluid dynamic
viscosity. Trulsson et al.[10] further investigated the rheology of submerged
granular flow in the visco-inertial regime, and proposed a combined
dimensionless number $K=\lambda I^{2}+I_{v}$ for successfully describing the
submerged granular flows in different flow regimes (by varying the viscosity
of interstitial fluid). However, $\lambda$ is obtained by fitting and has no
clear physical definition. Furthermore, the work of Lacaze et al [11] showed
that $K$ cannot describe the transient rheology of immersed granular flows
[12] with enough accuracy. Although previous works improve the understanding
of granular flows in diverse conditions, there is still work to be done to
translate these models into predictive tools for natural hazards [13]. More
effort is needed to establish a universal constitutive law suitable for the
complex granular flow where both particle interactions and hydrodynamic forces
are non-negligible. In this case, where the granular assembly goes from a
granular skeleton (where friction rules) to a dense suspension (where the
viscosity of the fluid offers the greatest shear resistance), a proper
rheology law is still lacking.

Inspired by the complex dynamics presented in granular column collapses
reported in Refs.[thompson2007granular, 16, 3], we establish a numerical model
using the Lattice-Boltzmann method (LBM) coupled with the discrete element
method (DEM) to describe the immersed granular system and study the transient
granular rheology based on granular column collapses. The numerical model is
validated by immersed granular column collapse experiments. The one-to-one
comparison between experiments and numerical simulations gives significant
data for up-scaling the microscopic mechanism into macroscopic constitutive
behaviors. Then, the rheology of immersed granular flow is investigated for
systems in viscous, inertial, and free-fall regimes. The three flow regimes
[12], which depend on the square root of the grain/fluid density ratio
$r=\sqrt{\rho_{s}/\rho_{l}}$ and the Stokes number
$St=[\rho_{s}(\rho_{s}-\rho_{l})gd^{3}]^{1/2}/(18\sqrt{2}\eta_{f})$, are
classified to describe the effect of the fluid on grains in submerged granular
flows.

$Experimental\,setup-$As presented in Fig. 1, the dimension of the transparent
plastic tank is 38 cm$\times$6.5 cm$\times$20 cm. Three positions are
considered for the vertical retaining gate, corresponding to three different
initial column lengths $L_{i}=$ 3, 6, and 9 cm to generate different sizes of
initial granular columns. Plastic beads are used in this study. Their density
is 1.18 g/cm3, frictional coefficient is 0.34$\pm$0.01, and radius is
0.245$\pm$0.004 cm. Each test is recorded by a high-resolution camera with a
frame rate of 100 fps. The particles are immersed in water for which the
dynamic viscosity is 0.001 Pa${\cdot}$s and the density is 1 g/cm3.

Figure 1: (a) Experimental setup of the submerged granular column collapse.
(b) Dimensionless runout length ($\mathcal{L}=(L_{f}-L_{i})/L_{i}$) as a
function of the aspect ratio $\mathcal{A}=H_{i}/L_{i}$ of granular column
collapse in fluid for experiments and simulation results,
$\mathcal{A}_{c}=1.6$ is the transition point.

First, the retaining wall is placed at the desired position. Plastic particles
are then gently poured into the reservoir delimited by the wall to generate
the initial granular column, after which we pour the liquid into the tank
until it reaches the desired level. Once the fluid surface and particles are
static, we measure the initial length $L_{i}$ and initial height $H_{i}$ of
the granular column. Then, the retaining wall is removed suddenly, and the
column collapses and propagates into the tank. When particles stop
propagating, we measure the deposit length, $L_{f}$, which is the final front
position, and the final peak height, $H_{f}$. In this work, the initial aspect
ratio, $\mathcal{A}=H_{i}/L_{i}$, of the granular column is varied within the
range of 0.3-5. More details about the experiments are presented in the
Supplementary Material.

$Simulation\,setup-$We use DEM with frictional contact interactions modeled by
a Hookean contact law with energy dissipation [14]. LBM is used to simulate
the fluid flow in the pore space and to calculate the momentum exchange
between the fluid and the particles [15]. As shown in Fig. 1(b), the
normalized run-out distance $\mathcal{L}=(L_{f}-L_{i})/L_{i}$ shows good
agreement between experiments and numerical simulations. The transition point
appears when the aspect ratio is at $\mathcal{A}_{c}=$1.6, which is similar to
Ref.[16], and smaller than the granular collapse in dry conditions
[lube2005collapses]. During the granular collapse process, shown in Fig. 2,
the profile of the granular assembly also shows good agreement.

We implement immersed granular column collapse simulations at the initial
aspect ratio $\mathcal{A}$=1.73 (the initial height is 10.3 cm and the initial
length is 6 cm) with different viscosities as shown in Fig. 3 to investigate
their rheological behaviors in different flow regimes.

Figure 2: Submerged granular column collapse in water with a aspect ratio of
1.73 at 0.3s: (a) Experiment. (b) Numerical simulation. (c) Apparent
frictional coefficient ($\mu$) distribution at y = 3 cm, area $I$, $\mu\leq
0.4$, area $II$; $0.4<\mu\leq 0.6$, area $III$; $\mu\geq 0.6$. (d) Granular
temperature ($T$) distribution. Figure 3: Granular-fluid flow regimes in the
(St,r) plane according to Ref.[12] implemented in the simulation: dry
condition (black $\bigstar$ in the free-fall regime), immersed condition with
different fluid viscosity (pink $\bigtriangleup$, $\eta_{f}=$1 g/(cm$\cdot$s)
and dark blue $\rhd$, $\eta_{f}=$0.2 g/(cm$\cdot$s) in the viscous regime; the
green $\bigcirc$, $\eta_{f}=$0.01 g/(cm$\cdot$s) and red $\square$,
$\eta_{f}=$0.001 g/(cm$\cdot$s) in inertial regime.)

The granular system is then discretized into several representative volume
elements (RVE) as shown in Fig. 2(b) as a black grid, with side 1 cm, the
particles in the surrounding four cubics of each cell are used to obtain the
macroscopic information such as the averaged stress $\sigma$, strain rate
$\dot{\gamma}$, solid fraction $\phi$, and granular temperature $T=\delta
v^{2}/D$, where $\delta v$ is the velocity fluctuation, and $D$ is the space
dimension. The averaged stress is calculated by the contact term
$\sigma_{c}=\frac{1}{V}\sum_{p\in V}f_{i}l_{j}$, where $f_{i}$ is the $i$
component of the contact force between colliding DEM particles, and $l_{j}$ is
the $j$ component of the branch vector, and $i,j$ represents the $x,y,z$
direction. This contact stress tensor needs to be corrected to account for
small REVs as discussed in Ref.[17]. However, although this tensor is widely
used for DEM studies, it will be shown later how it must be complemented by
the kinetic stress tensor to achieve a universal rheology law. The pressure
$P$ and the shear stress $\tau$ are given by
$P=-(\sigma_{11}+\sigma_{22}+\sigma_{33})/3$ and
$\tau=\sqrt{1/2\tau_{ij}\tau_{ij}}$, respectively, where
$\tau_{ij}=\sigma_{ij}+P\delta_{ij}$ is the deviatoric stress tensor. The
equivalent strain rate tensor is calculated through the coarse-graining
approach as described in Ref.[18].

As shown in Fig. 2(c), the apparent frictional coefficient increases spatially
along the arrow direction. When the granular collapse is in the dense quasi-
static regime (area $I$), the apparent frictional coefficient is close to the
microscopic frictional coefficient, where $\mu\leq 0.4$. In the area $II$,
with increasing of granular velocity, the apparent frictional coefficient
increases, where $0.4<\mu\leq 0.6$. Larger frictional coefficient appear
associated with large granular temperatures near the interface between the
granular material and fluid (area $III$), where $\mu>0.6$. This is clear by
comparing area $III$ with the region of high temperature shown in Fig. 2(d).
We consider this as a hint to introduce the kinetic stress tensor in our
analysis as will be shown later. In Fig. 4(a), we plot the relationship
between $\mu$ and $I$, where we find that the $\mu-I$ rheology is sufficient
to describe the constitutive relationship of systems in inertial regimes [dry,
$\eta_{f}=0.001$ g/(cm$\cdot$s), $\eta_{f}=0.01$ g/(cm$\cdot$s)]. However, as
we increase the fluid viscosity, to reach the viscous regime, the $\mu-I$
relationship of systems with $\eta_{f}=0.2$ g/(cm$\cdot$s) and $\eta_{f}=1$
g/(cm$\cdot$s) deviates from the others. Plotting the relationship between
$\mu$ and $I_{v}$ in Fig. 4(b) shows that the data have a better collapse for
systems with $\eta_{f}=0.2$ g/(cm$\cdot$s) and $\eta_{f}=1$ g/(cm$\cdot$s),
but cannot capture the behavior of systems with $\eta_{f}=0.01$ and
$\eta_{f}=0.001$ and the dry sample. Systems in different regimes result in
distinct rheological behaviors, which further indicates that a universal
rheology is needed to describe the transient granular flow among free-fall,
inertial, and viscous regimes.

Figure 4: (a) Relationship between $\mu$ and $I$. (b) Relationship between
$\mu$ and $I_{v}$. Data are obtained from the beginning to end of granular
collapse in fluid with different viscosities.

The rheology of the granular flow is usually represented by the microscopic
particle movement time scale $t_{f}$ divided by the macroscopic rearrangement
time scale $\mathcal{T}=1/\dot{\gamma}$, where $t_{f}$ can be seen as the time
for a particle to travel over a characteristic length, e.g. the particle
diameter, $d$.

Figure 5: (a) $\mu$ as a function of $\mathcal{G}$ for different times in
fluids of different viscosities. The symbols are for the same simulation cases
with the same symbols coding in Figure 4. (b) Plots of the relative kinetic
effect, ($(P_{T}-P)/P_{T}$ in the normal direction, color of blue, and
$(\tau_{T}-\tau)/\tau_{T}$ in the tangential direction, color of red), with
respect to the inertial number, $I$, for both dry (symbol $\bigstar$) and
submerged (fluid viscosity is 1 g/(cm$\cdot$s), symbol $\bigtriangleup$)
cases. (c) Apparent frictional coefficient that includes the kinetic effect
$\mu_{T}$ as a function of $\mathcal{G}_{T}$ (incorporating kinetic effects).

In an ideal condition, the equilibrium of a single particle settling in the
fluid is given by $(\pi/6)\rho_{s}d^{3}\frac{du_{p}}{dt}=(\pi/4)Pd^{2}-F_{d}$,
where $F_{d}$ is the hydrodynamic force in submerged condition [23] and
$u_{p}$ is the particle velocity. In the granular flow, the Reynolds number is
usually very low, hence the hydrodynamic force can be assumed as the Stokes
force $F_{d}=3\pi\eta_{f}du_{p}$. During the settling process of a single
particle, the particle velocity increases until the hydrodynamic force is
equivalent to the inertial force, and the particle reaches the maximal final
velocity $u_{f}=Pd/(12\eta_{f})$. In previous works [6, 9], at the inertial
regime (or dry granular flow), the drag force is neglected, which assumes that
particles travel with a constant acceleration $a_{c}=3P/(2\rho_{s}d)$, with a
deduced settling time $t_{f}=t_{dry}=\sqrt{4\rho_{s}d^{2}/(3P)}$ and a time
scale ratio $t_{f}/\mathcal{T}=(2/\sqrt{3})I$. The constant factor
$2/\sqrt{3}$ is usually ignored. In the viscous regime, it is assumed that the
particle travels with the maximum velocity $u_{f}=Pd/(12\eta_{f})$ for the
characteristic length $d$, the settling time is $t_{f}=t_{sub}=12\eta_{f}/P$,
and the time scale ratio $t_{f}/\mathcal{T}=12I_{v}$. However, when the
particle flows in the fluid where the inertial force is comparable to the
hydrodynamic force, also known as the viscous-inertial regime, the rheology
can be described by neither the inertial number $I$ nor viscous number $I_{v}$
individually. Hence, we derive the travel time in a transient condition as
shown in the Appendix, where
$t_{f}=\frac{12\eta_{f}}{P}+\frac{\rho_{s}d^{2}}{18\eta_{f}}[1-e^{-\frac{36\eta_{f}}{\sqrt{3P\rho_{s}}d}}]$,
and obtain the ratio between microscopic time scale and macroscopic
rearrangement time scale as

$\begin{split}\mathcal{G}=\frac{t_{f}}{\mathcal{T}}&=12I_{v}+\frac{I^{2}}{18I_{v}}\left[1-e^{-\frac{36}{\sqrt{3}}\frac{I_{v}}{I}}\right]\\\
&=12I_{v}\left[1+\frac{\mathcal{ST}^{2}}{216}-\frac{\mathcal{ST}^{2}}{216}e^{-\frac{36}{\sqrt{3}\mathcal{ST}}}\right]\\\
&=\frac{2I}{\sqrt{3}}\left[\frac{6\sqrt{3}}{\mathcal{ST}}+\frac{\sqrt{3}\mathcal{ST}}{36}-\frac{\sqrt{3}\mathcal{ST}}{36}e^{-\frac{36}{\sqrt{3}\mathcal{ST}}}\right],\end{split}$
(1)

where $\mathcal{ST}=I/I_{v}$. We define
$\mathcal{B}=1+\mathcal{ST}^{2}/216[1-e^{-36/(\sqrt{3}\mathcal{ST})}]$ and
$\mathcal{C}=6\sqrt{3}/\mathcal{ST}+\sqrt{3}\mathcal{ST}/36[1-e^{-36/(\sqrt{3}\mathcal{ST})}]$,
so that $\mathcal{G}=\mathcal{B}(12I_{v})$ or
$\mathcal{G}=\mathcal{C}(2I/\sqrt{3})$. Increasing the Stokes number leads to
$\lim_{\mathcal{ST}\to\infty}\mathcal{G}=2I/\sqrt{3}=t_{dry}/\mathcal{T}$,
while decreasing $\mathcal{ST}$ leads to $\lim_{\mathcal{ST}\to
0}\mathcal{G}=12I_{v}=t_{sub}/\mathcal{T}$. Hence, as granular materials flow
from a free fall regime to a viscous regime, the dimensionless number
$\mathcal{G}$ naturally transforms from an inertial number to a viscous
number. Furthermore, the transition of flow regime of granular assemblies is
redefined as shown in Fig.6(a). For $\mathcal{ST}\leq 1$, the granular flow is
characterized by viscous number individually: this is the viscous regime. For
$\mathcal{ST}>100$, the inertial number should be used instead: this is the
inertial regime. For $1<\mathcal{ST}\leq 100$, the granular flow in the fluid
is characterized by both viscous number and inertial number: this is the
visco-inertial regime. However, the modified Stokes number
$St_{M}=I^{2}/I_{v}$ cannot obtain a universal flow regime transition under
different conditions, as shown in Fig.6(b). These analytical results are
consistent with the work of Trulsson et al.[10], where they present the
fraction of the power dissipated by each force field such as the inertial
forces and hydrodynamic forces with different $\mathcal{ST}$. The hydrodynamic
force is dominant when $\mathcal{ST}\leq 1$, while the contact force (inertial
force) is dominant when $\mathcal{ST}>100$. Owing to
$\mathcal{G}=12[I_{v}+I^{2}(1-e^{-36I_{v}/(\sqrt{3}I)})/(216I_{v})]$, the
factor in dimensionless number $K=I_{v}+\lambda I^{2}$ is
$\lambda=(1-e^{\frac{-36I_{v}}{\sqrt{3}I}})/(216I_{v})$, which depends on $I$
and $I_{v}$. This accounts for the reason why $\lambda$ varies in different
cases[10, 11, tapia2022viscous].

Figure 6: Evolution of $\mathcal{B}=\frac{\mathcal{G}}{12I_{v}}$ and
$\mathcal{C}=\frac{\mathcal{G}}{2I/\sqrt{3}}$ in terms of (a)
$\mathcal{ST}=I/I_{v}$, (b) $St_{M}=I^{2}/I_{v}$.

As shown in Fig. 5(a), the rheology of immersed granular column collapses in
different flow regimes can be partly described by the dimensionless number
$\mathcal{G}$. With an increasing inertial number, the apparent frictional
coefficient increases dramatically, as shown in Fig. 4(a), which results in
the divergence of the relation between $\mathcal{G}$ and $\mu$ in Fig. 5(a),
which is significant for the case of dry grains. We include the dynamic
effect, to account for this inertia, by introducing the kinetic stress tensor
$\sigma_{k}$. The total stress $\sigma_{T}=\sigma_{k}+\sigma$ is the
combination of the kinetic part $\sigma_{k}=\frac{1}{V}\sum_{p\in
V}m^{p}\delta v_{i}^{p}\delta v_{j}^{p}$ [24, 25] and the potential part
$\sigma$. To explain why introducing the total stress tensor works, we
introduce the following quantities: $\mathcal{G}_{T}$, and
$\mu_{T}=\tau_{T}/P_{T}$ are calculated from the total stress $\sigma_{T}$,
$\tau_{T}$ and $P_{T}$ are the shear stress and pressure derived from the
total stress. When the inertial number is large, the kinetic stress increases
the pressure (resulting in a positive $P_{T}-P$ difference), while reducing
the shear stress (negative $\tau_{T}-\tau$), as shown in Fig. 5(b). The
decrease in pressure $P$ and increase in shear stress $\tau$ results in an
ever increasing value for $\mu$. Once the total stress tensor is considered,
such discrepancy disappears for the total friction coefficient $\mu_{T}$.
Furthermore, as shown in Fig. 5(c), $\mu_{T}$ for both dry and submerged
conditions with a broad range of different viscosities could be generalized as
a function of $\mathcal{G}_{T}$, with the following functional form

$\mu_{T}(\mathcal{G}_{T})=\mu_{c}+\frac{\mu_{F}-\mu_{c}}{1+\mathcal{G}_{0}/\mathcal{G}_{T}},$
(2)

where $\mu_{c}=0.34$ is the minimum apparent frictional coefficient,
$\mu_{F}=0.7$ is the maximal macroscopic frictional coefficient, and
$\mathcal{G}_{0}=0.054$ is a fitting factor. In addition, Eq.2 can be used to
quantify transient granular flows and predict their natural phenomena for the
complex systems often observed in debris flows and submarine landslides.

$Conclusion-$Using DEM-LBM simulations, we have proposed a general
constitutive relationship that is suitable for granular flow in different
conditions. The new dimensionless number $\mathcal{G}$ can accurately describe
the granular material flow in fluid, where the effect of fluid and grain force
changes under different confining pressures, fluid viscosities, and
macroscopic deformations. $\mathcal{G}$ naturally transforms into $I$, when
hydrodynamic effects are negligable, and it converges to $I_{v}$ when those
effects are significant. It is shown how the kinetic stress, which is
essential for transient granular systems, must be introduced into the
rheological relationship in order to achieve universality. The proposed
rheology law could be used to formulate constitutive models for large scale
prediction at larger scales than the ones explored in this study.

This work is supported by the National Natural Science Foundation of China
(NSFC major project grant NO. 12172305). We thank Westlake University
Supercomputer Center for computational resources and related assistance. The
simulations were based on the MECHSYS open source library
(http://mechsys.nongnu.org).

## References

  * Hutter  _et al._ [1995] K. Hutter, T. Koch, C. Pluüss, and S. B. Savage, The dynamics of avalanches of granular materials from initiation to runout. part ii. experiments, Acta Mechanica 109, 127 (1995).
  * Tegzes  _et al._ [2002] P. Tegzes, T. Vicsek, and P. Schiffer, Avalanche dynamics in wet granular materials, Physical Review Letters 89, 094301 (2002).
  * Yang  _et al._ [2020] G. C. Yang, L. Jing, C. Y. Kwok, and Y. D. Sobral, Pore‐scale simulation of immersed granular collapse: Implications to submarine landslides, Journal of Geophysical Research: Earth Surface 125 (2020).
  * Iordanoff and Khonsari [2004] I. Iordanoff and M. Khonsari, Granular lubrication: toward an understanding of the transition between kinetic and quasi-fluid regime, J. Trib. 126, 137 (2004).
  * Da Cruz  _et al._ [2005] F. Da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, Rheophysics of dense granular materials: Discrete simulation of plane shear flows, Physical Review E 72, 021309 (2005).
  * Jop  _et al._ [2006] P. Jop, Y. Forterre, and O. Pouliquen, A constitutive law for dense granular flows, Nature 441, 727 (2006).
  * GDR and MiDi [2004] GDR and MiDi, On dense granular flows, European Physical Journal E (2004).
  * Lacaze and Kerswell [2009] L. Lacaze and R. R. Kerswell, Axisymmetric granular collapse: a transient 3d flow test of viscoplasticity, Physical Review Letters 102, 108305 (2009).
  * Boyer  _et al._ [2011] F. Boyer, É. Guazzelli, and O. Pouliquen, Unifying suspension and granular rheology, Physical review letters 107, 188301 (2011).
  * Trulsson  _et al._ [2012] M. Trulsson, B. Andreotti, and P. Claudin, Transition from the viscous to inertial regime in dense suspensions, Physical review letters 109, 118305 (2012).
  * Lacaze  _et al._ [2021] L. Lacaze, J. Bouteloup, B. Fry, and E. Izard, Immersed granular collapse: from viscous to free-fall unsteady granular flows, Journal of Fluid Mechanics 912 (2021).
  * Du Pont  _et al._ [2003] S. C. Du Pont, P. Gondret, B. Perrin, and M. Rabaud, Granular avalanches in fluids, Physical review letters 90, 044301 (2003).
  * Forterre and Pouliquen [2008] Y. Forterre and O. Pouliquen, Flows of dense granular media, Annu. Rev. Fluid Mech. 40, 1 (2008).
  * Cundall and Strack [1979] P. A. Cundall and O. D. Strack, A discrete numerical model for granular assemblies, geotechnique 29, 47 (1979).
  * Galindo-Torres [2013] S. Galindo-Torres, A coupled discrete element lattice boltzmann method for the simulation of fluid–solid interaction with particles of general shapes, Computer Methods in Applied Mechanics and Engineering 265, 107 (2013).
  * Bougouin and Lacaze [2018] A. Bougouin and L. Lacaze, Granular collapse in a fluid: Different flow regimes for an initially dense-packing, Physical Review Fluids 3, 064305 (2018).
  * Yan and Regueiro [2019] B. Yan and R. A. Regueiro, Definition and symmetry of averaged stress tensor in granular media and its 3d dem inspection under static and dynamic conditions, International Journal of Solids and Structures 161, 243 (2019).
  * Goldhirsch and Goldenberg [2002] I. Goldhirsch and C. Goldenberg, On the microscopic foundations of elasticity, The European Physical Journal E 9, 245 (2002).
  * Quartier  _et al._ [2000] L. Quartier, B. Andreotti, S. Douady, and A. Daerr, Dynamics of a grain on a sandpile model, Physical Review E 62, 8299 (2000).
  * Kim and Kamrin [2020] S. Kim and K. Kamrin, Power-law scaling in granular rheology across flow geometries, Physical Review Letters 125, 088002 (2020).
  * Roux and Combe [2002] J.-N. Roux and G. Combe, Quasistatic rheology and the origins of strain, Comptes Rendus Physique 3, 131 (2002).
  * Goldhirsch [2003] I. Goldhirsch, Rapid granular flows, Annual review of fluid mechanics 35, 267 (2003).
  * Cassar  _et al._ [2005] C. Cassar, M. Nicolas, and O. Pouliquen, Submarine granular flows down inclined planes, Physics of fluids 17, 103301 (2005).
  * Todd  _et al._ [1995] B. Todd, D. J. Evans, and P. J. Daivis, Pressure tensor for inhomogeneous fluids, Physical Review E 52, 1627 (1995).
  * Subramaniyan and Sun [2008] A. K. Subramaniyan and C. Sun, Continuum interpretation of virial stress in molecular simulations, International Journal of Solids and Structures 45, 4340 (2008).

##  I Appendix

The settling velocity of a particle in fluid is expressed as

$u(t)=u_{f}(1-e^{-at}),$ (3)

with $a=18\eta_{f}/(\rho_{s}d^{2})$. The traveling distance of the particle is

$s_{1}(t)=\int_{0}^{t}u(t)dt=u_{f}t-\frac{u_{f}}{a}(1-e^{-at})$ (4)

The relation between different traveling distances and the time is shown in
Fig. 7(a), $s_{3}=a_{c}t^{2}/2$ is for the particle settling in dry conditions
with the constant acceleration $a_{c}=3P/(2\rho_{s}d)$, $s_{2}=u_{f}t$ is for
the particle traveling with the final velocity $u_{f}=Pd/(12\eta_{f})$, which
is obtained in the condition of inertial force ($(\pi/4)Pd^{2}$) equaling to
the Stokes force ($F_{d}=3\pi\eta_{f}du_{f}$), and $s_{1}$ is for the particle
settling in fluid with the Stokes force $F_{d}=3\pi\eta_{f}du_{p}$. The
settling time calculated from different conditions shows that $t_{sub}$ and
$t_{dry}$ are always smaller than $t_{f}$. The settling time can also be
expressed as $t_{f}=t_{sub}+\Delta t$, where $\Delta t$ is a time gap between
the particle settling at terminal velocity and the one starting from rest.

Figure 7: Single particle settling in the fluid: (a) Evolution of the
traveling distance in terms of time; $s_{1}$, particle travels in fluid under
the hydrodynamic force; $s_{2}$ , particle traveling with the final velocity
$u_{f}=Pd/(12\eta_{f})$, where the inertial force equal to the hydrodynamic
force; $s_{3}$ particle traveling with a constant acceleration due to the
inertial force. The green and blue dot lines represent that with decreasing of
$\eta_{f}$, $s_{1}$ is close to $s_{3}$, $s_{2}$ is close to $t=0$. (b)
Variation of time difference between $s_{1}$ and $s_{2}$ in terms of traveling
distance $d/d_{0}$.

As shown in Fig. 7(a), increasing the traveling distance, $\Delta t$ increases
from 0 to a final time difference $\Delta t_{0}$. We assume that, when the
settling velocity $u=\alpha u_{f}$, where $\alpha\approx 1$, the time
difference is not changing anymore. In this case, the settling time
$t_{0}=-\ln(1-\alpha)/a$, and the traveling distance is
$d_{0}=-u_{f}/a[\ln(1-\alpha)+\alpha]$. Due to $t_{1}=d_{0}/u_{f}$, the final
time difference $\Delta t_{0}=t_{0}-t_{1}=\alpha/a$. Hence, the time
difference could be given by $\Delta t=\Delta t_{0}(1-e^{-bd/d_{0}})$. We can
obtain

$\begin{split}\frac{d}{d_{0}}&=-\frac{ad}{u_{f}[\ln(1-\alpha)+\alpha]}\\\
&=-\frac{18\eta_{f}d}{\rho_{s}d^{2}}\frac{12\eta_{f}}{Pd[\ln{(1-\alpha)}+\alpha]}\\\
\end{split}$ (5)

and the settling time of the particle in the fluid is

$t_{f}=\frac{12\eta_{f}}{P}+\frac{\rho_{s}d^{2}\alpha}{18\eta_{f}}\left[1-e^{\frac{216\eta_{f}^{2}}{P\rho_{s}d^{2}}\frac{b}{\ln(1-\alpha)+\alpha}}\right].$
(6)

Due to $\lim_{\eta_{f}\to 0}t_{f}=t_{dry}=\sqrt{4\rho_{s}d^{2}/(3P)}$, one can
obtain
$\frac{b}{\ln(1-\alpha)+\alpha}=-\frac{\sqrt{\rho_{s}P}d}{6\sqrt{3}\eta_{f}\alpha}$,
and Eq.6 is

$t_{f}=\frac{12\eta_{f}}{P}+\frac{\rho_{s}d^{2}\alpha}{18\eta_{f}}\left[1-e^{-\frac{36\eta_{f}}{\sqrt{3P\rho_{s}}d\alpha}}\right]$
(7)

Hence, the ratio of the microscopic particle time scale to the macroscopic
rearrangement time scale is given by

$\begin{split}\mathcal{G}&=\frac{t_{f}}{\mathcal{T}}=\frac{12\eta_{f}\dot{\gamma}}{P}+\frac{\alpha\rho_{s}d^{2}\dot{\gamma}}{18\eta_{f}}\left[1-e^{-\frac{36\eta_{f}}{\sqrt{3P\rho_{s}}d\alpha}}\right]\\\
&=\frac{12\dot{\gamma}\eta_{f}}{P}+\frac{\alpha
P}{18\eta_{f}\dot{\gamma}}\frac{\rho_{s}d^{2}\dot{\gamma}^{2}}{P}\left[1-e^{-\frac{36}{\sqrt{3}}\frac{\eta_{f}}{P}\frac{\sqrt{P}}{\sqrt{\rho_{s}d\alpha}}}\right]\\\
&=12I_{v}+\frac{\alpha
I^{2}}{18I_{v}}\left[1-e^{-\frac{36}{\alpha\sqrt{3}}\frac{I_{v}}{I}}\right],\end{split}$
(8)

since $\alpha$ is close to 1, and the dimensionless number $\mathcal{G}$ could
be expressed as

$\begin{split}\mathcal{G}&=12I_{v}+\frac{\alpha
I^{2}}{18I_{v}}\left[1-e^{-\frac{36}{\alpha\sqrt{3}}\frac{I_{v}}{I}}\right]\\\
&=I_{v}\left[12+\frac{\mathcal{ST}^{2}}{18}-\frac{\mathcal{ST}^{2}}{18}e^{-\frac{36}{\sqrt{3}\mathcal{ST}}}\right]\\\
&=\frac{I}{\mathcal{ST}}\left[12+\frac{\mathcal{ST}^{2}}{18}-\frac{\mathcal{ST}^{2}}{18}e^{-\frac{36}{\sqrt{3}\mathcal{ST}}}\right].\end{split}$
(9)

When the inertial force is dominant, the hydrodynamic force is nearly nil, one
can obtain

$\begin{split}\lim_{\mathcal{ST}\to+\infty}\mathcal{G}&=\frac{I}{\mathcal{ST}}\left[12+\frac{\mathcal{ST}^{2}}{18}-\frac{\mathcal{ST}^{2}}{18}e^{-\frac{36}{\sqrt{3}\mathcal{ST}}}\right]\\\
&=\frac{I}{18}\left[1-e^{-\frac{36}{\sqrt{3}\mathcal{ST}}}\right]\mathcal{ST}\\\
&\approx\frac{I}{18}\left[\frac{36}{\sqrt{3}\mathcal{ST}}\right]\mathcal{ST}\\\
&=\frac{2}{\sqrt{3}}I=\frac{t_{dry}}{\mathcal{T}},\end{split}$ (10)

when the hydrodynamic force is dominant, one can obtain

$\begin{split}\lim_{\mathcal{ST}\to+0}\mathcal{G}&=I_{v}\left[12+\frac{\mathcal{ST}^{2}}{18}-\frac{\mathcal{ST}^{2}}{18}e^{-\frac{36}{\sqrt{3}\mathcal{ST}}}\right]\\\
&=12I_{v}=\frac{t_{sub}}{\mathcal{T}}\end{split}$ (11)