A spectral multidomain penalty method solver for the numerical simulation of granular avalanches Método penalizado de multidominios espectrales para la simulación numérica de avalanchas granulares

How to cite item Trujillo-Vela, M. G., Escobar-Vargas, J. A., & Ramos-Cañón, A. M. (2019). A spectral multidomain penalty method solver for the numerical simulation of granular avalanches. Earth Sciences Research Journal, 23(4), 317329. DOI: https://doi.org/10.15446/esrj. v23n4.77683 This work presents a high-order element-based numerical simulation of an experimental granular avalanche, in order to assess the potential of these spectral techniques to handle geophysical conservation laws. The spatial discretization of these equations was developed via the spectral multidomain penalty method (SMPM). The temporal terms were discretized using a strong-stability preserving Runge-Kutta method. Stability of the numerical scheme is ensured with the use of a spectral filter and a constant or regularized lateral earth pressure coefficient. The test case is a granular avalanche that is generated in a small-scale rectangular flume with a topographical gradient. A grid independence test was performed to clarify the order of the error in the mass conservation produced by the treatments here implemented. The numerical predictions of the granular avalanches are compared with experimental measurements performed by Denlinger and Iverson (2001). Furthermore, the boundary conditions and parameters such as lateral earth pressure coefficients and the momentum correction factor were analyzed in order to observe the incidence of these features when solving the granular flow equations. This work identifies the benefits and weaknesses of the SMPM to solve this set of equations, and thus, it is possible to conclude that the SMPM provides an appropriate solution to the granular flow equations proposed by Iverson and Denlinger (2001) and comparable predictions for the experimental data. ABSTRACT A spectral multidomain penalty method solver for the numerical simulation of granular avalanches


