{'fraction_non_alphanumeric': 0.039183291770573565, 'fraction_numerical': 0.013185785536159601, 'mean_word_length': 4.524539349061477, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Finite-Size Analysis of the Collapse of Dry Granular Columns

In this letter, we focus on the size effect of granular column collapses, which are potentially connected to the dynamics of complex geophysical flows, even if the link between microscopic structures of granular assemblies and their macroscopic behaviors is still not well understood. Using the sphero-polyhedral discrete element method (DEM), we show that the column radius/grain size ratio has a strong influence on the collapse behavior. A finite-size analysis, which is inspired by a phase transition around an inflection point, is performed to obtain a general scaling equation with critical exponents for run-out distances. We further link the size effect with the strong force network and formalize a correlation length scale that exponentially scales with the effective aspect ratio. Such a scaling solution shows similarities with the percolation problem of two-dimensional random networks and can be extended to other similar natural and engineering systems.


Introduction
Granular materials are omnipresent in natural and engineering systems, and the physics and mechanics of them are crucial for understanding geophysical flows, natural hazards, food processing, chemical engineering, and pharmaceutical engineering (Lube et al., 2019;Guyon et al., 2020). Granular materials can behave like a solid, a liquid, or a gas in different circumstances (MiDi, 2004), which increases the difficulty in capturing their macroscopic behavior. Besides, collective structures may form inside a granular system, and the existence and the size of such collective structures [such as bridging (Mehta, 2007), granular agglomerates , and contact networks (Zhang et al., 2014)] will introduce a strong size effect, which further increases the complexity of the problem. In recent decades, breakthroughs have been made to understand the basic governing principles, especially the constitutive relationships, of granular materials (MiDi, 2004;Trulsson et al., 2012;Man & Hill, 2021), where the behavior of granular materials or granular-fluid mixtures are considered to be described by dimensionless numbers expressed as the ratio between dominating stresses. Even though geophysical flows are rarely dry, investigating dry granular materials is still considered as an important start for more complex systems.
Due to the similarity and potential links between the collapse of granular columns and gravity-driven geophysical flows, previous research investigated the collapse of granular columns to analyze the post-failure behavior of granular systems (Thompson & Huppert, 2007;Lacaze & Kerswell, 2009). Lube et al. (2004) and Lajeunesse et al. (2005) independently determined relationships for both the normalized run-out distance R = (R ∞ −R i )/R i (where R ∞ is the final radius of the granular pile, and R i the initial radius of the granu-lar column), and the halt time of a collapsed granular column, which both scale with the initial aspect ratio, α = H i /R i (where H i is the initial height of the column), a parameter drawn out of dimensional analysis. Zenit (2005); Hinch (2005, 2007); Lacaze and Kerswell (2009) studied either 2D or axisymmetric granular column collapses using the discrete element method (DEM), and explored the influence of inter-particle properties, energy dissipation, and inertial effects on the run-out and deposition of granular columns. Farin et al. (2019) even linked the characteristics of granular column collapses to seismic signals. Warnett et al. (2014) and Cabrera and Estrada (2019) studied the collapse of granular columns with experiments and simulations, and argued that the relative size of a granular column, R i /d, where d is the average particle diameter, has strong influence on the run-out distance of a collapsed granular column, and to avoid significant size effects, R i /d must be larger than 75 for short columns and larger than 50 for tall columns.
Based on dimensional analysis and simulations, we investigated the collapse of axisymmetric granular columns and their resulting deposition with a wide range of inter-granular and particle/boundary frictional coefficients , and concluded that the normalized run-out distance, R, scales with an effective aspect ratio, where µ w is the frictional coefficient between particles and the boundary, µ p is the interparticle frictional coefficient, and β can be seen as the ratio between contributions from inter-particle frictions and those from particle/boundary frictions. They found that the collapse of granular columns can be classified into three different regimes: quasi-static, inertial, and liquid-like. The effective aspect ratio, which was used to describe the run-out behavior, was obtained from dimensional analysis and was shown to represent the ratio between the inertial stress and the frictional stress or, equivalently, the ratio between the kinetic energy and dissipated energy during the collapse of granular columns. In our study, as we change the frictional coefficient, the transition point, where the slope of R(α) curve changes on a log-log coordinate system, also changes accordingly. Moreover, after using α eff to re-scale the data, this transition point in R − α eff space becomes universal.
In this paper, after introducing the simulation set-up and the contact law in the DEM, we further conduct a systematic study of granular column collapses with the spheropolyhedral DEM to investigate the influence of relative column size on the run-out distance to obtain a scaling solution to describe granular column collapses with different sizes. This can not only further explore the physics of granular column collapses but also open a window for linking the behavior of granular columns to the mobility of geophysical flows. Inspired by the work of Roche et al. (2002Roche et al. ( , 2008; Bougouin et al. (2021), where they presented physical connections between granular dam breaks and pyroclastic flows, to strengthen the tie between granular column collapses and geophysical flows, we present an application of our results to various rock falls around the Soufriere Hills, Montserrat as a result of volcanic eruptions and the resulting pyroclastic flows and granular surges in Section 4.2 and Figure  3. Finally, we present our concluding remarks.

Discrete element method and simulation setup
We perform DEM simulations of the granular column collapses with Voronoi-based sphero-polyhedra particles (Pournin & Liebling, 2005;S. A. Galindo-Torres & Pedroso, 2010). The sphero-polyhedra technique allows for an easy and efficient definition of contact laws between particles. This is due to the smoothing of the edges of all geometric features by circles (in 2D) or spheres (in 3D). We note that, in our simulation, three types of contacts (vertex-vertex contact, edge-edge contact, and vertex-face contact) are considered. For these types of contacts, we implement a Hookean contact model with energy dissipation to calculate the interactions between particles. The DEM model used in this work has been thoroughly validated with experiments as shown in the supplemental material, where the shape of particles was accurately captured and the frictional coefficients of both inter-particle contact and particle/boundary contact were independently measured. (c) shows the method for measuring the run-out distance of a collapsed granular column. The x−axis is the radial position r, and the y−axis is the percentage of number of particles, P dm , located within (r − ∆r/2, r + ∆r/2) divided by the radial position.
We determine the final radius of a collapsed granular column when P dm (r) ≤ 5%P max dm ; (d) Initial state and (e) final state of granular column with Ri = 5.7 cm and Hi = 13 cm; (f) Initial state and (g) final state of granular column with Ri = 2.3 cm and Hi = 10 cm.
In the simulation, we create particles in a certain cylindrical domain with initial height, H i , and initial radius, R i [ Fig. 1(a)]. 20% of the particles are randomly chosen and removed from the simulation to create a granular packing with an initial solid fraction of φ s =0.8. The average particle size d is 0.2 cm. The coefficient of restitution of particle collisions is 0.1. The material properties are set to be the same as that of quartz sand (particle density 2.65 g/cm 3 ). Using 3D Voronoi algorithm to create sphero-polyhedral particles, we can make sure that, for each particle, the length scales in different directions are kept approximately the same ( as shown in the supplemental material). This paper focuses entirely on simulations with granular columns of circular cross-sections. More general cross-sections are currently under investigation and will be presented in a subsequent publication. The particle/boundary frictional coefficient µ w is 0.4, while the inter-particle frictional coefficients µ p are 0.2, 0.4, and 0.6.
We treat the relative column size (R i /d = 2, 2.5, 3.75, 5, 7.5, 10, 12.5, 15, 17.5, 20, 30) as a key parameter. Within one set of simulations, with the same R i /d, we varied the initial height H i and the inter-particle frictional coefficient µ p to obtain the general collapse behavior with a wide range of initial conditions. After the initial Voronoi-based particle packing is generated, we delete the confining tube in the simulation to release all the particles under gravity. After releasing particles to the horizontal plane, the granular material will form a pile of loosely packed grains [ Fig. 1

Determination of the run-out distance
In simulations, the measurement of the final packing radius is complicated due to the lack of a clear mathematical definition of the run-out distance. In cases with small particle/boundary and inter-particle frictional coefficients but large initial aspect ratios, the spread of particles is far-reaching and leads to sparse (single layer) coverage of the area, especially at the front edge. In these cases, it is difficult to determine the edge/boundary and hence the final run-out distance. Thus, we measured the final radius with a histogram of particle distribution for each simulation. Fig. 1(c) gives an example of how we measured the run-out distance.
In Fig. 1(c), the x-axis is the radial position r, and the y-axis is the percentage of number of particles located within (r − ∆r/2, r + ∆r/2) divided by the radial position, P dm (r) = (1/r)[N (r − ∆r/2, r + ∆r/2)/(Σ r N )], where ∆r is the bin width of the histogram, N (r − ∆r/2, r + ∆r/2) is the number of particles located between r − ∆r/2 and r + ∆r/2, and Σ r N is the total number of particles in one simulation. FIG. 1(c) shows the normalized particle number distribution, in this case, deposition morphology, of a simulation with µ w = 0.4, µ p = 0.4, R i = 4 cm, and H i = 32 cm. It shows that most particles locate within r ≤ 27 cm. Thus, we take the final deposition radius as R ∞ = 27 cm. After determining the final run-out distance of the collapsed granular column, we obtain the normalized run-out distance, R, the initial aspect ratio, α, and the effective aspect ratio , 3 Experiments Warnett et al. (2014) and Cabrera and Estrada (2019) examined the size effect in the granular column collapses with experiments. In this work, we present two sets of experiments to show that changing sizes could influence both the R − α relationship and the transition point. We perform column collapses of sand particles with grain sizes ranging from 1 mm to 2 mm (the average particle diameterd ≈ 1.5 mm). The experimental setup is the same as that in the simulations. The column diameters are 4.6 cm and 11.4 cm, respectively, resulting in relative column sizes R i /d ≈ 15.3 and 38 [ Fig. 1 We vary the initial height of the column so that the initial aspect ratio varies from ≈ 0.3 to ≈ 20. After randomly dropping sand grains into the cylindrical tube, we measure the initial height, H i , of the column, and then quickly lift the tube upward to let the granular column collapse. We measure the final radius of the sand pile in four different directions and take their averages as the final run-out distance, R ∞ . Then, the relationship between the initial aspect ratio, α, and the normalized run-out distance, R, can be obtained accordingly.
We plot the experimental relationship between the initial aspect ratio and normalized run-out distance in Fig. 2(a). For granular column collapses with both column sizes, as we increase the initial aspect ratio, the normalized run-out distance, R increases accordingly. We can then observe the clear transition point at α ≈ 1.2 for cases with R i = 5.7 cm and at α ≈ 2.0 for cases with R i = 2.3. The experiments show that the transition point becomes smaller as we increase the relative column size, and generally, the normalized runout distance of large columns is larger than that of small columns. The experiments show clear size effect for granular column collapses. It is an experimental proof that changing relative column sizes can not only influence the run-out distance but also affect the transition point on the R − α relationship. We note that the friction between sand particles and the confining tube while pulling it up might also influence the experimental results, but the influence is not significant, and it does not affect the transition we observe in the R − α relationship. In the supplemental material, we also present a set of experiments and corresponding simulations of the collapse of cubic particle assemblies. The agreement between experimental results and simulations shows that our DEM model could well capture the behavior of gravity-driven transient granular flows. the relationship between R and α eff with relative column size Ri/d = 2.5 and 17.5, respectively, to show that the critical inflection point (αc, Rc) changes as we change the relative size of the column; (d) shows the relationship between transitional aspect ratio, αc, and the relative column size, Ri/d. (e) plots the transitional normalized run-out distance, Rc against Ri/d; (f) shows the relationship between R(Ri/d) β 1 /ν and A = (α eff − αc∞)(Ri/d) 1/ν , with ν = 1.39 and β1 = 0.28.
The inset of the figure plots the same relationship in a log-log coordinate system.

Size effects of granular column collapses
After simulating granular column collapses with various column sizes and frictional coefficients, we obtain the relationship between R and α. Similar to what we have seen previously , for columns with the same relative size R i /d, less friction leads to larger run-out distances. As we change the x−axis to the effective aspect ratio, α eff , normalized run-out distance with the same R i /d collapse onto one curve in R − α eff space [ Fig. 2(b) and (c)]. The previous work indicates that the effective aspect ratio, which includes the influence of both initial aspect ratio and particle frictional properties, could describe the run-out behavior of the granular column collapse, and most importantly it shows that, for columns with the same relative column size, there is only one transition point for certain column size in the R − α eff relationship. The collapsing of the transition point of the R − α eff relationship indicates a possible phase transition related to the collapsing dynamics of granular columns. However, with various relative column sizes, we see that, in both Fig. 2(b) and (c), simulation results, especially the transition point in each R − α eff relationship, varies as we change R i /d. For granular columns with the same initial aspect ratio, granular columns with larger relative column size R i /d have longer run-out distances. In Fig. 2(b), we plot simulation results for all the system sizes, and simulation results with R i /d = 2.5 and R i /d = 17.5 in Fig. 2(c). The latter gives us a clearer picture of how changing relative system size could result in variations in the relationship between R and α eff and variations in the transition point (α c , R c ). This shows that the collapse of the granular column has a significant size effect, which was also observed by Warnett et al. (2014) and Cabrera and Estrada (2019). However, no quantitative studies has been made to universally include different frictional coefficients and boundary conditions to describe the run-out behavior and few researchers pointed out the physical nature of the initial aspect ratio α other than obtaining it from dimensional analysis .
Further, changing the relative column size also fundamentally influences the shape of R − α eff . For simulations with the same R i /d, a threshold α c of α eff exists to divide the R − α eff relationship into two groups [ Fig. 2(b)]. When α eff < α c , R approximately scales with α eff proportionally. When α eff > α c , R approximately scales with (α eff ) 0.5 with a rather sharp division between the two. Here, the transitional aspect ratio α c and the corresponding transitional normalized run-out distance R c can be seen as the critical aspect ratio and the critical run-out distance, respectively. Both α c and R c vary as we change the size of the granular column. For instance, in Fig. 2(c), the transition happens at α c ≈ 8, R c ≈ 5 when R i /d = 2.5, and happens at α c ≈ 2, R c ≈ 3.5 when R i /d = 17.5.
We show, in Fig. 2(b) and (c), that a critical transition point (α c and the corresponding R c ) exists and also varies with different relative size R i /d. Thus, the R−α eff curve is dictated by the position of α c and R c that Both α c and R c decrease with an increasing R i /d. In Fig. 2(d) and (f), based on the relationship between α eff and R in Fig. 2(b) and (c), we plot the relationship between α c and R i /d and the relationship between R c and R i /d. Thus, we could write α c and R c as functions of the relative column size, where, when R i /d approaches infinity, both α c and R c converge to α c∞ and R c∞ , as shown by where α c∞ = 0.2 and R c∞ = 0.732 are the fitted critical aspect ratio and the corresponding critical run-out distance, respectively, when the relative column size R i /d goes to infinity, and both parameters are fitted values. Also, a 1 = 16, b 1 = −0.75, a 2 = 5.4, and b 2 = −0.3 are fitting parameters. The fitted curves of α c and R c are the dashed curve in Fig. 2(d) and the dash-dot curve in Fig. 2(e).
The power-law decay of both α c and R c with respect to R i /d inspires us to perform a finite-size analysis of the run-out distance of the collapse of granular columns based on our working hypothesis of a potential phase transition. Since a general scaling analysis often works under thermodynamics limits, where the system size and its corresponding particle number is assumed to be infinity. Finite-size analysis then becomes a powerful tool for investigating real systems with finite sizes. As shown in Fig. 2(f), the normalized runout distance, R indeed exhibits excellent finite-size scaling, suggesting that the transitional aspect ratio, α c , is critical. All the normalized run-out distance data collapse nicely onto a master curve in the form (S. A. Galindo-Torres et al., 2015) with the limiting scaling of the normalized run-out distance being R ∼ α eff − α c∞ , where ν = 1.39 ± 0.14 and β 1 = 0.28 ± 0.04 are obtained to best collapse all the data. In Fig. 2(f), to better represent the results, we took A = (α eff − α c∞ )(R i /d) 1/ν as the x−axis. Since the system is axisymmetric, the scaling solution for the collapse of granular materials shows similar phenomena in the scaling solutions for connectivity of two-dimensional continuous random networks (S. A. Galindo-Torres et al., 2015), where the scaling parameters ν = 1.33 and β = 0.1389. This indicates that the flowing behavior of granular column collapse has strong correlations with the connectivity of grain contact networks. This also coincides with the suggestion made by Mehta (2007) that the transport of grains takes place percolatively, especially in the context of avalanches.
We then investigate the strong force network of columns with different relative sizes at the very beginning of collapses [ Fig. 3(a)-(d)]. In Fig. 3(a)-(d), we only plot the contact forces f that are larger than the mean contact force f linking the centroids of contacting particle pairs. The blue dashed lines represent the height that the strong force network can reach in the vertical direction, while the red dashed rectangles denote the region which the strong force network could occupy the whole cross-section. As we increase the relative size of the column, the height of the strong force network remains almost unchanged. However, the horizontal occupation of the strong force network is strengthened with the increase of the relative column size. This indicates that larger columns could form better-connected strong force networks in the horizontal direction at the beginning of the column collapse, and therefore might have the capacity to drive the system to flow further than smaller columns.
Equation 5 and the strong force network imply that a length scale exists during the collapse of granular columns in a form that the relative size of the length scale with respect to the column height has a power-law relationship with α eff − α c∞ of the form where ζ is the critical exponent for the length scale ξ, and ξ/H i denotes the degree of occupation of sheared granular media along the height of the initial granular column. To further investigate the length scale ξ associated with the collapse of granular columns, we choose four sets of simulations with different relative system sizes. We consider grain pairs along the vertical direction at the very beginning of the collapse when time t = 0.05 s, and define the corresponding correlation G(z), which describes the correlation between particle pairs with vertical distance z, based on the equation proposed by Barker and Mehta (Barker & Mehta, 1992;Mehta, 2007), by where ∆z i and ∆z j are the vertical displacements of particle i and j after one time step (e.g. ∆z i = v i ∆t where v i is the velocity of particle i and ∆t is the time step), z ij is the vertical distance between two particles, z is the vertical displacement between particles, and d ij = (x i − x j ) 2 + (y i − y j ) 2 is the horizontal distance between two particles. Here δ( ) is a rectangular function, where the function is equal to one when z − 0.5d < |z ij | ≤ z + 0.5d and equal to zero elsewhere, and Θ( ) is the Heaviside step function. This definition ensures that the averages run over all displacements of sphere pairs in the vertical direction. The correlation function G(z) shows how particles influence each other in the z−direction.
The displacement correlation function can be fitted with an exponential equation G(z) = Aexp(−z/ξ), where ξ can be seen as a correlation length scale associated with the displacement correlation among particles during the collapse of granular columns, and A is a fitting parameter. The occupation of correlated grains across the height of a granular column, represented by ξ/H t , where H t here is the height of the column at the time when we measure the correlation length scale (since ξ is measured at the very beginning of the collapse, H t ≈ H i ), indicates the collective motion of particles during a collapse of granular columns. Thus, we plot the relationship between ξ/H t and α eff − α c∞ in Fig. 3(e), and determine Compares simulation data to field data presented in Calder et al. (1999). The red markers denote our simulation results, while other markers represent the data compiled from Calder et al. (1999).
that the ratio between correlation length scale and the system height shows power-law decay as we increase α eff − α c∞ , which is consistent with the finite-size scaling of both run-out behavior of granular column collapses, and ζ = 0.89 in Eq. 6 best fits the power-law decay. When α eff approaches α c∞ , ξ/H t starts to deviate from the power-law relationship. This is because when the initial aspect ratio is too small, only several layers of particles are present along the height of a granular column.
Additionally, we calculate the correlation length scales, ξ, of the four cases shown in Fig. 3(a)-(d), and plot them against the ratio between the height of the fully occupied strong force region h s and the column height H t [ Fig. 3(f)]. We find that, as we increase the relative column size, not only does the height of the fully occupied strong force region increase, but the correlation length scale also increases accordingly. This implies that a larger system could help establish a better correlated moving mass, and as we increase the initial aspect ratio, columns with larger relative size, where well-correlated moving mass is easier to establish, would experience the phase transition from quasi-static collapse to inertial collapse earlier than smaller systems. To further analyze the phase transition associated with granular column collapses, we would have to conduct rigorous investigations on the mean-field theory, which we plan to do in future studies.

Application to geophysical flows
In order to present an application of granular column collapse results to a real geological situation, we also compare our simulation data to the rock falls around Soufriere Hills Volcano, Montserrat as a result of volcanic eruptions and the resulting pyroclastic flows and granular surges [presented in Calder et al. (1999)]. Although granular column collapses do not behave exactly the same as landslides, rock avalanches or pyroclastic flows, their underlying physics is similar since the mobility of collapsed materials is mainly governed by the ability of transforming potential energy into kinetic energy, by surpassing the dissipative action of friction, a competition of effects that is summarized in the proposed α eff as explained before. Roche et al. (2002Roche et al. ( , 2008 also argued that the dam-break experiments could reveal similar behavior as inertial geophysical flows and the results supported the pyroclastic flow modeling approach by Levine and Kieffer (1991). In Calder et al. (1999), R ∞ /H i , which indicates the ability of relocating materials, is plotted against a parameter, Ψ = ρgV /H 2 i , where V is the volume of the material being transported. In Fig. 3(g), we plot our simulation results along with the results presented in Calder et al. (1999), where the markers indicate our simulation results, while other markers represent the field data.
We see that most of simulation data do not fall into the region of field observations. This is mainly due to the length scale of simulations being much smaller than a typical geophysical flow. However, as we increase the system size of simulations, we do obtain granular collapses in the same regime as some of the pyroclastic flows and granular surges. The finite-size analysis early in this section implies that increasing the system size could lead to an asymptotic behavior of granular flows, where the transition point α c → α c∞ . We can also see this trend in Fig. 3(g) where, as we increase the system size of simulations, our simulation data merge into the field results. More intriguingly, this also shows that almost all the data acquired from dome-collapse pyroclastic flows [ in Fig. 3(g)] fit in the trend of simulation results; this type of pyroclastic flows is considered to have a flow behavior more similar to granular flows. The size effect of granular column collapses might be potentially connected to predict the behavior of geophysical flows with different length scales. More importantly, the scaling of run-out distances with respect to R i /d might help explaining the long-range run-out in some pyroclastic flows (Lube et al., 2019), where the initial radius of the pyroclastic "column" is much larger than the volcanic particles. The analysis of force networks within granular systems, as mentioned in Section ??, may be helpful to explain some phenomena of geophysical flows. However, it is still difficult to directly obtain the contact information within real geophysical systems.

Conclusions
Previous research concluded that we could combine the influence of initial and boundary conditions and the initial column aspect ratio to determine the normalized run-out distance with an effective aspect ratio, α eff , which introduced a universal relationship to link the behavior of granular column collapses in three different collapsing regimes (quasi-static, inertial, and liquid-like) . In this paper, we further investigated the size effect associated with the collapse of granular columns. Our research is performed with DEM simulations of Voronoi-based spheropolyhedron particles. We found that the transition point in R − α eff space, which distinguishes the inertial collapse regime from the quasi-static collapse regime, varies as we change the relative system size R i /d. Both R c and α c experience power-law decay with respect to R i /d, which implies possible finite-size scaling for the normalized run-out distance R.
Similar to the finite-size analysis of the jamming transition of granular materials (Liu et al., 2014), where the stress scales with the solid fraction φ − φ c∞ , we took the previously determined effective aspect ratio α eff as the key parameter, and discovered that both the runout distance of granular column collapse and the energy consumption of it follows strong finite-size scaling. Interestingly, the scaling parameters we discovered for the finite-size scaling of the run-out distance of granular column collapses are similar to those presented in percolation problems of two-dimensional random networks (S. A. Galindo-Torres et al., 2015), which may imply that the behavior of granular column collapses is strongly influenced by the contact network presented inside the column during the collapse. Additionally, to better understand the scaling of granular column collapses, we further analyze the correlation length scale at the very beginning of the collapse. Simulation results show that, as we increase α eff −α c∞ , the length scale ξ/H t shows a power-law decay with the critical exponent ζ ≈ 0.89. This study is based on our previous work where we introduced a physics-based dimensionless number (ratio between inertial effects and frictional effects) to describe the behavior of granular column collapses, and we further expand our analysis to include the size effect, which is crucial to applications in engineering, such as chemical engineering, food processing, civil engineering. We also associate the size effect to the strong force networks at the beginning of the column collapse to show that a larger column size results in a well percolated strong force network, which may contribute to relatively larger run-out distances.
However, the granular system we are studying is still axisymmetric, and our preliminary results show intriguing phenomena when the cross-section of a granular column is no longer axisymmetric. Also, we note that the system we are dealing with has no interstitial fluid, and only Voronoi-based particles are considered. In natural and engineering applications, most granular systems are wet, saturated, or submerged in liquid. A rigorous investigation on the fluid-solid interaction and its influence on the bulk behavior, which is beyond the scope of this study, is need to thoroughly describe the behavior of granular-fluid systems. Further investigations related to granular-fluid systems, the influence of cross-section shapes, and the influence of particle shapes will be conducted and presented in future publications.

Data availability Statement
The data associated to this work is archived and published as Man (2021) in Mendeley Data, which can be accessed freely after registration. Part of the data in Fig. 3(g) is available through Calder et al. (1999). The simulation and the corresponding model were described in MechSys: Multi-physics Simulation Library (2021) and S. Galindo-Torres (2013). The software can be obtained from http://mechsys.nongnu.org/.