3D DEM investigation of granular column collapse: Evaluation of debris motion and its destructive power

tool quantitatively predicting the destructivepowerofagiven ﬂ owatany locationofinterestalongitspath.Thiscanbeusefulfor thedesignofen-gineering works for natural hazard mitigation. To this end, also the distribution of the linear momentum of the ﬂ ow over depth wascalculated. It emerges that the distribution isinitially bilinear, due to the presence of anup-permost layer of particles in an agitated loose state, but after some time becomes linear. This type of analysis showcases the potential of the Distinct Element Method to investigate the phenomenology of dry granular ﬂ ows and to gather unique information currently unachievable by experimentation. © 2014 Elsevier B.V. All rights reserved.


Introduction
Long run-out granular flows (e.g.rock and debris avalanches) can travel distances several times larger than the initial size of their source topography, sweeping away populated areas located far away from the landslide source (Crosta et al., 2005;Carrara et al., 2008).Many theories and assumptions have been proposed in the attempt to explain the apparent high mobility of granular material, including the air cushion trapped at the base of moving mass (Shreve, 1968), basal rock melting (Goren and Aharonov, 2007;De Blasio and Elverhøi, 2008), sand fluidization (Hungr and Evans, 2004), destabilization of loose granular material at the failure plane (Iverson et al., 2011), acoustic fluidization (Melosh, 1979;Collins and Melosh, 2003) or grain segregation-induced friction decrease (Phillips et al., 2006;Linares-Guerrero et al., 2007).However, no experimental evidence has been found to validate these theories and assumptions so far.
It has been recognised that debris flows and dry granular flows may behave similarly: for instance, they can sustain shear stresses with very slow deformation due to lasting, frictional grain contacts, and they can flow rapidly, sustaining inelastic grain collisions (Iverson, 1997).Thus, research has mainly focused on small scale laboratory experiments and numerical simulations of dry granular materials (Kerswell, 2005;Saucedo et al., 2005;Mangeney et al., 2007Mangeney et al., , 2010;;Roche et al., 2011).Although these studies make important simplifications of the problem, they are useful in elucidating the mechanical behaviour of granular flows under simple, well controlled conditions (Crosta et al., 2009).
With regard to the numerical simulations, the Distinct Element Method (DEM) (Cundall and Strack, 1979) has been widely used to simulate dry granular avalanches (Cleary and Campbell, 1993;Staron and Hinch, 2005;Lacaze et al., 2008;Utili and Nova, 2008;Utili and Crosta, 2011a,b).The application of DEM to the simulation of granular flows, as first proposed by Cleary and Campbell (1993) in 2D, proved to be useful to enhance our understanding of the behaviour of dry granular flows.
In the last two decades, it has become increasingly popular to shed light on fundamental mechanical characteristics of landslides (Staron and Hinch, 2007;Tang et al., 2009).As stated by Zenit (2005), the use of DEM to simulate landslides is very powerful, since all the numerical data are accessible at any stage of the simulation, including quantities which are difficult, or impossible, to be obtained directly from laboratory experiments, such as the distribution of energy and momentum inside the granular flow and their variation over time and space.
In this paper, we illustrate the results of 3D DEM analyses in plane strain conditions of granular flows originated by the quick release of prismatic granular steps.The paper is organised as follows: in Section 2 the DEM contact model is illustrated; in Section 3 a dimensional analysis of the problem is carried out; then in Section 4 the numerical results obtained are presented.The evolution of the granular column profile over time for various column aspect ratios and intergranular friction is presented together with the illustration of the flux of the kinetic energy and of the linear momentum of the granular system through vertical planes located at various distances from the initial position of the column and the distribution of linear momentum within the granular flow.Finally, in Section 5, conclusions on the capability of the DEM to model granular flows are illustrated.

