Application of an Inertia Dependent Flow Friction Model to Snow Avalanches: Exploration of the Model Using a Ping-Pong Ball Experiment

Abstract

Snow avalanches are catastrophic phenomena because of their destructive power. Therefore, it is very important to forecast the affected area of snow avalanches using numerical simulations. In our study, we focus on applying a numerical model to snow avalanches. The inertia-dependent flow friction model, which we call the “I-dependent” model, is a promising numerical model based on granular flow experiments and includes the local inertial effect. This model was introduced in previous studies as it predicts the shape and velocity of the granular flow accurately. We numerically investigated the particle diameter effect of the I-dependent model, and found that the smaller the particle diameter is, the faster the flow front velocity becomes. The final flow shape is similar to a crescent shape when the particle diameter is small. We applied this model to the ping-pong ball flow experiment, which imitated a snow avalanche on a ski jump slope. Comparing between the experimental and simulated results, the flow shape is better reproduced when the particle diameter is small, while the numerical simulation using a real ping-pong ball diameter did not show the clear crescent shape. Moreover, the relative error analysis shows that the best fit between experimental and simulated flow front velocity occurs when the particle diameter is larger than the actual size of a ping-pong ball. We conjecture that this discrepancy is mainly caused by aerodynamic effects, which, in this case, are large due to the low density of ping-pong balls. Therefore, it is necessary to explore the granular features of ping-pong balls or snow avalanches by conducting experiments, as done in previous experimental studies. Through such efforts, it may be possible to apply this I-dependent model to snow avalanches in the future.

1. Introduction

Snow avalanches are ubiquitous phenomena in snowy mountainous regions. They sometimes damage infrastructure and cause fatalities due to their energetic power. Therefore, it is very important to forecast the spreading and runout distance of snow avalanches using numerical simulations. We focus on developing a numerical model for snow avalanches that can predict the extent of avalanche effects.

Dry snow avalanches consist mainly of two parts. The first part is the basal dense layer where granular dynamics are dominant. The second part is the dilute powdery layer on top of the dense layer. This part is also called the suspension layer because particles are suspended in the air by turbulence. We focus on the granulometric characteristics of the dense layer where the basal and internal frictions are dominant because the dense layer is thought to be a major part of such avalanches, determining approximate deposit shapes.

Physical characteristics of snow avalanches were studied by observations (e.g., [1,2,3]) and experiments (e.g., [4,5]). Observational works revealed a correlation of run-out distance and terrain morphology [1], while other scientists studied snow entrainment and its effect on the mass and velocity of the avalanche flows [2,3]. Laboratory experiments can easily reproduce granular flows and investigate their characteristics, while it is difficult to observe individual particles [4]. Moreover, laboratory experiments cannot show the three-dimensional characteristics of avalanches [6]. Therefore, ping-pong ball flow experiments on a ski jump slope with 300,000 ping-pong balls were carried out (Figure 1a, [6,7,8]). The velocity distribution of the flows was analyzed and has revealed that the flow velocity increases with the number of ping-pong balls and the velocity is proportional to 1/6th of the power of the number of balls [7]. Moreover, it was found that the ping-pong ball flows resulted in a crescent shape (Figure 1b), like many snow avalanches in nature [9]. This experimental result implies that the crescent shapes of both the snow and ping-pong ball flows are the outcome of granular flow characteristics. Although the mass entrainments and depositions are important aspects of the snow avalanche, herein we focus more on the granular features of the snow avalanches for example, crescent shapes.

To simulate the characteristics of snow avalanches, such as the flow thickness, flow speed, and run-out distance, models have an ingenious friction law. RAMMS [10] is probably the best-known model among many numerical models for snow avalanches. This model is based on the Voellmy friction law [11]. Some other models based on the work by [12,13] use a Mohr–Coulomb-type failure criterion, characterized by a single bed friction angle for the bed friction law. The Savage–Hutter model is also incorporated into a numerical model called Titan2D [14]. In Titan2D, the Savage–Hutter model was extended to be used with an arbitrary 3D-terrain model. Titan2D is widely used for simulating debris flows, volcanic pyroclastic flows, and snow avalanches [15,16,17,18].

