Effect of viscosity on the avalanche dynamics and flow transition of wet granular matter

The dynamic behaviour of granular flows is important in geo-mechanics and industrial applications, yet poorly understood. We studied the effects of liquid viscosity and particle size on the dynamics of wet granular material flowing in a slowly rotating drum, in order to detect the transition from the avalanching to the continuous flow regime. A discrete element method (DEM) model, in which contact forces and cohesive forces were considered, was employed to simulate this flow behaviour. The model was validated experimentally, using glass beads in a wooden drum and water–glycerol mixtures to tune the liquid viscosity. The DEM simulations showed comparable results to the experiments in terms of average slope angle and avalanche amplitude. We observed that the avalanche amplitude, flow layer velocity and granular temperature decrease as the liquid viscosity increases. This effect is more pronounced for smaller sized particles. The increase in viscous forces causes the flowing particles to behave as a bulk, pushing the free surface towards a convex shape. In addition, avalanches become less pronounced and the granular flow transitions from the avalanching regime to the continuous regime. The avalanching flow regime is marked by intermittent rigid body movement of the particulate bed and near-zero drops in the granular temperature, while no rigid body movement of the bed occurs in the continuous flow regime. We identified the avalanching-continuous flow transition region as a function of a dimensionless granular Galileo number. © 2020 Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences. Published by Elsevier B.V. This is an open access article under the CC BY license (http:// p s F r T L T R 2 Introduction Granular flows have been extensively studied by researchers over the past years (Duran, 2012; Jarray, Magnanimo, Ramaioli, & Luding, 2017; Taghizadeh, Hashemabadi, Yazdani, & Akbari, 2018; Xu, Xu, Zhou, Du, & Hu, 2010). Understanding the dynamics of flowing granular matter is important in both industrial applications and the management of geo-physical phenomena, such as landslides and avalanches (Pudasaini & Hutter, 2007; Wu, 2015). A granular avalanche occurs when the slope of a pile of grains exceeds its maximum angle of stability m, whereafter the surface angle decreases until the angle of repose r is reached. Rotary drums filled with ∗ Corresponding authors. E-mail addresses: j.h.kasper@student.utwente.nl (J.H. Kasper), a.jarray@utwente.nl (A. Jarray). w

Avalanche dynamics are often characterized by the evolution of the slope angle and particle velocity. The avalanche amplitude Â = Â m − Â r is found to be a well-defined quantity, that increases with particle size (Balmforth & McElwaine, 2018;Jaeger, Liu, & Nagel, 1989). Courrech Du Pont et al. (2005) observed that the velocity profile during an avalanche decreases exponentially in pile depth direction (perpendicular to the pile surface) and that no https distance-to-contact point vector h measure of fluctuation in surface angle I mass moment of inertia of a particle k elastic spring constant m particle mass N number of bins N a number of adjacent particles N p number of particles in the drum n case number n unit normal vector q torque vector due to rolling resistance R drum radius r particle radius r * effective radius at contact between two particles T g granular temperature t a total avalanche duration t s sampling duration U fluctuation velocity vector u particle velocity vector V bond liquid bond volume v flow layer velocity v c characteristic velocity of a particle in the flow layer w drum width x, y (particle) coordinates x particle position vector steady state develops. Li et al. (2018) found that the variance of the particle velocity remains nearly zero when the particles are at rest, while it reaches a peak and gradually decreases as an avalanche occurs. The addition of interstitial liquid between the particles induces viscous and capillary cohesive forces and alters the dynamics of granular avalanches (Courrech du Pont, Gondret, Perrin, & Rabaud, 2003;Hornbaker, Albert, Albert, Barabási, & Schiffer, 1997), resulting in a qualitatively new behaviour. Several studies have demonstrated the importance of the liquids' properties on granular flow and pile stability (Chou, Liao, & Hsiau, 2010;Chou, Yang, & Hsiau, 2019;Courrech du Pont et al., 2003;Finger & Stannarius, 2007;Kanule, Ng'etich, & Rotich, 2019;Liao, 2018;Liu, Yang, & Yu, 2011;Louati, Oulahna, & De Ryck, 2017;Roy et al., 2019;Samadani & Kudrolli, 2001;Shi, Roy, Weinhart, Magnanimo, & Luding, 2020). Chou et al. (2019) examined the effect of low viscosity liquids ( < 10 mPa s) on creeping granular flows, and observed that the thickness of the flowing layer slightly increases with liquid viscosity. Samadani and Kudrolli (2001) showed that the interstitial liquid causes particles to clump together, yielding an increase in the angle of repose. Shi et al. (2020) demonstrated that liquid-induced cohesion can either decrease or increase the packing fraction, depending on the inter-particle friction. Liquid-induced cohesion has also been shown to promote the collective motion of particles, by keeping them in contact (Roy et al., 2019). Liu et al. (2011) showed that by increasing the liquid surface tension, the flow of wet particles transitions from a steady continuous surface flow to a discrete avalanching flow. Other studies by Tegzes et al. (1999), Tegzes, Vicsek, and Schiffer (2003) and Tegzes, Vicsek, and Schiffer (2002) investigated the transition from intermittent avalanches to continuous flow as a function of the liquid content, and identified three fundamental regimes of wet particles: the granular, the correlated, and the viscoplastic regime.
Despite these efforts, it is not yet fully understood how the interstitial liquid affects particle contacts, and in turn the flow of granular avalanches. A micromechanical interpretation is lacking, especially for highly viscous liquids ( > 100 mPa s). Whether the viscous forces enhance or reduce the volume fraction of wet particles also remains unanswered. Furthermore, the viscosity-induced transition from discrete avalanching (also named slumping) to continuous flow is currently poorly understood.
This work focuses on the effects of liquid viscosity and particle size on the dynamics and flow behaviour of wet granular avalanches. We used a discrete element method (DEM) model to simulate avalanches. The simulations provided extra information and allowed a wider range of research conditions than we could achieve experimentally. We used a limited set of experiments to validate our DEM model, which yielded comparable results in terms of average slope angle and avalanche amplitude, without any calibration of the simulation parameters. We examined the evolution of the slope angle, particle velocity, volume fraction and granular temperature, as functions of liquid viscosity. The transition from the avalanching to the continuous flow regime was then characterized using dimensionless quantities that compare particle and liquid properties. The obtained results have important implications for geo-mechanics and landslide control.