DEM model and input parameters
The DEM open source code ESyS-Particle (Weatherley et al., 2011) was employed to run the simulations here presented.Preliminary runs of granular column collapse tests were run employing non-linear elasticity according to the Mindlin no-slip solution (Johnson, 1985) and linear-elastic contacts.No significant difference was found between these analyses.Hence, since linear elastic contacts require less computational time, linear elasticity was adopted.Accordingly, the normal contact force at the contact plane is linearly proportional to the relative normal displacement between two particles (see Fig. 1(a)), expressed as: in which K n is the normal contact stiffness and U n is the normal relative displacement between two spheres in contact.The response along the tangential direction (see Fig. 1(b)) is calculated incrementally, as: in which F t n and F t n − 1 are the tangential forces calculated at the current and previous simulation time steps; K s is the shear stiffness, and dU s is the incremental tangential displacement.The maximum tangential force is limited by the Mohr-Coulomb criterion (see Figure 1

(b)).
A moment-relative rotation law is also present.Here, the rolling law was employed with the only aim of accounting for the shape effect of non-spherical particles.The law accounts for moments arising from the fact that the line of action of the normal contact force in the case of non-spherical particles no longer passes through the centre of mass of the particles and hence generates rotational moments (Belheine  The particle size distribution adopted in the numerical simulations is plotted as the red curve.2009)).Without this law, unrealistically low angles of repose would be obtained for the granular assembly.In the model adopted, the magnitude of the elastic rolling moment (M e ) is proportional to the relative rotational angle (see Figure 1(c)), and is calculated incrementally as: where M e n and M e n − 1 are rolling moment values calculated at the current and previous simulation time steps; k r = βK s r 2 is the rolling stiffness, with β being the coefficient of rolling stiffness, r being the average particle radius; Δθ r is the relative rotational angle between two particles in one iteration step.The maximum rolling moment that can be exchanged is defined as: in which η is the so-called coefficient of plastic moment.
The particle size distribution (PSD) is one of the most important factors controlling landslide initiation and soil permeability.The PSD of granular flows varies hugely at different locations (see for instance Casagli et al. (2003)).In addition, the grain size distribution may vary significantly within the same landslide mass at different depths (Crosta et al., 2007).Fig. 2 shows examples of particle size distributions from 7 cases of landslides in the Northern Apennines (Casagli et al., 2003) and 6 cases of rock avalanches in Val Pola in the Alps (Crosta et al., 2007).It can be observed that the grain size ranges from 0.001 mm to 1000 mm, with a large percentage of fine and medium sized grains and a small amount of coarse grains.Large discrepancies can be observed between the various site investigations.
According to Fig. 2, grains with diameters ranging from 0.1 to 10 mm were widely observed in different locations.However, in DEM simulations, due to computational limitations, we used a much narrower particle size distribution with the ratio of maximum to minimum particle sizes equal to 2, as shown in the red curve in Fig. 2. The input parameters are listed in Table 1.
In this study, we carried out numerical experiments of granular columns analogous to the columns of the experimental investigation of Lube et al. (2005).Initially particles were generated within a prismatic space bounded by two rigid and frictionless walls (A and B in Figure 3(a)), a rough base (D in Figure 3(a)), and a periodic boundary in the out-of-plane direction (C in Figure 3(a)) to impose plane strain conditions.The horizontal base (D) of the domain (see Figure 3(a)), is made of particles of the same PSD as the granular column that were kept fixed at all times to simulate a non-erodible base of the same roughness as the flowing material.To generate the granular column, first particles were randomly created in space (with the grain position allocated according to standard algorithms of random number generation), then they were left to settle under gravity until a dense column was obtained.In all the simulations, the granular flow was initiated by the instantaneous removal of wall B.
Once the flow has stopped, the final run-out length (L f ) and deposit height (H f ) can be measured (see Figure 3(b)).The presence of single Fig. 4. The angle of repose of a granular deposit for different rolling resistances.
particles detached from the front of the granular mass sometimes makes the calculation of L f a non-straightforward exercise.Zenit (2005) calculates L f considering only particles remaining in contact with each other, disregarding individual loose particles detached from the deposit mass centre.In this study, in order to track the position of the front in a way consistent amongst the various simulations (the number of loose single particles moving ahead of the front depends on the aspect ratio), we implemented an algorithm which identifies the front as the boundary between 99% and 1% of the flow mass, i.e. 1% of the mass is travelling ahead of the boundary.This guarantees that the front is not confused with the position of loose single particles jumping ahead the flow.