Although constructing a numerical model to capture major characteristics of granular flows is a challenge, some studies have successfully reproduced key morphological features such as levees or channel structures as observed in nature and the laboratory, using a depth-averaged model with a velocity-dependent friction law (e.g., [19,20]). Maeno et al. [21] applied an inertia-dependent flow friction model (I-dependent model) to study the slumping of a granular mass on a slope in the laboratory, and showed that the model can work well to reproduce these phenomena. The characteristic of the I-dependent model is that the friction is controlled by the inertia number I, and I itself varies with shear rate, particle diameter, and pressure, as explained in the next section. This type of friction model incorporates the parameters related to the particles and may be useful to simulate granular flows in various situations. Intriguingly, the simulation results of the I-dependent model reproduce the crescent shape, while Titan2D does not (Figure 2). This suggests that the I-dependent model could be useful to reproduce the ping-pong ball experiments as well. Originally, the model implements inertia-dependent flow friction laws derived from experimental works such as [19,22,23]. With this background, we have started to think about using the I-dependent model for simulating snow avalanches. Some encouraging examples were presented by [24] and a Swiss group [25,26]. They modeled the relation between shear and normal stresses. Nishimura and Maeno [24] measured cohesion based on their chute flow experiments using real snow particles. They proposed using the Bingham relation for shear and normal stresses among other shear–normal stress relations. Being inspired by their study, the Swiss group repeated experiments with their larger slope using both wet and dry snow particles. They also measured cohesion and suggested a new relation for the shear and normal stresses similar to the I-dependent model. Counting the Swiss shear–normal stress relation in the numerical simulation, they successfully reproduced the crescent shape feature of real snow avalanches [26]. Their outcomes encouraged us to work with the I-dependent model by applying it to real snow avalanches. Recently, Issler et al. [27] pointed out that I-dependent rheology has many features that are observed in the dense part of snow avalanches. With these encouraging examples, we conducted a comparison between our model results and the ping-pong ball avalanche results.

In this paper, we discuss the effects of particle diameter in the model which was not clearly discussed in [21] and compare it with the experimental results of the ping-pong ball experiment. We explain the I-dependent friction model and its numerical model in Section 2. The method of verifying the flow model and searching for the appropriate parameter values by comparison with the ping-pong ball experiment results is also shown in Section 3. The simulation results are shown in Section 3. The effect of particle diameter and the model characteristics are discussed in Section 4. Concluding remarks are presented in Section 5.

2. Model and Method

The numerical model of granular flow proposed by [21] featured a dense flow layer at the bottom of the granular flow. Their equations of mass and momentum conservation use a depth-averaged assumption with the latter conservation including a source term which contains the basal shear stress, T. The basal shear stress T includes the friction coefficient μ, which is a dimensionless parameter defined as the ratio of shear stress and normal stress ($\mu =\tau /\sigma$, $\tau :$ shear stress; $\sigma :$ normal stress). μ depends on the inertia number I. Therefore, we call the numerical model of [21] the “I-dependent model”, while other researchers refer to this model as the “$\mu \left(I\right)$” model (e.g., [23]). In the I-dependent model, the relationship between the friction coefficient μ, and the inertia number I was derived by [28] as follows,

$$\mu \left(I\right)={\mu }_S+\frac{{\mu }_2-{\mu }_S}{\frac{I_0}{I}+1}$$

(1)

where I₀ is a constant that is derived based on the experimental results. The value of I₀ for glass beads is approximately 0.3 according to [22], and we use this value even though our material is different. For the constant friction model, ${\mu }_2$ and ${\mu }_S$ are defined as

$${\mu }_2=tan{\delta }_2, {\mu }_S=tan{\delta }_S$$

(2)

where ${\delta }_2$ and ${\delta }_S$ are the upper and lower limits of slope angles in Figure 3 of [29] expressed as θ (°). In fact, the relationship of Equation (1) was derived based on the experimental results of [29] (the first researcher to conduct this experiment) and other following studies (e.g., [22,28,30]). Numerically, a granular flow can be treated as an incompressible and depth-averaged flow. Thus, the inertia I is expressed in its depth-averaged form $⟨I⟩$ using the depth-averaged velocity $\widehat{V}$, particle diameter d, gravity acceleration $g$, slope inclination angle θ, and flow thickness h as