Introduction
The study of granular avalanches provides a basis for understanding a variety of mass movements, such as rock avalanches, snow avalanches, debris flows and pyroclastic flows Denlinger and Iverson (2001). Some authors , Savage and Hutter, 1989, Wang et al., 2004, Pastor et al., 2015, Koch et al., 1994 have made efforts to perform laboratory experiments and numerical simulations of granular avalanches. The continuum mechanics framework supplies a helpful theoretical standpoint to represent multiple types of flows by using the laws of conservation of mass and momentum. A reduced form to use the conservation laws applied to granular flows is by means of a shallow water approach. The bidimensional equations that describe the movement of a granular mass taking into account the stresses, caused by either the grains or the fluid.
The granular flow equations are prone to develop discontinuities due to its nonlinear nature. Often, the solution to these equations can be obtained by mean of numerical methods such as finite differences, finite volumes, and finite elements, known as low-order techniques. Results using methods of low order are reported by Gray et al., (1999), Hutter and Koch, (1991), Hutter and Greve, (1993), Hutter et al., (2005), Ouyang et al., (2013), Savage and Hutter, (1989), Tai et al., (2002), Wieland et al., (1999), Denlinger and Iverson, (2001), Denlinger and Iverson, (2004), Koschdon and Schafer, (2003) , Vollmoller, (2004), Wang et al., (2004), where the numerical solution shows a similar structure when compared with experimental data. However, the techniques mentioned above present significant numerical error, such as artificial diffusion and spurious oscillations that affect the stability, accuracy, and conservative properties of the solution as is shown by Fletcher, (1991). In fact, Diamessis and Redekopp, (2006) mentioned that numerical stability resulting from low-order finite-difference schemes is caused by the well-known inherent artificial dissipation, which might generate a significant error in the longtime integration. This issue has been handled with a series of robust methods known as shock-capturing schemes such as ENO, WENO, NOC, TVD, HLLC, among others , Tai et al., 2002, Wang et al., 2004, Toro, 2013. However, it is still necessary to perform significant efforts for developing numerical techniques that allow solving the systems of equations as well as possible to describe these natural phenomena (Castro-Orgaz et al., 2015, Iverson, 2014. In order to minimize the error associated with numerical diffusion generated by using low-order techniques, researchers as Hesthaven and Gottlieb, (1996), Reed and Hill, (1973), Diamessis and Redekopp, (2006), Giraldo and Warburton, (2008) developed element-based high-order discontinuous techniques. Benchmark cases with an exact solution were used to evaluate the accuracy and convergence of spectral methods by solving the inviscid shallow water equations (Escobar-Vargas et al., 2012). Based on the evidence from the high-order precision of these methods O(10 -11 ), the potential of an element-based spectral method for solving a set of equations that describe the movement and deposition of granular avalanches is evaluated, taking into account its lack of numerical diffusion. This work was performed with the aim of test the potential of high-order methods for more applied cases since these techniques have been implemented in idealized cases.
The numerical method implemented is based on a collocation approach of multiple non-overlapping quadrilateral subdomains, which are connected by a penalty term that ensures the stability of the solution by imposing a weak continuity at the subdomain interfaces (see Appendix A). This technique is called spectral multidomain penalty method, henceforth SMPM Gottlieb, 1996, Escobar-Vargas et al., 2012). In order to maintain the stability of the spatial discretization technique, a strong-stability preserving Runge-Kutta was implemented for the temporal derivative (Appendix B), simultaneously, a spectral filter was imposed to control the numerical oscillations (Appendix C). The potentiality of the SMPM is evaluated in two ways. The first one is by means of a grid independence test that shows the error in mass conservation expression that is a function of hydraulic diffusivity D and the initial pore pressure ratio  such as is shown in Iverson and Denlinger, (2001). k a p / is the lateral earth pressure coefficient; and  is the momentum correction factor that comes from the integration of the advective terms of the Navier-Stokes equations (Zhou, 2004). The structure of equations (2-5) is similar to the one of inviscid shallow water equations, where the source terms S x and S y take into account the stresses generated by fluid and solids. Details of these features are shown in equations 6 and 7.
where  int is the internal friction angle of displaced materials,  bed is the friction angle of bed material, n f is the fluid volume fraction (i.e., the porosity),  is dynamic viscosity of the fluid, and  x ,  y are the components of local angle on the bed in the x and y directions, which are described in Section 3.1. The lateral earth pressure coefficient (k a p / ) that appears on Equations (3-7) can be represented by This coefficient associates the horizontal and vertical stresses (Savage and Hutter, 1989;Iverson and Denlinger, 2001). Equation 8 shows that k a p / has a discontinuous nature, indicated by the plus/minus sign, which jumps from active to passive states (i.e., ∂ ∂ > u x / 0 or ∂ ∂ < u x / 0 for the x direction and for the y direction). This feature makes the solution to be prone to become unstable. A better formulation of this coefficient is given by a regularized version of equation 8 [Tai and Gray, 1998]. The discontinuity is avoided by using the following expression: where A is assumed to be a monotonically decreasing functions (Equation 10), which depends of the gradient of the velocity.
In Equation 10,  determines the steepness of the transition between active and passive coefficient, and ∂ ∂ u x / 0 is chosen as the center of the transition (i.e., when ∂ ∂ = u x / 0) (Tai and Gray, 1998). This regularized lateral earth pressure coefficient is applied in both the x and y directions and yields a more stable solution (see section 5.2). Another version of the lateral earth pressure coefficient was implemented in Tai et al., (1999), Gray et al., (1999), Tai et al., (2001), Wieland et al., (1999), Gray, (2003, where the relationship between vertical and horizontal stresses in the y direction is a function of the lateral earth pressure coefficient in the x direction. introduced by the numerical scheme with the spatial refinement. The grid independence test shows that the error in mass conservation does not improve as much as is reported in the literature for the case of a granular avalanche. The second means of evaluating the SMPM involves comparing the numerical simulations and the experimental measurements performed by Denlinger and Iverson, (2001). The comparison shows that the SMPM provides comparable results for depths of the flow concerning the experimental measures and other numerical techniques. Finally, a discussion of boundary conditions, lateral earth pressure coefficients and the momentum correction factor were raised due to its incidence when solving the granular flow equations.
The paper is organized as follows. In section 2, the granular flow equations are presented. In section 3, the numerical setup for the case of study, such as initial and boundary conditions, and treatments are described for zones with or without material. The analysis of a grid independence test is performed in section 4; and in section 5, the comparison between experimental and numerical simulation is performed where boundary conditions, lateral earth pressure coefficients and momentum correction factor are analyzed. A discussion and conclusions are presented in section 6 and 7, respectively. In Appendices A, B, and C we present the numerical techniques with which the equations are discretized and stabilized.

Governing equations of the granular flow
Many granular avalanche models have been proposed based on the 2D shallow water approach as shown by Savage and Hutter, (1989), Iverson and Denlinger, (2001), Denlinger and Iverson, (2004), Hutter et al., (2005), Wang et al., (2004), where the main differences are focused on evaluating some features of the stresses. In this sense, the model proposed by Iverson and Denlinger, (2001) was selected as a sample in order to evaluate the potential of the SMPM for solving granular flow equations (see Appendix A). The temporal derivative is handled by a strong-stability preserving Runge-Kutta (See Appendix B). The model can be presented in conservative form as where q denotes the conservative variables vector. F(q) and G(q) represents the horizontal fluxes in the x and y direction, respectively, and S(q) is the vector that contains the source terms, which dissipate the energy (Escobar-Vargas et al., 2012). By decomposing this system of equations we have In this set of equations, h is the flow depth, m hu x = and m hv y = , where u and v are the averaged velocities in x and y directions. g g x = sin , g y = 0 and g g z = cos  are the gravitational acceleration components in x, y and z directions, respectively, g m s = 9 81 2 .
/ is the gravitational constant and  is the / is the ratio of the basal pore fluid pressure to the total normal basal stress.  is evaluated from an analytical summation In addition, a constant lateral earth pressure coefficient (k nat ) was employed assuming that internal and bed friction angles are equal (   int bed � = ). The main importance of implementing this constant coefficient is to perform the evaluation of other parameters in a more controlled manner by avoiding the discontinuous nature of equation 8 or even the softened function 10, by using equation 11 .
The set of equations 1 -5 belongs to a family of models in which one set of equations is employed, and the stresses are split into two components, solid and fluid. This characteristic makes a significant difference since it avoids the difficulties from a real two-phase model, where a set of equations is needed for each phase and its coupling (Pudasaini, 2012). The fraction of the fluid is a constant parameter, meaning that the proportion of fluid is homogeneous, and it does not vary in time. The role of the fluid fraction in the configuration used here (i.e., granular ow) seems to be irrelevant since the dynamic viscosity of the air is a small value (i.e., 2x10 -5 ). Hence, stresses produced by the fluid are weak regarding the stresses caused by the solids.

Granular flow case
Geometrical basis and initial conditions of the granular flow Denlinger and Iverson (2001) presented a small-scale experiment of a dry granular avalanche. They built a 1.2 m long rectangular flume with an inclined bed at 31.4° with respect to the horizontal surface. There is a curved section with a 10 cm radius of curvature between the inclined and the horizontal surface ( Figure 1). Dry quartz sand with a grain diameter of 0.5 mm was used. The initial setup of the sand is wedge-shaped and is located behind a vertical gate positioned at a 62.5 cm slope down from the top of the flume, which is discharged suddenly when the gate opens ( Figure 1a) . The flume and gate width is 20 cm. The geometry of the bed varies only in the x direction. This variation can be written as Based on Equation 12, the geometrical configuration of the flume is split into three parts: the slope (i.e., 0 ≤ ≤ x x l ), the curved zone (i.e.,), and the horizontal zone (i.e, x x x l r ≤ ≤ ) (Figure 1a). In Equation 12, b(x) describes the elevation of the flume, b x , is the maximum height of the curved zone, at left side of the curved part, x l is any point of x, on the curved zone (r x ), x r is any point of x on r x , and m 0 is the slope of the inclined surface.
It is possible to obtain the slope angle from the geometrical configuration of the basis (i.e., ) and the curvature (i.e., ∂ ∂  x x / ) that appears in equations 6 and 7. The initial condition was built as a wedge (Figure 1), which is described by the following equation: The initial condition (h i ) is split into three points, left (x l ), central (x c ) and right point (x r ). The wedge that represents the initial condition is constructed between the left and right points. m 1 and m 2 are the slope of each side of the triangle. Table 1 shows the properties of the materials used in the numerical simulation presented in this work     .

Boundary conditions
An appropriate implementation of boundary conditions is one of the most important aspects to obtain a robust numerical approximation of any partial differential equation. The relevance of boundary conditions is presented in Section 5.1, where it is shown that different solutions are obtained by using different treatments on boundaries. In order to analyze the influence of the boundary conditions in the solution, two types of boundary conditions were implemented and are presented below. Figure 2 shows the representation of two boundary conditions. Reflective boundary conditions are imposed in both cases, on boundaries perpendicular to the flow direction (i.e., the velocity vector is u n ). In this case, the right boundary must be far enough from the mass movement to avoid any influence in the dynamic of the numerical experiment. Thus the mass can stop by the frictional force solely and not by the effect of an obstacle. Also, slip boundary conditions were imposed on boundaries parallel to the flow direction, u n · = 0 (Figure 2a). In the second case, "non-slip" boundary conditions were implemented on boundaries parallel to the flow direction, which consists of imposing a velocity value equal to zero, u = 0 (Figure2b).

Flow direction
Flow direction

The tracking of the free surface
To solve the granular avalanche equations (Equation 1), it is necessary to define a method to treat regions without material (i.e., h = b). Shallow water equations are defined by a height such that h H o = +, where H 0 is an average depth and  is the displacement of the free surface (Kundu and Cohen, 2008). From solving the system of equations, the velocities u and v are computed by means of conservative variables with Equation 14.
The use of average depth H 0 allows the appearance of undesired oscillations below the average height of the flow that cannot appear in dry zones (h = b), when a rigid bed exists. In this work, the average depth is H 0 = 0 for the granular flow equations, then h = . In areas without material (i.e., h = 0) the solution of the set of equations in the conservative form generates a mathematical indetermination (i.e., division by zero to compute u and v , see Equation 14). In order to avoid this problem, a tolerance value was defined (e.g., Γ = − 10 4 ). If h < Γ, a trivial solution is imposed (i.e., u v = = 0) (Bradford and Sanders, 2002;Pudasaini and Hutter, 2007). Otherwise, if h > Γ, the velocities are computed without imposing any restriction (Pudasaini and Hutter, 2007).

Stability condition and mass conservation
The stability condition for the SMPM is given by the Courant number, which is defined by Giraldo and Warburton (2008) as follow: Conservation of mass was used as the main criterion by which to evaluate the conservative properties of the method. The error introduced by the numerical method was calculated with respect to the initial mass. The mass of the system can be calculated as where  represents the domain. The metric to evaluate mass conservation is based on the relative error, defined with respect to the initial values of mass as In equation 17, R M is the relative error of the mass, and M 0 , M t are the corresponding values for the mass at the initial time (t = 0) and final time (t) of the simulation. Thus, the loss or gain of mass is specified at each time step. The stopping criterion was focused on the experimental data since the last measurement available to compare is at 1.5 seconds. In addition, the variation of the position of the centroid of the mass and its velocity were kept below O(10 -2 ) as a secondary criterion. By another hand, the error in energy conservation was not verified in this work because the source terms in the momentum equations are responsible for the energy dissipation.

Grid independence test
The grid independence test was performed to determine the optimal degree of grid refinement at which the solution does not change when the mesh is reduced. The error in mass conservation (equation 17) was chosen as a reference quantity to assess the differences of the solutions among multiple grids. The mesh refinement was handled with two different strategies. The first option is developed by dividing the domain into multiple subdomains; this is called h-refinement. The second option keeps the numbers of subdomains fixed while increasing the polynomial degree; this is called p-refinement (Boyd, 2001). Figure 3 shows the results of error in mass conservation for the test case with the two types of boundary conditions (slip and non-slip) described in Section 3.2. This was performed in order to observe whether the error is increased with respect to the slip boundary conditions when the non-slip boundary conditions are imposed. The simulation was performed for 9 cases of h-refinement, where the subdomains were kept square (Figure 3), and 19 cases of p-refinement (i.e., 2  N  20). It is possible to observe that the error in mass conservation does not improve when the refinement of the mesh is increased. It is noticeable that the non-monotonicity of the grid independence test disappears when p-refinement is greater than 16th degree, independently of the h-refinement for both boundary conditions. The non-monotonic behavior presented in the grid independence test makes more difficult the choice of an optimal spatial discretization. In addition, the error does not significantly decrease, as is reported with benchmark cases (Escobar-Vargas et al., 2012).
In addition, the implementation of the collapse of a cylinder and a Gaussian bump of sand for two different average depths (H 0 = 1 and H 0 = 0) (Figure 1b and 1c, respectively) were performed. Considering the Gaussian bump (Figure 4a), it was obtained a very low numerical error O(10 -9 ) both with and without average depth, whereas the cylinder generates a greater numerical error (Figure 4b). The tracking of the free surface provides an additional error caused by the zero-average depth (H 0 = 0). Otherwise, an average depth greater than zero (H 0 = 1) helps to give a lower error. If the value of the tolerance for the ʽʽwet/dry'' condition is 10 -4 , the numbers after the fifth decimal place are neglected when the velocities are computed. These velocities are inputs for the next time-step, which affects the accuracy. Although it is natural that hyperbolic equations develop sharp fronts, with these set up cases it was possible to observe that discontinuities might not mature. This is because of the physics of the modeled phenomena included by means of the source terms since these can be stopped soon. Figure 5 presents the comparison between the experimental data (D-I2001)  and the numerical prediction given by the SMPM with a grid that showed the lower error in mass conservation obtained from the grid independence test (Figure 3). Also, slip boundary conditions and the constant lateral earth pressure coefficient were imposed.

Boundary conditions analysis
The numerical prediction is akin to the experiment, particularly in the deposit of the granular flow. The maximum heights of the deposit given by the numerical prediction are approximated to the maximum heights of the experiment with differences of 2 mm, approximately. However, the movement of the material is not well captured, because the "adhered" material on the lateral walls does not appear in the numerical prediction as was observed in the experiment. Figure 6 presents the comparison between the experimental data (D-I2001) and the numerical prediction given by the SMPM with non-slip boundary conditions and the constant lateral earth pressure coefficient. Nonslip boundary conditions allow us to represent the tail of the granular flow while the granular mass is going downhill. The "adhesion" of the material on the lateral walls given by the numerical prediction is similar to the experimental data (the tail of the flow). However, a deposit zone is generated caused by these boundary conditions (the tail of the flow is still present), which is not the case in the experiment. Denlinger  solved this set of equations by means of a method called Harten-Lax-van Leer-Contact (HLLC) (Toro, 2013). A quantitative comparison among the numerical methods with respect to the experimental data is presented for five different stages ( Table 2). The SMPM allows us to get a comparable solution with the HLLC method, since the error in mass conservation has the same order of magnitude for the case of study, with respect to the experimental data.
It is possible to see from the Table 2 that the SMPM has generated less differences with respect to the experiment in the first three stages (i.e., while the material is descending the inclined plane), whereas the HLLC method produced better results in the last two stages (i.e., when the material is forming the deposit, and finally is stopped).

a) b) Experiment (D -I2001)
Numerical prediction (SMPM)  Another comparison between numerical results and experimental data is presented in Figure 7. The continues and discontinuous lines represent two ways to compute the centroid of the mass and the centroid of the velocity field, such as is shown in Pudasaini and Hutter, (2007, p. 149). The asterisk represents the centroid of the experiment computed from the cloud of points as It is possible to observe that the numerical prediction matched in a good way the position of the center of the moving mass computed from the experimental data, however, the centroid of the velocity field shows difficulties that the prediction fits the experiment. V is referring to the total volume.