Drum setup
A schematic representation of the experimental setup is shown in Fig. 1. The setup consists of a drum, that was partially filled with granular materials and rotated at a speed of = 0.5 rpm.  This rotational speed proved low enough to observe intermittent avalanching flow in the dry and low viscosity cases. The drum is an updated version of the one used previously by Windows-Yule et al. (2016) and Jarray, Shi, et al. (2019). It is composed of a cylinder with a radius of R = 60.5 mm and a width of w = 22 mm. The bases of the drum consist of two clear circular plexiglass (PMMA) walls of 5 mm thickness to allow optical access. These were coated with fluorinated ethylene propylene (FEP) to prevent wet glass particles from sticking to the walls. The side of the drum is made of poplar wood.
In each experiment, the drum was filled to circa 35% of its volume with monodisperse borosilicate glass particles, with a density of = 2500 kg/m 3 and a radius of r = 1.25 mm. We studied a dry case, as well as several cases where the particles were mixed with 4 cm 3 of water-glycerol mixtures, leading to a pendular state (Gabrieli, Lambert, Cola, & Calvetti, 2012;Iveson, Beathe, & Page, 2002;Jarray, Shi, et al., 2019;Urso, Lawrence, & Adams, 1999). Here, a small amount of liquid is shared by neighbouring particles, but between them is mostly air. The surface tension and contact angle of the mixtures stay near-constant as a function of glycerol concentration, while the liquid viscosity strongly varies (Glycerine Producers' Association, 1963;Vicente, André, & Ferreira, 2012). Therefore, the induced viscous forces strongly increase with glycerol concentration, while the capillary forces remain relatively constant. Table 1 displays the properties of the mixtures. The case n = 0 indicates dry conditions. A subset of cases were reproduced with both experiments and simulations (see Section "Numerical model"), allowing validation. For liquids with > 370 mPa s, particles start to stick to the front wall of the rotating drum, obstructing accurate image post-processing, thus only simulations were performed. Table 1 also shows the values of the capillary number Ca, defined as with liquid viscosity , surface tension , contact angle ˛ and characteristic velocity v c . The capillary number is the ratio of the viscous force to the capillary force. Since the velocity of the particles in the flowing layer during an avalanche is mainly induced by gravity, the characteristic velocity was taken equal to the free-falling speed of a particle, after falling a distance of 2r, i.e. v c = 4gr, which yields an approximate high-end value of Ca. In Table 1, we used r = 1.25 mm for the calculation of Ca.
Image post-processing Images of the rotating drum were recorded using a Canon Legria HFG40 camera, operating at 50 FPS. The images captured by the camera were post-processed using Fiji ImageJ (Rueden et al., 2017;Schindelin et al., 2012) to obtain the slope angle Â and the avalanche amplitude Â, as depicted in Fig. 1. In the first processing step, the irrelevant parts of the image sequence (i.e. the background and the drum walls) were removed and the light and contrast were adjusted. Subsequently, the Trackmate (Tinevez et al., 2017) package was used to detect the particles based on the difference of Gaussian distributions (Lowe, 2004). The detected particles were stored in tabular form and imported into the data visualization program Paraview (Ayachit, 2015). Finally, an in-house Python script was employed to determine the slope angle and the avalanche amplitude.

Numerical model
We employed the software package LIGGGHTS (Kloss, Goniva, Hager, Amberger, & Pirker, 2012) to simulate wet particles in a rotary drum, using the discrete element method. Particle trajectories are calculated by solving Newton's equations of motion. The movement of an arbitrary particle i, that is in contact with N a adjacent particles, is described as with mass m i , position vector x i , contact force vector F c ij , cohesive force vector F l ij , gravity vector g, mass moment of inertia I i , angular position vector Â i , distance-to-contact point vector h ij and additional torque vector q ij (to account for rolling resistance). Note that variables printed in bold represent vectors.

Contact model
A spring-dashpot model (Di Renzo & Di Maio, 2005;Kloss et al., 2012;Silbert et al., 2001) was used to compute the force F c ij between two colliding particles i and j. F c ij is nonzero only if there is a positive normal overlap ı n between the particles, which is defined as where r is the particle radius and n ij the unit normal vector. When ı n > 0, F c ij is computed as

G Model
with tangential overlap ı , unit tangential vector ij , elastic spring constant k and visco-elastic damping constant . The subscripts n and denote normal and tangential components respectively. The Hertz-Mindlin model (Kloss et al., 2012;Mindlin & Deresiewicz, 1953), based on Hertz theory in normal direction and the Mindlin no-slip improved model in tangential direction, was employed to calculate the elastic and damping constants. k n is given by with modulus of elasticity E, Poisson's ratio and effective radius r * = r i r j / r i + r j . Furthermore, n is computed as with coefficient of restitution (CoR) e and particle density . As for the tangential components, k is defined as where G is the shear modulus and is defined as A constant directional torque model (Ai, Chen, Rotter, & Ooi, 2011;Kloss et al., 2012) was implemented to account for rolling resistance, and contributes with an additional torque vector q ij to Eq. (3), which is given by with coefficient of rolling friction (CoRF) RF . The vectorÂ rel is the relative rotational velocity between particles i and j. In the case of particle-wall contact, Eqs. (6)-(10) are taken in the limit of radius r j going to infinity.

Liquid cohesion model
A composition of models, as suggested by Easo and Wassgren (2013), was used to define the cohesive force between particles i and j. It is assumed that a surface liquid film exists on each particle, allowing for a liquid bridge to form upon contact with another particle. The bridge ruptures when x i − x j exceeds the rupture distance d 0 , which is given by with contact angle ˛ and liquid bond volume V bond (Lian, Thornton, & Adams, 1993). It is assumed that V bond equals 5% of the combined liquid film volume, which distributes evenly over the two particles after bridge rupture. The cohesive force vector F l ij between the particles is defined as with capillary force vector F cap ij and viscous force vector F vis ij .
The capillary force acts in normal direction and is defined by Soulie, Cherblanc, El Youssoufi, and Saix (2006) as with liquid surface tension and r max = max r i , r j . Parameters a, b and c are given by The viscous force consists of a tangential and normal component and is described by Nase, Vargas, Abatan, and McCarthy (2001) as with liquid viscosity . The value for the normal overlap is bounded by to prevent limitless development of the cohesive forces. In the case of particle-wall contact, F cap ij and F vis ij are set equal to zero. To summarize, the total force vector F ij can be decomposed in normal and tangential components as with total normal force vector F n ij and total tangential force vector F ij . Furthermore, the tangential component of the contact force with coefficient of friction (CoF) F , to ensure the satisfaction of Coulomb's law.

DEM simulation parameters
The model parameters in the DEM simulations were set to mimic the physical experiments. The material properties of the particles and drum were set to resemble glass and poplar wood respectively (Gray, 1972;Green, Winandy, & Kretschmann, 1999;Oberg, Jones, Horton, & Ryffel, 2000;Serway & Jewett, 2004;Sigmund Lindner, 2018), and are summarized in Tables 2 and 3 .
In all simulations, the drum was filled to 35% of its volume with monodisperse particles and rotated at a constant rotational speed of = 0.5 rpm, proven low enough to observe intermittent avalanches in the dry and low viscosity cases. The particles had a radius of  (-) 0.01 r = 1.25, 2.0 or 3.0 mm. We performed simulations under dry conditions, and simulations with seven water-glycerol mixtures, as presented in Table 1. We used a reduced CoF for wet particles to mimic lubrication effects in the granular system. In each simulation, we tracked the particle positions and velocities for 50 seconds, to examine the evolution of the slope angle and velocity profiles. The software package ParaView (Ayachit, 2015) and in-house python codes were used to evaluate the slope angle Â, the avalanche amplitude Â and the packing fraction , defined as the volume of all particles divided by the volume of the packing. Furthermore, we evaluated the granular temperature T g (Campbell, 2006;Zhou & Sun, 2013), defined as where N p is the total number of particles. The fluctuation velocity of particle i is given by where u i is the velocity of particle i, andū i the mean velocity of the particles surrounding i, within a spherical region of radius 7r i /3. Thus, U i indicates a particle's velocity with respect to its surrounding particles. The velocity profiles of the particles in the bed were obtained using a Gaussian Kernel point volume interpolation scheme (Jarray, Magnanimo, & Luding, 2019;Price, 2012;Schroeder, Lorensen, & Martin, 2004), with a kernel radius of 8 mm for the cases with r = 1.25 mm and r = 2 mm, and a kernel radius of 11 mm for the cases with r = 3 mm.

Model validation
To validate our numerical model, we first compare the simulated and measured time-averaged surface angle Â and time-averaged avalanche amplitude Â . Fig. 2(a) shows Â as a function of the interstitial liquid viscosity . Both experiments and simulations indicate that has no significant effect on Â in the investigated viscosity range. However, wet particles yield a higher Â in comparison to dry particles, due the formation of capillary bridges. In all wet cases, simulations yield a slightly lower averaged slope angle than the experiments. In the DEM model, the liquid bridge volume was set to always distribute evenly over both particles after bridge rupture. Thus, we propose to attribute the observed difference between simulations and experiments to the occurrence of in-homogeneous liquid migration between particles in the experiments. Fig. 2(b) shows the averaged avalanche amplitude as a function of . We observe a decreasing trend of Â with liquid viscosity in both the numerical and experimental cases.
Comparable results are obtained from the experiments and numerical simulations, which confirm the validity of the adopted model. Since experimental research did not allow investigation of the avalanche behaviour for high viscosity liquids ( > 370 mPa s), as discussed in Section "Experimental method", and to get more insights into the micro-mechanics of granular avalanches, we will henceforth focus on the results obtained from the DEM simulations.

Simulation results
Time evolution of the granular flow Fig. 3(a) shows a typical evolution of the surface angle during several successive dry avalanches. Maximum stability and repose angles Â m and Â r respectively correspond to the local maxima and minima in the "sawtooth"-like curve. The avalanches in this case are very pronounced; rigid body movement of the particulate bed is alternated with particles flowing down in the surface layer. This is further illustrated by the evolution of the granular temperature T g , shown in Fig. 3(b). We observe peaks in T g during the avalanches, which indicate large relative movements between particles, that coincide with drops in Â. In between avalanches, T g drops to nearly zero, while the particulate bed is lifted again as a rigid body. Qualitatively similar behaviour is observed for wet granular avalanches with an interstitial liquid of low viscosity ( = 1.01 mPa s, n = 1), shown in Figs. 1(c)-(d) by the continuous blue line. However, for high viscosity ( = 1410 mPa s, n = 7), the evolution of the surface angle does not show a steady periodic trend. Smaller fluctuations occur more frequently, while distinct avalanches disappear, implying a more continuous flow. We observe a similar trend in the development of T g . In comparison to the dry case n = 0 and low viscosity case n = 1, peaks for the case n = 7 are significantly smaller. Furthermore, T g does not drop as significantly with the increase of Â, indicating that there is always relative movement between particles, and the particulate bed never behaves as a rigid body.
An interesting feature of the avalanches shown in Fig. 3(a) and (c) is the rather large variation in the amplitudes Â of successive avalanches. A slower rotational speed of the drum is likely to reduce these variations and provide an even more periodic trend (Balmforth & McElwaine, 2018;Fischer et al., 2008).
Appendix "Additional data for other particle sizes" depicts the time evolution of Â and T g for simulations with r = 1.25 mm and r = 3 mm. In these cases we observe qualitatively similar behaviour to those with r = 2 mm shown in Fig. 3. Fig. 4(a)-(d) shows typical velocity profiles for avalanches with low to high liquid viscosity. In all cases, a high velocity flow layer near the surface is observed, while below this layer particles are at rest with respect to the drum, i.e. only move due to its rotational speed . The dry case n = 0 yields a larger flow layer velocity v than the wet cases, for which v decreases with liquid viscosity. This result has also been observed by Brewster, Grest, and Levine (2009), who found that cohesive forces limit the development of the velocity in the flow layer. Fig. 4(e) shows the particle velocity along the pile depth d, which is along the inward-pointing normal to the free flowing surface, starting from its middle. In all cases, the velocity of the flowing layer decreases slowly with d initially, near the surface; then more rapidly, in the middle of the flowing layer; but then declines, until the velocity of the particles is purely due to the rotational speed of the drum. The slope of the velocity curve decreases with viscosity, and is steepest for the dry case. However, the thickness of the flow layer appears to be independent of liquid viscosity. The velocity profiles for particle radii of r = 1.25 mm and r = 3 mm, presented in Appendix "Additional data for other particle sizes", show qualitatively similar behaviour to those with r = 2 mm G Model PARTIC-1401;No. of Pages 12 J.H. Kasper et al. Particuology xxx (xxxx) xxx-xxx  shown in Fig. 4. The flow layer velocity was observed to increase in magnitude with particle size, yet the overall trends remain similar. Fig. 5(a) shows the time-averaged avalanche amplitude Â as a function of liquid viscosity, for several particle sizes. Note that as increases, Â becomes less a measure of avalanche amplitude, and more a measure of surface angle fluctuation in time, as the flow transitions from the avalanching regime to the continuous regime. The avalanche amplitude is shown to increase with particle radius r. This trend is consistent with experimental research on dry (Balmforth & McElwaine, 2018) and wet (Courrech du Pont et al., 2003) granular avalanches, in which larger avalanches are observed as the ratio of particle over drum radius r/R increases. For all particle sizes, Â is observed to decrease with viscosity and saturate at high values of , in agreement with statements by Courrech du Pont et al. (2003).

Avalanche amplitude
We can exclude the effect of the particle size on the wet avalanche amplitude by considering its dimensionless form Â wet / Â dry . As a result, the normalized avalanche amplitude data collapse on a single viscosity-dependant master-curve, plotted in a dashed line in Fig. 5(b). The fitting curve is described by where a and b are fitting constants equal to 0.4 and 0.47, respectively. Furthermore, c is a characteristic viscosity equal to 295 mPa s. We infer that for < c the decay of Â wet / Â dry is large, while for > c , the dimensionless amplitude saturates.

Slope curvature
The curvature of the surface can be used to differentiate between avalanching and continuous flow behaviour. For granular flows in the avalanching regime, the surface angle is generally observed to be constant along the surface (Fischer et al., 2008;Jaeger et al., 1989). For continuous flows however, the surface is observed to be "S-shaped" if the cohesion is low and > 5 rpm (Jarray, Magnanimo, et al., 2019), or even slightly convex if cohesive forces are large (Brewster et al., 2009). Significant fluctuations in the surface angle (i.e. the deviation between the local surface angle and the average angle of the whole surface) are indicative of continuous flow behaviour.
We can quantify these fluctuations using the following method. First, we approximate which particles are part of the surface layer by dividing the drum into N bins and considering the top particle in each bin. Polyfitting a linear function through the coordinates of all surface particles yields a flat averaged surface (see inset of Fig. 5(c)). Fluctuations in surface angle are evaluated by measur- ing the average deviation h, of the distance between the surface particles and the polyfitted surface. We thus define h as in which (x i , y i ) denote the coordinates of particle i, and (x fit i , y fit i ) denote points on the polyfitted surface. Fig. 5(c) shows that indeed the time-averaged normalized surface deviation h wet / h dry , slightly increases with liquid viscosity, and consistently decreases with increasing particle size.

Granular temperature
The time-averaged granular temperature T g is plotted in Fig. 6(a), and shows a similar trend as the time-averaged amplitude Â . T g decreases with viscosity, but increases with particle size. As r/R increases, T g is evaluated over a larger local volume, thus larger relative movements are indeed expected. The granular temperature of the dry cases is of an order of magnitude greater than for the wet cases. The liquid-induced cohesive bonds between particles limit relative particle motion. Particles clump together, and the number of inter-particles collisions reduces. This effect gets magnified as increases and the cohesive bonds become stronger. Interestingly, the curves for r = 1.25 mm and r = 2 mm are quite smooth, while for r = 3 mm some outliers appear. This is likely due to interlocking of particles, caused by the high particle over drum size ratio (r/R ≈ 0.05). Liu, Specht, and Mellmann (2005) pointed out that the avalanching regime is marked by intermittent rigid body movement of the particulate bed, during which T g drops to near-zero (see Fig. 3). As the liquid viscosity increases, the time intervals of rigid body movement become shorter, and the flow transitions from the avalanching regime to the continuous regime. In complete absence of these intervals, the flow can be characterized as fully continuous. In order to investigate the transition from avalanching flow to continuous flow, we evaluate avalanching time ratio t a /t s . Here t a is the total avalanche duration (i.e. the duration for which T g ≥2 ·10 −9 m 2 /s 2 , corresponding to non-rigid body movement of the bed), and t s is the sampling duration. Fig. 6(b) depicts t a /t s as a function of viscosity, where t s was set to the final 30 s of the simulation. We find that t a /t s increases logarithmically with liquid viscosity. For particle sizes r = 1.25 mm and r = 2 mm, continuous flow develops beyond = 520 mPa s and = 940 mPa s, respectively. The corresponding average avalanche amplitudes are below 2.5 • . However, for r = 3 mm, continuous flow does not seem to develop within the range of viscosities examined in our research.

Dimensional analysis and flow transition region
In order to estimate the relative importance of gravitational forces with respect to the viscous forces acting on wet particles, we introduce the granular Galileo number: The Galileo number is an especially relevant parameter, as the displacement of particles during avalanches is dominated by gravitational forces. In Fig. 7(a), we plot t a /t s as a function of the square-root of the inverse of Ga g . From Fig. 6(b), we found that continuous flow develops only for cases with r = 1.25 mm and r = 2 mm, beyond = 520 mPa s and = 940 mPa s, respectively, where t a /t s ≈ 1. We are able to define a region in which the transition from avalanching flow to continuous flow takes place, marked by the blue zone in Fig. 7(a). This region roughly corresponds to a granular Galileo number Ga g ≈ 1, i.e. the viscous forces start to surpass the gravitational forces at this point. We also observe that the data set collapses onto a single hill-type curve, that tends to t a /t s = 1 at high values of 1/Ga g . Having established the avalanching-continuous flow transition zone in terms of the Galileo number, we subsequently explore the effect of this number on the avalanche amplitude in Fig. 7(b). For 1/Ga g < 10 −2 , the flow behaviour is almost unaffected by the viscous forces, but rather governed by capillary and gravitational forces, resulting in a relatively high avalanche amplitude. As 1/Ga g increases, the flow becomes more dilated and the avalanche amplitude decreases, particularly for 1/Ga g > 0.5. Average surface angle and volume fraction Fig. 8(a) shows the time-averaged surface angle Â , with standard deviations, as a function of , for several particle radii. For the dry case, Â increases with particle size, in agreement with previous observations (Balmforth & McElwaine, 2018). For the wet cases, Â shows a slightly decreasing trend with r, which is in agreement with observations made by Nowak, Samadani, and Kudrolli (2005). Again, some outlying points appear for simulations with r = 3 mm, possibly the result of interlocking. We observe no clear dependency of the average angle on the liquid viscosity .
Finally, the time-averaged volume fraction of the particulate bed is plotted as a function of the interstitial liquid viscosity in Fig. 8(b). Here, is defined as the summed volume of the particles over the volume occupied by the packing. We observe that wet particles yield a higher volume fraction than dry particles, indicative of a more compressed system (Shi et al., 2020). Overall, it seems that viscosity only has a weak effect on . However, decreases with increasing particle size. This is likely due to the increasing particle over drum radius ratio r/R, which negatively impacts packing efficiency.

Conclusions
Using discrete element method simulations, we have investigated how the dynamics of wet granular avalanches in a rotary drum are affected by the viscosity of the interstitial liquid and particle size. Simulations were validated against a limited set of experimental measurements and good agreement was obtained.
Discrete avalanches were observed for the dry case and the wet cases with relatively low viscosity. As viscosity increases, the flow behaviour transitions from the avalanching to the continuous regime. This transition is indicated by a decrease in avalanche amplitude, flow layer velocity and granular temperature and an increase in curvyness of the flow surface. The avalanching regime is marked by intermittent rigid body movement of the particulate bed, and near-zero drops in the granular temperature. For continuous flow, rigid body movement of the bed does not occur and the granular temperature assumes finite values well above zero. Based on the evolution of the granular temperature and avalanche amplitude, we were able to identify an avalanching-continuous flow transition region, as a function of the dimensionless granular Galileo number Ga g .
As the particle size increases, the avalanche amplitude, flow layer velocity and granular temperature increase as well, indicating an increase in relative importance of inertial effects with respect to frictional forces and viscous forces.
Our findings suggest that the effects of liquid viscosity and particle size should always be considered when studying granular avalanches, particularly for 1/Ga g > 0.5, whether large scale investigation, pilot scale experiments or numerical simulations are of interest.
The analysis presented in this paper could be extended to correlate the flow behaviour (i.e. the angle of repose or the avalanche Figs. 9 and 10 show the time evolution of the surface angle Â and the granular temperature T g for cases with r = 1.25 mm and r = 3 mm, respectively. Similar to cases with r = 2 mm, shown in Fig. 3, discrete avalanches are observed for dry cases and low viscosity cases ( = 1.01 mPa s). Highly viscous liquids ( = 1410 mPa s) yield significant reductions in surface angle fluctuations and fluctuations in T g , and avalanches are less pronounced. This effect appears less dominant for larger particles however, as the r = 3 mm case with high viscosity still shows distinct avalanches.
Figs. 11 and 12 show the velocity profiles for avalanches with low to high liquid viscosity for cases with r = 1.25 mm and r = 3 mm, respectively. The velocity of the flowing layer increases with particle size, while its thickness remains relatively constant.