$$⟨I⟩=\frac{5}{2h}\frac{\widehat{V}d}{\sqrt{ghcos\theta }}$$

(3)

The derivation of Equation (3) was explained in detail in [21,22]. The governing equations such as mass conservation and momentum conservation are based on the shallow water theory. We explain the governing equations by referencing Titan2D, which is also a shallow water model. The global equations are

$$\frac{\partial }{\partial t}\left(\overrightarrow{U}\right)+\frac{\partial }{\partial x}\left(\overrightarrow{F}\left(\overrightarrow{U}\right)\right)+\frac{\partial }{\partial y}\left(\overrightarrow{G}\left(\overrightarrow{U}\right)\right)=\overrightarrow{S}\left(\overrightarrow{U}\right)$$

(4)

where vectors $\overrightarrow{U}$, $\overrightarrow{F}\left(\overrightarrow{U}\right)$, and $\overrightarrow{G}\left(\overrightarrow{U}\right)$ are as follows:

$$\overrightarrow{U}=\left[\begin{array}{c}h\\ h{V}_x\\ h{V}_y\end{array}\right], \overrightarrow{F}\left(\overrightarrow{U}\right)=\left[\begin{array}{c}h{V}_x\\ h{V}_x^2+C_{0}g{h}^2\\ h{V}_x{V}_y\end{array}\right], \overrightarrow{G}\left(\overrightarrow{U}\right)=\left[\begin{array}{c}h{V}_y\\ h{V}_x{V}_y\\ h{V}_y^2+C_{0}g{h}^2\end{array}\right].$$

(5)

where V_x and V_y are the flow velocities in the x and y directions, respectively. The first equation is a mass conservation equation, and the second and the third equations are momentum conservation equations in the x and y directions. The x-axis represents a downslope direction on a plane parallel to the slope in our simple terrain similar to a ski jump hill (Figure 1). The y-axis represents the width of the slope. C₀ is a constant showing a clear difference between the I-dependent model and Titan2D, being defined as $C_0=1/2$ in the I-dependent model, and as $C_0=k_{ap}/2$ in Titan2D, where $k_{ap}$ denotes an earth pressure coefficient. Moreover, the source term of Equation (4) “$\overrightarrow{S}\left(\overrightarrow{U}\right)$” shows the features of each numerical model because this term includes friction effects. The source term of the I-dependent model is

$$\overrightarrow{S}\left(\overrightarrow{U}\right)=\left[\begin{array}{c}0\\ ghsin\theta -T_x/\rho \\ -T_y/\rho \end{array}\right]$$

(6)

with $T_x=\frac{\mu \rho gh cos\theta V_x}{\sqrt{V_x^2+V_y^2}}$ and $T_y=\frac{\mu \rho gh cos\theta V_y}{\sqrt{V_x^2+V_y^2}}$, where $\rho$ denotes the bulk density and $T_x$ and $T_y$ denote the basal shear stress components defined as T = ($T_x, T_y$). The source term of Titan2D is provided here as a reference.

$$\overrightarrow{S}\left(\overrightarrow{U}\right)=\left[\begin{array}{c}0\\ g_{x}h-h{k}_{ap}sign\left(\frac{\partial V_x}{\partial y}\right)\frac{\partial \left(gh\right)}{\partial y}sin{\varphi }_{int}-\frac{V_x}{\sqrt{V_x^2+V_y^2}}max\left(g+\frac{V_x^2}{r_x},0\right)htan{\varphi }_{bed}\\ g_{y}h-h{k}_{ap}sign\left(\frac{\partial V_y}{\partial x}\right)\frac{\partial \left(gh\right)}{\partial x}sin{\varphi }_{int}-\frac{V_y}{\sqrt{V_x^2+V_y^2}}max\left(g+\frac{V_y^2}{r_y},0\right)htan{\varphi }_{bed}\end{array}\right]$$

(7)

where $g_x$ and $g_y$ are gravity accelerations in the x and y directions, and ϕ_int and ϕ_bed are the internal and bed frictions. As seen in these equations, the source term of Titan2D clearly consists of the gravity term, internal friction term, and basal friction term, while the source term of the I-dependent model consists of the gravity term and shear stress terms of T_x and T_y, neglecting the internal friction term (see [31] for another I-dependent model that retains internal friction). This includes $\mu$ of Equation (1). According to these equations, ${\delta }_2$ and ${\delta }_s$ of the I-dependent model are physically different from ϕ_int and ϕ_bed of Titan2D. We need to be careful of this difference when we discuss these two models. The stopping criterion is applied using the threshold value of the shear stress ${\sigma }_c$ defined as