Analysis of lateral earth pressure coefficient (k a p / )
The lateral earth pressure coefficient (k x y a p , / ( ) ) relates the vertical stresses with the horizontal ones of the granular flow. Figures 5b and 6b show the solution generated with a constant k nat (equation 11), whereas Figure 8b presents the solution given by the slip boundary conditions and a regularized k x y a p , / ( ) (equations 9-10). This specific coefficient tends to produce asymmetric solutions caused by nonphysical oscillations, and hence, steep localized gradients from the velocity can be produced. Since the derivative of the velocity is an input of the lateral earth pressure coefficient, this has repercussions in the behavior of this coefficient and to the next time step. These oscillations can be mitigated either by decreasing the order of the polynomial or by handling the spectral filter (Appendix C). The use of a regularized k x y a p , / ( ) allows us to represent in a better approach the physics of the processes that occur in these types of flows, by means of the smooth connection between the two states of stress (dilation or compression of the material). By another hand, the use of a constant k nat indicates that neither dilation nor compression zones are presented, but this behavior does not occur in nature. For example, experimental data presented by Pudasaini and Hutter, (2007) showed that velocities can change from one point to another. However, a constant k nat was used because it makes it possible to obtain a more stable solution, whereas with a regularized k x y a p , / ( ) the solution is more likely to be unstable, even more with the discontinuous coefficient k x y a p , / ( ) .