Dimensional analysis
The physical parameters ruling the motion and the depositional morphology of the granular column include (Zhao et al., 2012): grain properties, geometrical properties of the domain, flow velocity and duration time.In Table 2, all the parameters are listed together with their units expressed in terms of fundamental dimensions, i.e. mass (M), length (L) and time (T).Depending on their role in the simulation, parameters can be categorised as either independent or dependent.We considered the grain properties and the domain geometry as independent input parameters and final runout distance (L f ), deposit height (H f ), flow velocity (v f ) and flow duration time (t) as the dependent ones.
The governing parameters of the mechanical laws ruling particle interactions are the normal and shear particle contact stiffness (K n , K s ), particle friction (ϕ μ ) and the coefficients of rolling resistance (β, η).Also porosity has a significant influence on the mechanical behaviour of granular flows (Craig, 1997).However, since the primary purpose of this study is to explore the capability of DEM in modelling granular flows, the influence of porosity has not been considered here.The relationship between the independent and dependent quantities can be expressed by a general functional relationship of the form: Performing dimensional analysis (Palmer, 2008), and taking out porosity which is the same in all the tests carried out (all columns were generated with the same initial porosity), Eq. ( 5) can be rewritten as: L i is the normalised maximum final deposit height; V ½ ¼ v f = ffiffiffiffiffiffiffi ffi gH i p is the normalised flow velocity (v f can be chosen in a variety of ways, one of which is the velocity of propagation of the flow front, i.e. v f = dL/dt); ½T ¼ t= ffiffiffiffiffiffiffiffiffiffi H i =g p is the normalised flow duration time; a = H i /L i is the initial column aspect ratio; ε = ρ s gH i /(K n /D) is a characteristic compressive strain of the granular column and [S] = H i /D is the model-to-particle size ratio.Note that [-] indicates a dimensionless variable.

Calibration of the angle of repose
In studies on non-cohesive granular flows, the material internal friction angle (ϕ) is often approximated by the angle of repose.In simple terms, it can be said that a slope with an angle over the horizontal larger than the angle of repose is unstable, whereas it is stable if the opposite holds true.Typical values for the angle of repose range from 25°for flows made of smooth spherical particles to 40°for rough angular particles (Carrigy, 1970;Pohlman et al., 2006).Here, we assumed as baseline reference value, an angle of repose of 31°which is representative of coarse quartz sand (Lube et al. (2005)).
With regard to the calibration of the DEM parameters, concerning the rolling stiffness, β, in Modenese (2013) it is shown that its influence on the angle of repose is negligible, so it was not investigated here.Then, two contact parameters remain to be determined: the intergranular friction angle, ϕ μ , and the coefficient of rolling resistance, η.In Zhao (2014), it is shown that some arbitrariness exists since the same angle of repose can be obtained from several combinations of them.Here, we opted to choose an intergranular friction angle of ϕ μ = 30°, which is a typical value for quartz sand grains, and to calibrate η against the angle of repose observed when the flow comes to a stop.In Fig. 4, the angles of repose obtained for various values of η are shown.For η N 0, a linear relationship between η and the angle of repose exists.For η = 0, an unrealistically low angle of repose is obtained.This is a wellknown result reported in Calvetti and Nova (2004), Pöschel and Buchholtz (1993) and Rothenburg and Bathurst (1992)).In all the subsequent analyses, a value of η = 0.1, corresponding to an angle of repose for the flow of 31°, was adopted.

Kinematics of motion
Depending on the initial aspect ratio of the granular column, different depositional morphologies are obtained (see Figure 5).As it can be expected, the upper surface of the granular flow gradually lowers until the angle of repose is reached.Our results confirm that the flow profiles at various times and the final deposit depend strongly on the initial aspect ratio (Lube et al., 2005).This highlights the need for presenting profiles in dimensionless form.To this end, the final profiles of the granular assembly, normalised by the final runout distance and the deposit height, were plotted in Fig. 6 together with the experimental profiles obtained by Lube et al. (2005).The profiles from our 3D DEM analyses exhibit a good agreement with the experimental results for all the aspect ratios analysed.
A key feature of the DEM is the possibility of obtaining data on the motion of single particles (e.g.trajectories, velocities, etc.) inside the granular flow, allowing performing a detailed analysis of the potential heterogeneities in the flow.For instance, from our analyses, it emerges that in granular columns with small aspect ratios (e.g. a = 0.93), only a small portion of material is involved in the motion.To compare the kinematic fields obtained from columns of different aspect ratios, it is convenient to plot the measured velocities in dimensionless form.To this end, we define ffiffiffiffiffiffiffi ffi gH i p as the characteristic velocity of the system.This represents the maximum velocity that a single particle in free-fall from a height corresponding to the centre of the column (H i /2) reaches at the end of the fall.In Fig. 7(a), particles with velocities smaller than 1% of the characteristic velocity are plotted in grey, whilst particles moving at higher velocities are plotted in red.From the figure, it can be noted that failure of the granular assembly occurs approximately along a plane, identified by a dashed line in the figure.According to the Rankine's theory of earth pressure (Rankine (1857), assuming that the failure of the granular column occurs when conditions of active thrust are in place, the inclination of the failure plane to the horizontal (θ f ) can be estimated as: where ϕ is the internal friction angle of the granular material.From Fig. 7(b), it emerges that the inclination angle of the active failure plane is approximately 61°for all the aspect ratios tested.This value is very close to the theoretical value predicted by Eq. ( 7), which is 60.7°.