$${\sigma }_c=\mu \rho gh\cos \theta .$$

(8)

If the norm of the basal shear stress $\left|T\right|$ is smaller than ${\sigma }_c$, the velocity becomes 0.

Numerical simulations are implemented using the finite volume method developed by [21]. They were based on the code of [32], but with the source term calculation from the Runge–Kutta method.

The parameters for simulating snow avalanches were considered based on the experimental and observational results. Firstly, particles in real snow avalanche flows were measured by [33]. They measured the particle diameter of snow avalanche debris in snow deposits. Graphs of the particle diameter are shown in Figure 3. These histograms of wet and dry snow particles have a peak between 0 and 10 cm. The graph of the wet avalanche particles shows a wider variation and many particles are 10–20 cm in diameter. Snow avalanche particles increase in size [33] or are fractured by collision [34] during the flow. The observed snow in Figure 3 was both wet and dry, and the peaks of both particles are within the range of 0–10 cm. Moreover, the median of the grain size distributions, which are also taken from the field observation of the Swiss avalanche sites, was within 5–16 cm [35], and the order was similar to that measured by [33]. Therefore, we used a particle diameter from 0.1 mm to 10 cm as an input parameter value.

According to [21], their simulation results with δ_s = 22–26° and δ₂ = 28–32° were in agreement with the experimental results carried out using glass beads. Even though our experiments involved ping-pong balls or snow, and not glass beads, we used their values to define the input values of the friction angles δ_s = 20–26° and δ₂ = 28–34°, because the best estimates of the ping-pong balls or snow material are not known.

3. Results

Here, the numerical simulation results are shown: one investigating the particle diameter effect of the I-dependent model and the other a comparison with the ping-pong ball experiment [8], where the effects of friction angles were also considered.

3.1. Particle Diameter and Flow Characteristics

The diameter of snow avalanche particles is variable as shown in Figure 3. Therefore, we investigated the dependency of flow characteristics on the particle diameter.

Firstly, we have conducted numerical experiments using various size particles and slopes to examine how particle diameter d affects the flow shape. In these numerical experiments, we used slopes which have steep angles (47°, 37°, 26°, and 15°) at the releasing position, and thus, the release heights (corresponding to the “total vertical drop” in [9]) vary as 2190, 1640, 1090, and 550 m respectively. The steepness gradually decreases with distance from the releasing position, and finally, the slopes become flat (Figure 4). A cylindrical pile, which has a 50 m radius and a 10 m height, is set at the release position and the granular flow starts from rest. The model parameters are chosen as ${\delta }_s=20°$ and ${\delta }_2=30°$.

The final resting shape becomes elongated in the direction perpendicular to the flow direction and tapers downstream as the particle diameter becomes smaller (Figure 5). In other words, the smaller the particle diameter is, the clearer the crescent shape becomes visually. Our analysis of the width (elongation perpendicular to the flow direction) and tail length (stretch length in the flow direction) is shown in

Table 1. Width and tail length of the final resting shape of simulated granular flows. The width and tail length are measured using the simulation results shown in Figure 5.

	Particle Diameter
	10 cm	1 cm	1 mm	0.1 mm
Width (m)	525	555	585	615
Tail length (m)	155	345	525	580

. As shown in Table 1, both the width and tail length become larger as the particle diameter decreases.

The flow front positions are plotted against time for the two slopes, which have slope angles of 47° and 26° at the release position (Figure 6). The flow front velocity is faster when the particles are smaller. The comparison of flow front positions shows that the relationship between the velocity and particle diameter is maintained even when the slope inclinations and release heights change (Figure 6a,b).

3.2. Comparison with the Ping-Pong Ball Experiment

Secondly, numerical simulations were implemented in order to verify whether this model can reproduce results from ping-pong ball experiments. To simulate ping-pong ball flows on the ski jump slope, we used the slope data of [7] (Figure 7).