Analysis of momentum correction factor ()
The momentum correction factor  (see equations 3 and 4) is a parameter derived from depth-integrating the Navier-Stokes equations, undertaken to obtain the shallow water equations (Zhou, 2004, Trujillo-Vela, 2015. The physical meaning of  is that it modifies the profile of the velocity of the flow. For most of the cases this factor is taken as the unity by assuming that the velocity profile is almost uniform through the depth of the avalanche Iverson, 2001, Savage andHutter, 1989). However, this assumption may not be the most appropriate after of the release, in the proximity of obstacles and near to the deposit where both height and vertical velocity may change considerably and, in some cases, abruptly (Pudasaini and Hutter, 2007). For this reason, the effect of the momentum correction factor is analyzed in this section.  (Figure 9b). If  = 10.5 at t  0.32 s, the flow presents more blunt fronts, and the arrival of material at the zone of deposit is sooner, as is shown in Figure 9d. Figures 10a and 10c show experimental results of the profile of the free surface at the center of the flume at t = 0.53 s and t = 0.93 s, respectively. The best predictions at t = 0.53 s is given by  = 1.15, because this value of  presents a distribution of material similar to the distribution in the experiment. At the same time, if  = 1.0, a significative amount of material has descended the slope and the generation of the deposit begins prematurely, with respect to the experiment. At t = 0.93 s as shown in Figure 10c. The experimental measurement shows that the deposit is almost fully formed (Figure 10c). A good approximation for the tail of the flow and the deposit is reached when  = 1.10 ̶ 1.15 (Figure 10b and d), although the tail is longer than in the experiment. Also, Hutter et al., (2005) show this behavior. On the other hand, there is a numerical oscillation on the left side of the deposit when  = 1.0 ̶ 1.10caused by the slope change in this zone and the selected . Figure 11b shows that the best prediction of the deposit at t = 1.50, which is obtained by  = 1.0 because the tail is less noticeable (i.e., thinner) unlike if other values for the momentum correction factor are used, and the experiment does not show a tail. In addition, the maximum heights of the deposit were predicted better when  = 1.0.
The analysis of the momentum correction factor offers information about changes in the shape of the granular flow prediction by imposing different values of this factor. Therefore, it is possible to obtain a better prediction in the process where the modeler is most interested in, either during the movement or the deposit, by using   1.0 or  = 1.0, respectively. Hence, we propose that future research examine this parameter in order to generate a momentum correction factor that can change depending on what period of movement is being examined when the shallow water approach is used.