Influence of the column aspect ratio
For rock avalanches and granular flows, the assessment of the final run-out distance is of primary importance, as it determines the extent of the regions affected by the avalanche or landslide.In Fig. 8(a), our simulations have been compared with previous numerical simulations (Staron and Hinch, 2005;Zenit, 2005;Crosta et al., 2009) and experimental observations (Balmforth and Kerswell, 2005;Lajeunesse et al., 2005;Lube et al., 2005) from the literature.We performed three sets of simulations: T1 and T2 refer to simulations run for the same values of ε (ranging from 6.5 × 10 −7 to 6.5 × 10 −6 ), but with η = 0.1 and η = 0 respectively, whilst T3 refers to simulations run for ε ranging from 6.5 × 10 −5 to 6.5 × 10 −4 and η = 0.1.Some of the T1 and T2 simulations were run with gravitational acceleration scaled up 100 times and K n = 3 × 10 7 N/m, whereas some others with unscaled gravity but particles 100 times softer, i.e.K n = 3 × 10 5 N/m.In fact, in light of dimensional analysis, scaling can be introduced by scaling either the value of gravitational acceleration or particle stiffness.The fact that the obtained results are aligned in consistent trends in Fig. 8 can be seen as a verification of the correctness of the performed dimensional analysis.Comparison between T1 and T2 simulations allows examining the influence of particle shape on run-out, whereas comparison between T1 and T3 allows examining the influence of ε.The values of the independent variables employed in our simulations are reported in Table 3 in terms of the dimensionless groups identified in Eq. ( 6).In the table, the values of the parameters employed in other tests reported from the literature in Figs. 8 and 9 for comparison purposes are listed as well.In Fig. 8(a), it can be observed that the final normalised run-out distance obtained from our simulations matches well from a qualitative viewpoint both the experimental observations of Lube et al. (2005) from tests run in plane strain conditions, and the 2D FEM numerical analyses of Crosta et al. (2009).Also it emerges that if rolling resistance is not employed (simulations T2), unrealistically large run-out distances are obtained since particle angularity tends to reduce run-out.In other words, a flow of spherical particles is more prone to sliding than a flow of particles of any non-spherical shape.Equally, if 2D DEM simulations are employed (Staron and Hinch, 2005;Zenit, 2005), unrealistically long run outs are obtained.Presumably, this is due to the fact that the 2D kinematics of particle interaction is too different from the real 3D kinematics.Comparison between the simulation series T1 and T3 shows that the characteristic strain of a granular column, ε, has no influence on the observed flow behaviour, at least for the range of values here employed.In conclusion, simulations T1 and T3 show that if particle shape is accounted for, albeit by means of a very crude approximation (i.e.employing a moment-relative rotation law with spherical particles which avoids simulating the real non-spherical geometry of the particles), the obtained run-out and final heights of the simulated flows are in good agreement with the available experimental data.Now let us look at how the front of the flow propagates over time.In Fig. 9, the position of the front is plotted against time for various column aspect ratios.According to previous literature (Lube et al., 2005;Crosta  analyses accounting for the effect of particle shape provide results that are in agreement with experiments.

Analysis of the energy contributions in the flow
The potential energy of the column at any time is: where m i and h i are the mass and height of a single particle i, respectively, and N is the total number of particles of the column.The kinetic energy of the system at any time is calculated as: where v i and ω i are the translational and angular velocities respectively of a generic particle i, and I is its moment of inertia (for a spherical particle I = 2mR 2 /5).
A part of potential energy gets dissipated rather than being transformed into kinetic energy, due to unelastic particle collisions (e.g.unelastic particle rebounds and frictional sliding).In light of the principle of energy conservation, the energy dissipated in the flow at any given time can be calculated as: where E 0 is the total energy of the system, which can be calculated from the initial potential energy of the column before particles start to move, as: where At the particle level, energy is mainly dissipated via frictional sliding at particle contacts and viscous damping along both the normal and tangential directions of contacts when damping is present, but also by the relative rotation between particles once the plastic limit of the rolling moment is reached (see Figure 1(c)).In Fig. 10, the entire temporal evolution of the energy components of the flows, from start of column collapse until end of motion, is plotted.After the instantaneous removal of the confining gate at T = 0, particles start to fall downwards, with potential energy being progressively transformed into kinetic energy.The kinetic energy exhibits a peak at around T = 1.0 (see the dashed line in the figure), when both the rate of potential energy loss and the rate of cumulative energy dissipation (see the slopes of the curves in the figure) reach their maximum values.After this time, both potential and kinetic energy decrease.The flow comes to a stop at about T = 4.0.Considering now columns of different aspect ratios, the dissipated  energy in terms of percentage of the initial total energy of the columns has been plotted against their aspect ratio in Fig. 11.From our simulations, it emerges that the higher the aspect ratio, the larger the proportion of energy dissipated during the flow.
Then, we investigated how much of the kinetic energy of the flow is due to the translation of particles and how much is due to their rotations.In Fig. 12, the two sources of kinetic energy are plotted separately against time.Unfortunately, viscous damping coefficients strongly depend on the type of material making the granular flow.Furthermore, a reliable experimental determination of damping coefficients for granular flows is currently not yet available.Therefore, we decided to run two cases only of simulations: one without viscous damping and the other one with a viscous damping coefficient of 0.3, in order to obtain a rough indication of the potential influence of viscous damping on the variables of interest but without trying to model any specific natural granular flow.In the simulations with viscous damping, the damping coefficient of 0.3 was employed in both the normal and tangential directions of particle contacts.This value corresponds to a coefficient of restitution of around 0.5 (Tsuji et al., 1992).In Fig. 12, it can be noted that the kinetic energy stemming from particle rotations remains negligible at all times (i.e. less than 0.5% of the total kinetic energy) in both cases, with and without the presence of viscous damping.In light of this finding, in the next section, where we investigate the spatial distribution of the momentum inside the flow, we have concentrated our attention on the linear momentum since we know the angular momentum to be negligible.Also, from Fig. 12, it emerges that the curves obtained for flows with and without damping, are similar.This implies that the presence of damping decreases the magnitude of the kinetic energy of the system, of the same proportion, throughout the duration of the flow.

Linear momentum
The evolution over time of the linear momentum of the flow, p ! ¼ , is plotted in Fig. 13, in the case of the presence of viscous dissipation and in its absence.Given the imposition of periodic boundary conditions in the y direction, we expect p y to be negligible at all times, as it is shown in Fig. 13(a).This confirms the effective presence of plane strain conditions in the x-z plane for the simulated flows.Analogously to the energy analysis of the flow, it is convenient to normalise p ! , for the sake of generality in presenting our results.To this end, we introduce a scalar quantity, p 0 , with: This quantity can be thought of as an approximate average of the momentum of the flow: being the characteristic velocity of the flow.In Fig. 13(a), the components of the linear momentum along the x, y and z axes, normalised by p 0 , are plotted against dimensionless time.A small difference between the curves for the case with and without damping is noted with damping having the effect of reducing the amount of momentum as it is expected.From the figure, it emerges that the momentum in the vertical direction, p z , exhibits a higher peak than the momentum in the direction of flow propagation, p x .The two peaks occur at different times.To better investigate which one is dominant and when, it is convenient to make a relative comparison between the two components.To this end, the square of each component over the square of the magnitude of the vector ( p ! ) is plotted in Fig. 13(b) and (c).The use of the squares allows for plotting the components as percent, since From the figure, it can be noted that at the beginning of the flow, the vertical component, p z , dominates due to the fact that the motion of the particles is mainly gravity driven free fall.Then, during the propagation phase, the horizontal component, p x , becomes dominant, stabilising itself at around 90%. Finally, when the flow is coming to rest, a surge of vertical component appears.This is due to the presence of decelerating particles exhibiting bouncing in the vertical direction, especially near the flow forefront.In comparison with Fig. 13(b), the start point of the chart in Fig. 13(c) looks shifted ahead in time, since it takes some time for the flow front to start moving ahead after gate removal.
Let us now examine the influence of the column aspect ratio on linear momentum.This is an important aspect in order to ascertain how general the observed trends on the linear momentum of the flow are.In Fig. 14, the evolution of the vertical and horizontal components of the linear momentum is plotted (as percentage over the magnitude of the vector p ! ) against dimensionless time for columns of various aspect ratios.Similar trends amongst the curves can be noted.However, in the initial phase of the flow (for T b 1.1), the vertical component of momentum increases with the aspect ratio and obviously the opposite is true for the horizontal one.Instead, after T = 1.1, the vertical component  of momentum decreases with the aspect ratio.A possible explanation for the observed aspect ratio dependent trend could be that the gravity driven free fall of particles, which gives rise to the particle vertical motion, increases with the height of the column whereas the friction between particles, which opposes the horizontal motion of the particles, is independent of the column aspect ratio.In Fig. 14(b) the components of the momentum are plotted against run-out distance.

Flux of kinetic energy
To assess the vulnerability of existing structures hit by a granular flow/avalanche and to design engineering works for the protection of existing structures, two quantities are of interest: energy and momentum.The kinetic energy of the particle flow can be seen as an upper bound on the destructive energy that could be unleashed on the structure impacted by the flow.The amount of kinetic energy transferred from the flow to the structure depends on flow and structure interact during the time the structure is impacted by the flow.Hence, the amount of energy released by the flow on the structure is a function of the characteristics of both flow and structure (for instance the relative stiffness between the two).Also the flow-structure interaction is likely to change over time, for instance due to the development of irrecoverable (plastic) deformations in the structure.So, it is not possible to predict the transfer of energy (and equally of momentum), unless a specific flow and a structure of interest are modelled.Here, however, we provide an analysis of the linear momentum and kinetic energy of the flow, measured at various distances from the initial position of the columns, in order to identify an upper bound on the maximum energy that may be imparted to structures knocked by granular flows, under the simplifying assumption of disregarding the effects of any structureflow interactions.
Considering an imaginary vertical section perpendicular to the direction of horizontal propagation of the flow (see the vertical plane depicted in Figure 15) the flow mass transiting through such a section is a function of time.This is nil until the flow front reaches the position, then it gradually increases, then decreases and eventually becomes nil again once the flow has come to a stop (see Figure 15).In the following analysis, five locations along the flow path, shown in Fig. 16, are considered.Each location is identified by a letter (see Figure 16).A convenient measure of the maximum energy that can be transferred from the flow to the impacted structure is what we here call the flux of kinetic energy.This is evaluated as follows: for a given location of interest, the kinetic energy of the particles passing through the vertical plane oriented in the direction perpendicular to the flow is recorded during a specified time interval Δt and their total kinetic energy, ΔE, is evaluated.The flux of kinetic energy through the plane is given by the ratio, ΔE/Δt, a quantity with the dimensions of power.It is convenient to normalise this quantity, so that a comparison of fluxes between columns of various sizes can be carried out.A simple way of doing so is by normalising both numerator, ΔE, and denominator, Δt.Thus, we define the normalised flux, P, as: Substituting E 0 from Eq. ( 11) into Eq.( 13) and rearranging, Eq. ( 13) can be written as: Eq. ( 14) makes clear that 2 M ffiffiffiffiffiffiffi g 3 H i p is the multiplying factor to normalise the measured flux ΔE/Δt.The flux, P, represents the maximum energy that can be transferred from the flow to the impacted structure.In fact, if all the energy were to be transferred away from the flow, the flow would be suddenly deprived of all its kinetic energy and therefore it would come to a stop, which is evidently an unrealistic scenario.In reality, only a portion of the energy of the flow is lost in the interaction with the structure that will cause the flow to slow down but not to stop.So, P can be thought of, as an upper bound on the maximum destructive power that the flow may impart on the impacted structure.In Fig. 17(a), the flux of kinetic energy at the selected locations is plotted versus dimensionless time.It can be noted that the section with highest flux is B. Obviously the fluxes in the case damping is present are smaller than the case of undamped flow.An interesting finding is the fact that the peak takes place at a different time, with the time of peak for the damped system shifting progressively ahead of the peak time for the undamped one.Also, the difference of value between the peaks for the damped and undamped systems increases with the distance of the section investigated from the column initial position reaching up to 50% of the peak value for the undamped system.Fig. 17(b) illustrates the evolution of the height of the granular mass at different locations.From the figure, it can be concluded that the further the location is away from the slope source region, the lower the height of the final granular mass is.
Fig. 18 illustrates the evolution of the flux of kinetic energy for different column aspect ratios.From the figure, it can be observed that the position of the section where the flux is highest depends on the aspect ratio.For instance, in the case of small aspect ratios (e.g. a = 0.93), the flux of destructive energy at location A is largest, with only a small amount of particles travelling to locations further down section B. For intermediate aspect ratios (e.g. a = 3.26, 5.91), the largest flux takes place at location B. For large aspect ratios (e.g. a = 9.27), the largest flux occurs again at section A.

Distribution of flux of kinetic energy and linear momentum over depth
To able to design protective structures as effectively as possible, the spatial distributions of kinetic energy and momentum over the depth of the considered sections are also needed.Considering section B, we have split the depth of the flow, h(t), into 5 parts and calculated the flux of kinetic energy through the section for each part so as to obtain a vertical profile of the flux of kinetic energy (see Figure 19(a)).Looking at the figure, it emerges that the profile of the flux is initially  bilinear, with the maximum flux at middle height of the column, then the flux evolves into a linear profile whose amplitude progressively reduces over time until becoming nil.The bilinear distribution points out to the presence of an uppermost layer of particles in an agitated loose state, which, after some time, consolidates so that a linear distribution is obtained.The presence of this layer of loose material is confirmed by the calculation of the profile of mass rate normalised by the total flow mass in the section (see Figure 19(c)).
A convenient measure of the maximum momentum that could be transferred by the flow to the impacted structure is the flux of linear momentum over time.Analogously to the flux of kinetic energy, we can evaluate the linear momentum of the particles passing through the vertical section of interest, Δp, during a specified time interval Δt.The flux of linear momentum through the whole section or parts of it, is given by Δp/Δt, a quantity with the dimensions of force (so we call it F).It is convenient to normalise this quantity so that comparison of fluxes between columns of various sizes can be carried out.A simple way of doing this is by normalising both numerator and denominator.So, we define the normalised linear momentum of the flow in the unit of time as: with p 0 being an average linear momentum for the flow here defined in Eq. ( 12).Substituting p 0 from Eq. ( 12) and rearranging, Eq. ( 15) can also be written as: Eq. ( 16) makes clear that 1/Mg is the multiplying factor to normalise the momentum going through the plane in the unit of time.In terms of the value of the flux of the momentum at various sections in time, similar figures as those obtained for the flux of kinetic energy (Figures 17 and  18) are obtained which are not reported here for sake of brevity.To locate critical sections and times of interest one of the two fluxes, either the flux of kinetic energy or of momentum, is enough.However, it is of interest to practitioners appointed with designing engineering works for the mitigation of the flow hazard, to know the value of the flux of momentum over depth in order to have an idea of the distribution of the pressure that can act on the structure.In Fig. 19(b), the distribution of the horizontal component of the linear momentum along the direction of flow propagation, F x , over depth is plotted against time.As it can be expected, the same shape of the profile as the profile of the flux of kinetic energy is found.

Conclusions
This paper presents a numerical investigation of the behaviour of dry granular flows generated by the collapse of prismatic columns via 3D Distinct Element Method (DEM) simulations under plane strain conditions.This type of analysis showcases the potential of the Distinct Element Method to investigate the phenomenology of dry granular flows and to gather unique information currently unachievable by experimentation.By means of dimensional analysis, the governing parameters of the problem were identified.Then, the influence of key variables of the problem was analysed.The main results are summarised as below: (1) Different regimes of granular motion have been observed, depending on the initial aspect ratios.The DEM results qualitatively match the FEM analyses by Crosta et al. (2009) and the experimental results by Lube et al. (2005).The granular material slides along a plane which approximately corresponds to the active failure plane of the column in agreement with Rankine's theory.
(2) Quantitative relationships between the column aspect ratio and normalised runout distance and deposit height were established.
Using the rolling resistance model, i.e. assigning a momentrelative rotation contact law, the DEM simulations give rise to runout distances and deposit height which match well the available numerical and experimental results.(3) A detailed analysis of how energy is dissipated by granular flows was performed from which emerges that most of the energy of the columns is dissipated by inter-particle friction, with frictional dissipation increasing with the column aspect ratio.Also, the translational and rotational components of the kinetic energy of the flows, associated to particle rotational and translational motions respectively, were monitored during the simulations.It is found that the rotational component is negligible in comparison with the translational one; hence in order to calculate the destructive power of a granular flow slide, only the translational contribution of the kinetic energy is relevant.(4) A methodology is presented to calculate the flux of kinetic energy over time carried by the granular flow through any vertical section of interest.This can be related to the energy released by landslide induced granular flows impacting against engineering structures under the simplifying assumption of neglecting all structure-flow interactions.This represents the first step towards achieving a computational tool quantitatively predicting the destructive power of a given flow at any location of interest along its path.This could be useful for the design of engineering works for natural hazard mitigation.To this end, also the distribution of the linear momentum of the flow over depth was calculated.It emerges that the distribution is initially bilinear, due to the presence of an uppermost layer of particles in an agitated loose state, but after some time becomes linear.

Fig. 2 .
Fig. 2. Particle size distribution of landslides occurred in Northern Apennines (Italy) (after Casagli et al. (2003)) and rock avalanches in the Alps (Italy) (after Crosta et al. (2007)).The particle size distribution adopted in the numerical simulations is plotted as the red curve.

Fig. 3 .
Fig. 3. Model configuration: (a) initial sample; (b) final deposit (A: fixed smooth back wall; B: movable front gate; C: periodic boundary; D: fixed coarse floor; L i : initial column length; H i : initial column height; L f : runout length; H f : deposit height).

Fig. 7 .
Fig. 7. Active failure state of the granular sample.(a) Motion of solid grains.(b) Inclination angle of the slope failure plane.

Fig. 8 .
Fig. 8. Relationship between aspect ratio and (a) normalised run-out distance and (b) final deposit height.The red symbols are experimental results; blue symbols are numerical results, while the numerical results of this study are coloured black.

Fig. 9 .
Fig. 9. Normalized granular spreading length versus normalized time from experimental and numerical tests.

Fig. 11 .
Fig. 11.Total energy dissipated during the flow normalized by the initial potential energy versus column aspect ratio.

Fig. 12 .
Fig. 12. Evolution of kinetic energy over time (a = 3.26) (E k_trans and E k_rot are the translational and rotational sources respectively of the kinetic energy of the system).

Fig. 13 .
Fig. 13.Evolution of linear momentum over normalised time for flow with and without damping: a) components of momentum normalised by p 0 ; b) components of momentum as a percentage of the total momentum against normalised time; c) components of momentum as a percentage of the total momentum against normalised run-out; x is the direction of flow propagation, y is the out of plane direction, and z is the vertical direction.

Fig. 15 .
Fig. 15.Schematic view of granular flows past a structure.

Fig. 17 .
Fig. 17.Evolution of the flux of destructive energy and height of debris mass.(a) Flux of kinetic energy at different locations (a = 3.26).(b) Height of flow at different locations (a = 3.26).

Fig. 18 .
Fig. 18.Evolution of the flux of destructive energy for different granular columns.

Fig. 19 .
Fig. 19.Profile of kinetic energy, momentum and mass trough plane B at different times (a = 3.26).(a) Profile of kinetic energy flux (b) Profile of normalised momentum.(Δp x is the rate of momentum through the portion of section considered; (c) Profile of granular mass moving (Δm is the mass flowing through the portion of section considered.

Table 1
Parameters of the granular flow simulations.

Table 2
Parameters of the granular flow problem.

Table 3
Properties for experimental and numerical tests.