The height difference between the release position and the flat plane at the bottom of the slope was approximately 70 m. The slope inclination was 18° at the release position and 36° at the steepest part. In total, 28 simulations were run (four δ_s cases and seven δ₂ cases) as shown in

Table 2. Conditions of simulations for ping-pong ball flows.

	δ2 = 28°	δ2 = 29°	δ2 = 30°	δ2 = 31°	δ2 = 32°	δ2 = 33°	δ2 = 34°
δs = 20°	case 01	case 05	case 09	case 13	case 17	case 21	case 25
δs = 22°	case 02	case 06	case 10	case 14	case 18	case 22	case 26
δs = 24°	case 03	case 07	case 11	case 15	case 19	case 23	case 27
δs = 26°	case 04	case 08	case 12	case 16	case 20	case 24	case 28

. In our simulation, a cylindrical pile, which has a 5 m radius and a 5 m height, was released. In order to study the effect of particle diameter on our numerical model, we used five different sizes of particles—10 cm, 3.8 cm, 1 cm, 1 mm, and 0.1 mm—though the diameter of the real ping-pong ball used in the experiment was 3.8 cm. Therefore, we have carried out $4\times 7\times 5=$140 calculations for simulating the ping-pong ball experiment.

The final shapes of the ping-pong ball deposits for ${\delta }_{\mathrm{s}}$= 20° and ${\delta }_2$= 29° (case 05) are shown in Figure 8 as an example. This variation of the flow shape for different particle sizes is common for all other cases; a crescent shape is clearer for smaller particle sizes. In order to compare the simulated results and experimental results, we calculated the flow front velocity from the simulation as it was the only parameter that could be obtained from both ping-pong ball flows and simulated flows. To have an intuitive image of the difference between experimental flows and simulated flows, we show the graphs of the flow front velocity, including experimental and simulated results for ${\delta }_{2 }$ = 29° as an example (Figure 9). Our simulation results show that the flow is decelerated for the first moment and accelerated after deceleration, showing a smooth and convex parabolic curve. Plots of the raw data show a decrease in velocity around the first 10 m (Figure 9), although the smoothed curve ignores this reduction in velocity. Therefore, the qualitative feature of the velocity is the same among experimental and simulated results. These profiles also show that the velocity distributions of smaller particle flows are always larger than those of large particles, which agrees with the simulated results using the various particle diameters, as shown in Section 3.1.

Despite the similar qualitative features, the simulated velocities show a large difference from the experimental profiles with ${\delta }_2$= 29°. This discrepancy becomes smaller when both ${\delta }_2$ and ${\delta }_s$ are increased. As the simulated flow of smaller particles is faster than that of larger particles, the larger the particles become, the better the velocity fits with the experimental results.

When we focus on the real ping-pong ball size (3.8 cm), the velocity distributions closely resemble the experimental results in the case of δ₂ = 32° and δ₂ = 33° (Figure 10). The maximum velocity of the experimental results is especially similar to that in the simulated cases of {δ_s = 26°, δ₂ = 32°} and {δ_s = 20°, δ₂ = 33°}.

3.3. Error Analysis

To investigate the best fit and convergence conditions of the experimental and simulated velocity distributions, we carried out an error analysis using these velocity data. The error was calculated by interpolating the velocities with a distance resolution $\Delta x$= 1 m as a standard across all simulations. We applied spline interpolation to the simulated velocity. The relative error is calculated as

$$E_r=\frac{\sum \left|v_e^i-v_s^i\right|}{\sum v_e^i}$$

where $v_e$ is the velocity of the experimental result and $v_s$ is the velocity of the simulated results at each distance. The index $i$ denotes each position of the distance where the interpolated values were chosen with a resolution, $\Delta x$= 1 m.

The relative error of velocity for each particle diameter (10 cm, 3.8 cm, 1cm, 1 mm, and 0.1 mm) was plotted on the given ${\delta }_s$ and ${\delta }_2$ area (Figure 11). The sets of parameters which give the best fit for each particle diameter are {d, ${\delta }_s$, ${\delta }_2$, E_r} = {10 cm, 20°, 32°, 0.1486}, {3.8 cm, 24°, 32°, 0.1665}, {1 cm, 26°, 32°, 0.1860}, {1 mm, 26°, 34°, 0.2442}, {0.1 mm, 26°, 34°, 0.4080}. According to these numbers, we could see that the larger the particle diameter is, the better the velocity of the experimental and simulated results fit. In order to clarify how the minimum relative error of the velocity varies with the particle diameter, these minimum error values were plotted against the particle diameter (Figure 12).