Discussion
The results obtained in this work allow us to emphasize some important aspects regarding the solution of granular flow equations. The first aspect is the importance of performing a grid independence test for any numerical method employed in order to obtain an optimal grid. The grid independence test is important because the order of the error can vary significantly for each numerical method and test case. Therefore, this test helps us to find refinement of the grid that presents low error in the mass conservation. Another aspect is the importance of specifying details of the numerical setup, such as initial and boundary conditions, stabilization techniques (Appendix C), and "wet-dry" treatments, since most of the times a part of this information is omitted.
Obviating the error caused by the difficulty of solving hyperbolic systems, the following most important issue that causes the error is the tracking of the free surface for zones without material (i.e., wet-dry cells for shallow water equations), as was shown in Section 3.3. None of the boundary conditions used in this work present a completely satisfactory prediction about the distribution of material near lateral boundaries, since the slip boundary condition solely allow us to represent adequately the deposit of material and the non-slip conditions allow us to predict correctly the movement of the mass. This highlights the importance of exploring different types of boundary conditions, with which a great representation of the movement and deposit of the mass can be obtained.
As shown in Section 2, the use of a discontinuous lateral earth pressure coefficient (k x y a p , / ( ) ), such as was proposed by Savage and Hutter, (1989), makes prone to get unstable solutions. Furthermore, as shown in Section 5.2, the use of a constant lateral earth pressure coefficient promotes the getting of more stable solutions but is not the best choice from the physical standpoint, because the material suffers thinning and thickening during the movement. Therefore, the regularization of k x y a p , / ( ) as proposed by Tai and Gray, (1998), is required. Finally, the tests performed using the momentum correction factor () allow us to conclude that the most suitable value of the momentum correction factor is  = 1.0 after release and close to the deposit. On the other hand, to obtain a better representation of the distribution of material during movement, to mitigate the steep slope of the propagation front and to avoid the premature generation of the deposit it is necessary to use   1.0. This parameter shows the limitation of the depth-integrated models to represent a wave propagation front when the unity of the  is selected, as is commonly performed.