The graph in Figure 12 is a log–log plot of the relative error against the particle diameter. We expected to have a convergence of the relative error of velocity, although it does not show a flattening of the curve at the bottom. This means that the error can be smaller if the particle diameter is larger and the convergence of the relative error does not appear within the calculated range.

4. Discussion

In the proceeding section, we have shown the results of the numerical simulations with the I-dependent friction model. Furthermore, we have compared the simulation results with the experimental result of a ping-pong ball avalanche experiment on a ski jump hill. Based on these results, we discuss the reproducibility of the shape and flow front velocity in terms of the effect of particle diameter d and two model parameters, the upper and lower limits of slope angles allowing steady flows,${ \delta }_2$ and ${\delta }_s$.

4.1. Shape Reproducibility

The flow shape obtained in the numerical experiments has the following tendency—the more elongated and trailed off the flow deposit shape becomes, the more clearly the crescent shape emerges (Figure 5). This is probably due to the inertia number definition $I=\frac{T_p}{T_{\gamma }}$, which is originally a ratio of the timescale $T_p$ and $T_{\gamma }$, where $T_p$ is the confinement timescale and $T_{\gamma }$ is the typical time of deformation [22]. Let us assume that there are two particle layers, one above the other. $T_{\gamma }$ is the travel time of distance d laterally within a moving flow, and the $T_p$ is the time for a particle to drop into the hole below. This is to say that all particles, whether piling up, resting, or proceeding, are controlled by the inertia number I. Then, whether or not this particle behavior produces a crescent shape depends on the particle diameter d. Based on the equations, the friction coefficient μ increases with the inertia number I (Equation (1)), and it increases with the particle size d as well (Equation (3)). The flow becomes more mobile as the particle size d decreases. The final resting shape width becomes larger and the tail becomes longer as the flow spreads.

The comparison of the flow deposit shape between the ping-pong ball experiment result and the simulation results showed that the flow deposit shape is clearly a crescent shape with the particle size d = 1 and 0.1 mm. However, when we used the true particle size of the ping-pong ball (d = 3.8 cm), the simulated results did not show the clear bending shape as observed in the experiment (Figure 1b). There might be a significant discrepancy between the numerical model and the ping-pong ball flow. We will discuss this point in more detail later (Section 4.2).

4.2. Velocity Reproducibility

The velocity of the simulation is also dependent on the particle size d. The smaller the particle size is, the faster the flow becomes (Figure 6) for fixed values of ${\delta }_s$ and ${\delta }_2$. Based on the simulations with various slope angles, this tendency is qualitatively the same even if the slope angle and the release height change (Figure 6a,b). The inertia number effect also appears in the relationship between the flow velocity and the particle diameter. Mobility or deformability is explained as the competition of particle movement going forward and dropping in the gap between granular flow, and thus, the flow becomes more mobile and faster when the particle size d decreases. Ultimately, the simulated flow velocity is higher than the ping-pong ball velocity when the particle size d is small. In order to find the best fit condition, we had to compensate for the reduced friction from smaller particles (due to the lower value of $⟨I⟩$ at a given velocity and flow depth) by increasing both ${\delta }_s$ and ${\delta }_2$. Therefore, in the cases where the particle size is d = 1 and 0.1 mm, the minimum relative error emerges at the largest value ${\delta }_s$ and ${\delta }_2$, which is shown at the top right cell in Figure 11d,e.

Drawing a graph of the relative error against the particle diameter (Figure 12), the best fit particle size (10 cm) is found to be larger than the real particle size of a ping-pong ball. There are several possible reasons for this discrepancy: (1) the ping-pong ball flow is not suitable for verifying dense granular flow behavior, (2) the rheology is ill-posed mathematically as [36] suggested, and (3) the value of the parameter I₀ for ping-pong ball flows differs substantially from 0.3 and needs to be calibrated or extracted from laboratory experiments with ping-pong balls.