Conclusions
The solution of granular flow equations with a spectral method such as the SMPM, have not been reported yet. Based on the results reported for the first time in this work, the SMPM solves appropriately the granular flow equations proposed by Iverson and Denlinger (2001) and it provides comparable predictions with the experimental data and others' numerical.
We observed from the results that the error in mass conservation increases significantly when a discontinuity is imposed as an initial condition or the sharp fronts are developed during the movement of the mass and, also when the treatment for the tracking of the free surface is used. By another hand, it is possible to use the advantageous property of hp-refinement of the SMPM, instead of h-refinement, as occurs with low-order techniques.
Slip and non-slip boundary conditions modify the mass distribution considerably. Numerical solutions with slip boundary conditions provide better predictions about the deposit, whereas non-slip boundary conditions provide better predictions about the movement of the granular flow.
Instabilities in the solution of the granular avalanche equations introduced by the discontinuous lateral earth pressure coefficient proposed by Savage and Hutter (1989) can be mitigated by the regularization of this coefficient, such as shown in Tai and Gray (1998). Also, it is possible to obtain a stable solution with a constant lateral earth pressure coefficient.
The best prediction of the beginning of the flow and the deposit are obtained when the momentum correction factor is imposed equal to the unity ( = 1.0), whereas, the best prediction of the movement of the flow is obtained if the momentum correction factor is greater than the unity ( > 1.0).

Appendix A. Spectral multidomain penalty method
In order to assess the potential of high-order discontinuous methods for real applications and in order to reduce the numerical error introduced by low-order techniques in this paper, we performed the solution of the granular flow equations (equation 1) with a method that combines an exponential convergence, a weak artificial dissipation, and an adaptability based in classic finite elements-volumes techniques (Boyd, 2001). The spectral multidomain penalty method (SMPM) is based on a collocation approach of multiple nonoverlapping quadrilateral subdomains, which are connected by a penalty term that ensures the stability of the solution by imposing a weak continuity at the subdomain interfaces Gottlieb, 1996, Escobar-Vargas et al., 2012). The formulation of SMPM will be presented in two parts: the subdomain interior and the interfacial treatments.
Subdomain interior. In the interior of each subdomain any function q(x, y, t) can be approximated by using N-th order Lagrange interpolating polynomials on a Gauss-Lobatto-Legendre grid (GLL) (Hesthaven and Gottlieb, 1996).
where q(x i , y j , t) is the value of function at the discrete point (x i , y j ), and l i (x), l j (y) are the i-th and j-th Lagrange interpolating polynomials based on the GLL nodes in the and y directions, respectively. The spatial derivatives in the x-direction in the global coordinate system are approximated as where ∂ ∂  / x represents the mapping from the local coordinates system  ∈ −     1 1 , , given by GLL points into the global coordinates system x ∈. D ik is the Legendre spectral differentiation matrix, which is computed following the procedure suggested by Costa and Don (2000).
Interfacial treatment. The penalized form of the granular flow equations in a collocation point situated on a subdomain boundary requires an additional penalization term, which allows us establishes the continuity between adjacent points (Escobar-Vargas, 2012).
where  are weights associated to the GLL points. For a uniform polynomial approximation, N, in each subdomain can be used a unique value of  = + ( ) 2 1 / N N and ∆ ∆ ∆ s x y = ( ) , which depends on the orientation of the interfaces in the subdomains.