An issue pointing towards possibility (1) as an explanation is the transition of the front of the ping-pong ball flow to a semi-dilute state, as ping-pong balls were jumping and colliding with each other during the descent. A picture of the ping-pong ball experiment (Figure 1a) shows that the flow front became semi-dilute and two “eyes-like” features appeared at the flow front. Aerodynamic drag on the light ping-pong balls becomes non-negligible at moderate velocities; consequently, the ping-pong ball flow is slower than the simulated dense flow with the same particle diameter.

The second possibility, that the rheology is ill-posed in a certain condition, was recently studied by [36]. They showed that the I-dependent friction model is mathematically ill-posed for sufficiently large or sufficiently small values of I/I_0. Some of our simulations are in the unstable region, but no signs of numerical instability have been observed (perhaps due to a sufficient degree of code-inherent numerical viscosity). Therefore, it appears unlikely that ill-posedness of the model equations can explain the mismatch between the best-fit particle diameter and the real particle diameter.

The third possibility, that the material difference needs new experiments to calibrate the model, should always be checked when we have new material. Values used in our simulations were based on the experimental result with glass beads. It would be most valuable to have directly determined values of the friction parameters ${\mu }_2$ and ${\mu }_s$ or the corresponding angles ${\delta }_2$ and ${\delta }_s$ (Figure 13) and the shearing scale I₀. Unfortunately, our attempts at measuring the angle ${\delta }_s$, which can be viewed as a repose angle, for ping-pong balls, have failed so far.

Obtaining values of these parameters for snow avalanches presents even larger challenges: snow deforms, melts, sticks, and makes larger grains [38]. It remains to be seen whether experiments similar to those of Pouliquen and Forterre [19] can be carried out with real snow. The chute experiments reported in [39] represent a first step in this direction.

Moreover, it will be important to understand the dynamics of granulation because the granulation changes the effective particle size of the flow. Further experimentation along the lines suggested by [34] should be most valuable in this respect.

We conclude from this discussion that the discrepancy between the real and best-fit particle sizes is most likely due to the transition of the flow front to a semi-dilute regime and the non-negligible effect of aerodynamic drag on the ping-pong balls. Further work is needed to estimate the extra retarding force of aerodynamic drag and to include it in the model as an extra source term.

Furthermore, it is essential to validate this model with observed data in order to apply this model to the realistic snow avalanches. Avalanche velocity was estimated based on the seismic signals (e.g., [18,40,41]), and the Doppler radar measurements [42,43,44]. We can use these observational results for the validation of the model. Gauer [45] revealed the scaling behavior of the front velocity, which is proportional to the square root of the fall height. It might be interesting to see if the I-dependent rheology can reproduce such velocity features of snow avalanches.

5. Conclusions

We have presented a numerical model using an I-dependent friction model [21], which was developed based on the experimental approaches of granular flows done by Pouliquen and collaborators (e.g., [19,29,37]). After conducting a parameter study with a variable particle diameter d, we clarified that the smaller the particle is, the more the granular flow becomes mobile, resulting in a more widely spreading faster flow.

The deposit shape of the flow is so dependent on the particle diameter d that the real ping-pong ball diameter could not reproduce the crescent or bending shape of the flow. Unexpectedly, the simulated flow front velocity with smaller particles does not agree with the experimental flow front velocity even though the crescent shape reproduced well with the smaller particles. This discrepancy was mainly caused by the difference in properties of the various materials. The ping-pong balls in the flow jumped and collided with each other, since most of them were affected by aerodynamic drag. Ping-pong balls have a different property than glass beads, but might be similar to the snow avalanche property, because snow avalanches are sometimes dilute, sometimes dense, and sometimes in an intermediate regime.

Considering the application of the I-dependent model to the snow avalanches, we should study the snow avalanche property in macro and micro scales. The mean particle diameter d in a snow avalanche changes during the event. Flow experiments using real snow are expected to provide more realistic values for the two critical slope parameters ${\delta }_2$ and ${\delta }_s$. Additionally, I₀ should be chosen on the basis of experiments, rather than using the value 0.3 appropriate for glass beads.

With these experiments and improvement of the numerical code instability, it may be possible to apply this model to practical simulations in order to predict the run-out distance and affected area of snow avalanches. For the next step, it will also be necessary to improve the usability of the model.