Keywords: Granular avalanches; Numerical simulation; Spectral multidomain penalty method; Parameter analysis. Palabras clave: Avalancha granular; simulación numérica; método penalizado de multidominios espectrales; análisis de parámetros. 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), 317- 329. 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 Método penalizado de multidominios espectrales para la simulación numérica de avalanchas granulares ISSN 1794-6190 e-ISSN 2339-3459 https://doi.org/10.15446/esrj.v23n4.77683 Este trabajo presenta la simulación numérica de una avalancha granular a escala experimental basada en elementos de alto orden, con el fin de evaluar el potencial de estas técnicas espectrales para resolver las leyes de conservación en geofísica. La discretización espacial de estas ecuaciones se desarrolló a través del métodos penalizado de multidominios espectrales (SMPM). Los términos temporales se discretizaron utilizando un método Runge-Kutta de preservación de estabilidad fuerte. La estabilidad del esquema numérico se garantiza con el uso de un filtro espectral y un coeficiente de presión lateral de tierras constante o regularizado. El caso de prueba es una avalancha granular que se genera en un canal rectangular de pequeña escala con gradiente topográfico. Se realizó una prueba de independencia de malla para aclarar el orden del error en la conservación de la masa producida por los tratamientos previamente mencionados. Las predicciones numéricas se comparan con las mediciones experimentales realizadas por Denlinger & Iverson (2001). Además, se analizaron las condiciones de contorno, y los parámetros como el coeficiente de presión lateral de tierras y el factor de corrección de momentum para observar su incidencia al resolver las ecuaciones del flujo granular. Este trabajo identifica las bondades y debilidades del SMPM para resolver este conjunto de ecuaciones y es posible concluir que el SMPM proporciona una solución adecuada de las ecuaciones propuestas por Iverson & Denlinger (2001). Además, las predicciones son comparables a los datos experimentales y resultados numéricos dados por otros esquemas. RESUMEN Record Manuscript received: 05/02/2019 Accepted for publication: 20/09/2019 EARTH SCIENCES RESEARCH JOURNAL Earth Sci. Res. J. Vol. 23, No. 4 (December, 2019): 317-329 Mario Germán Trujillo-Vela1, Jorge Alberto Escobar-Vargas1, Alfonso Mariano Ramos-Cañón2 1Departamento de Ingeniería Civil, Pontificia Universidad Javeriana, Bogotá, Colombia. 2Departamento de Ingeniería Civil y Agrícola, Universidad Nacional de Colombia, Bogotá, Colombia. * Corresponding author: mario.trujillo@javeriana.edu.co G EO PH YS IC S 318 Mario Germán Trujillo-Vela, Jorge Alberto Escobar-Vargas, Alfonso Mariano Ramos-Cañón Nomenclature h (m) Depth of the flow. u (m/s) Vector velocity. u v, (m/s) Average velocities. m mx y,� Conservative variables. q Vector of conservative variables. F q G q( ) ( ), Fluxes of the conservative variables. S q( ) Vector of source terms. S Sx y, Source terms. Ho (m) Average depth.  (m) Tolerance value  (m) Displacement of the free surface.  Ratio of the basal pore fluid pressure. Pbed (kg/m s3) Total basal normal stress. k x y a p, /( ) Lateral earth pressure coefficient for x and y directions. ka p/ Lateral earth pressure coefficient. knat Constant lateral earth pressure coefficient.  Steepness of the transition between active and passive coefficients. gx, gy, gz (m/s2) Gravitational components.  Momentum correction factor. b (x) (m) Topographical basis in function of x. hi (x) (m) Depth of material as initial conditions. m0 , m1, m2, (m) Steepness of the topographical and initial conditions. xl, xc, xr , (m) Left, central and right points in the x axis. rx (m) Radius of curvature. x , x (°) Components of local angle. n f Fraction of the fluid.  (Pa-s) Dynamic viscosity of the fluid.  (kg/m3) Density of the mixture. f (kg/m3) Density of the fluid. s (kg/m3) Density of the soil. bed (°) Basal friction angle. int (°) Internal friction angle. CFL Courant number. t Time-step. U Magnitude of the velocity in a cell. V (m3) Total volume of the moving mass. S Grid spacing. M Mass of the system. M0, Mt (m 3) Corresponding values for the mass at the initial and final times of the simulation.  Computational domain. RM Relative error of the mass. Xc exp( ) Position of the centroid of the experimental moving mass. Uc exp( ) Velocity in the x-direction of the experimental centroid of the moving mass. 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 (Denlinger and Iverson, 2001, Savage and Hutter, 1989, Gray et al., 1999, 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 long- time 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 (Denlinger and Iverson, 2001, 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 (Hesthaven and 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 319A spectral multidomain penalty method solver for the numerical simulation of granular avalanches expression that is a function of hydraulic diffusivity D and the initial pore pressure ratio  such as is shown in Iverson and Denlinger, (2001). ka 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 Sx and Sy take into account the stresses generated by fluid and solids. Details of these features are shown in equations 6 and 7. S g h sgn u g u x h n u h n h x x z x bed f f = − ( ) −( ) + ∂ ∂     − + ∂ 1 32 2 λ θ φ µ ρ µ ρ tan uu x∂ 2 − ∂ ∂     ∂ ∂ −( )  + ∂ ∂ sgn u y hk y g h n h u y a p z int f / sin1 2 2λ φ µ ρ (6) S g h sgn v g v y h n v h n h y y z y bed f f = − ( ) −( ) + ∂ ∂     − + ∂ 1 32 2 λ θ φ µ ρ µ ρ tan vv y∂ 2 − ∂ ∂     ∂ ∂ −( )  + ∂ ∂ sgn v x hk x g h n h v x a p z int f / sin1 2 2λ φ µ ρ (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 (ka p/ ) that appears on Equations (3-7) can be represented by k cos a p int bed int / / cos tan �= ± − +( )    −2 1 1 1 1 2 2 1 2 2    (8) This coefficient associates the horizontal and vertical stresses (Savage and Hutter, 1989; Iverson and Denlinger, 2001). Equation 8 shows that ka p/ has a discontinuous nature, indicated by the plus/minus sign, which jumps from active to passive states (i.e., ∂ ∂ >u x/ 0 or ∂ ∂ v y/ 0 or ∂ ∂ Γ, 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: CFL t U g S = +( )        max ∆ ∆ h (15) where t is time step size, U u v= +2 2 is the magnitude of average velocities into the cell, gh is average geopotential height into each cell and ∆ ∆ ∆S x y= +2 2 is the grid spacing. The sum of U gh+ is known as the characteristic velocity at which the wave travels (LeVeque, 2002). For SMPM, the maximum number of Courant that ensures the stability is CFL 1 (Escobar- Vargas et al., 2012). 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 M hd= ∫ Ω Ω � (16) 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 R M M MM t= −( )0 0 (17) In equation 17, RM is the relative error of the mass, and M0, Mt 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 322 Mario Germán Trujillo-Vela, Jorge Alberto Escobar-Vargas, Alfonso Mariano Ramos-Cañón 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 (H0 = 1 and H0 = 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 (H0 = 0). Otherwise, an average depth greater than zero (H0 = 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. Comparison with experimental data Boundary conditions analysis Figure 5 presents the comparison between the experimental data (D- I2001) (Denlinger and Iverson, 2001) 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. 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. Non- slip 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. a) b) Figure 3. Mass conservation error with a) slip and b) non-slip boundary conditions. (Subx and Suby represents the subdomain numbers in x and y directions, respectively). b)a) Figure 4. Mass conservation error in function of time with and without average depth and without a topographic gradient, a) Gaussian bump and b) cylinder. 323A spectral multidomain penalty method solver for the numerical simulation of granular avalanches Denlinger and Iverson, (2001) 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) Figure 5. a) Granular flow experiment. b) Numerical prediction given by SMPM for spatial refinement with low numerical error (N = 8, Subx = 108 and Suby = 18), slip boundary conditions, and constant lateral earth pressure coefficient. a) b)Experiment (D -I2001) Numerical prediction (SMPM) Figure 6. a) Granular flow experiment. b) Numerical prediction given by SMPM for spatial refinement with low numerical error (N = 8, Subx = 108 and Suby = 18), non-slip boundary conditions, and constant lateral earth pressure coefficient. 324 Mario Germán Trujillo-Vela, Jorge Alberto Escobar-Vargas, Alfonso Mariano Ramos-Cañón Table 2. Euclidean norm of the error for slip and nonslip SMPM, and for HLLC, regarding the experimental data. Time (s) SMPMslip SMPMnonslip HLLC 0.01 0.2227 0.2197 0.4428 0.32 0.1634 0.1583 0.1796 0.53 0.1694 0.1733 0.2108 0.93 0.5699 0.6237 0.4904 1.50 0.4761 0.5739 0.4157 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 X X c exp c exp c exp ( ) ( ) ( )= ∑ ∑ = h x h U t i i i , ∆ ∆ (18) 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. Figure 7. Comparison between numerical results and experimental data of the centroid of the moving mass and the velocity field. a) The position of the centroid of the descending mass at each time step. b) Centroid of the velocity field. V is referring to the total volume. Analysis of lateral earth pressure coefficient (ka 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 knat (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 knat 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 knat 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 (Denlinger and Iverson, 2001, Savage and Hutter, 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. Figures 9 to 11 show a longitudinal section at the center of the flume for the experiment (in each Figures a and c) and the numerical results (in each Figures b and d) at five different times (0.10, 0.32, 0.53, 0.93, and 1.50 seconds). Heights are shown on the ordinate and the longitudinal length is shown on the abscissa. The numerical results were computed for five values of  (1.0, 1.05, 1.10, 1.15, 1.20). Figures 9a and 9c show the profile of the free surface at the center of the flume of the experimental results at t = 0.10 s and t = 0.32 s, respectively. Figures 9b and 9d show numerical predictions at the same times. The best prediction of the front at t = 0.10 s is given by  = 1.0, since the numerical results represent the front as steep as in the experiment. Results with   1.0 show very blunt fronts (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 325A spectral multidomain penalty method solver for the numerical simulation of granular avalanches a) b)Experiment (D-I2001) Numerical prediction (SMPM) Fig. 8. a. Granular flow experiment. b. Numerical prediction given by SMPM when the spatial refinement presents low numerical error (N = 8, Subx = 108 y Suby = 18), slip boundary conditions and regularized lateral earth pressure coefficient (k x y a p, /( ) ). a) b) c) d) Experiment (D-I2001), = 0.10 s Experiment (D-I2001), = 0.32 s Numerical prediction SMPM, = 0.32 sNumerical prediction SMPM, = 0.10 s Figure 9. x versus h plots for different  at times t =0.10 s and t =0.32 s, and the comparison with the experimental measurements for the same time. 326 Mario Germán Trujillo-Vela, Jorge Alberto Escobar-Vargas, Alfonso Mariano Ramos-Cañón 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. a) b) c) d) Experiment (D-I2001), = 0.53 s Experiment (D-I2001), = 0.93 s Numerical prediction SMPM, = 0.93 sNumerical prediction SMPM, = 0.53 s Figure 10. x versus h plots for different  at times t = 0.53 s y t = 0.93 s, and the comparison with the experimental measurements for the same time. a) b) Experiment ( D-I2001 ), = 1.50 s Numerical prediction SMPM, = 0.50 s Figure 11. x versus h plots for different  at time t = 1.50 s, and the comparison with the experimental measurements for the same time. 327A spectral multidomain penalty method solver for the numerical simulation of granular avalanches 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). Acknowledgements The authors acknowledge Richard Iverson and Roger Denlinger for providing valuable experimental and numerical data for comparison with our numerical results. Fundación CeiBA supported this work within the scholarship of doctoral studies of the first author. References Blackburn, H. M., & Schmidt, S. (2003). Spectral element filtering techniques for large eddy simulation with dynamic estimation. Journal of Computational Physics, 186(2), 610-629. Boyd, J. P. (2001). Chebyshev and Fourier spectral methods. Courier Dover Publications, Mineola, New York, 2nd edition, 670 pp. Bradford, S. F. & Sanders, B. F. (2002). Finite-volume model for shallow-water flooding of arbitrary topography. Journal of Hydraulic Engineering, 128(3), 289-298. Canuto, C., Hussaini, M. Y., Quarteroni, A., & Thomas Jr, A. (2012). Spectral methods in fluid dynamics. Springer Science & Business Media. Berlin, Germany, 566 pp. Castro-Orgaz, O., Hutter, K., Giraldez, J. V., & Hager, W. H. (2015). Nonhydrostatic granular flow over 3-D terrain: New boussinesq-type gravity waves? Journal of Geophysical Research: Earth Surface, 120(1), 1-28. Costa, B. & Don, W. S. (2000). On the computation of high order pseudospectral derivatives. Applied Numerical Mathematics, 33(1), 151-159. Denlinger, R. P. & Iverson, R. M. (2001). Flow of variably fluidized granular masses across three-dimensional terrain: Numerical predictions and experimental tests. Journal of Geophysical Research, 106(B1), 553- 566. Denlinger, R. P. & Iverson, R. M. (2004). Granular avalanches across irregular three-dimensional terrain: 1. theory and computation. Journal of Geophysical Research, 109 (F1). Diamessis, P. J., Lin, Y.-C., & Domaradzki, J. A. (2008). Effective numerical viscosity in spectral multidomain penalty method-based simulations of localized turbulence. Journal of Computational Physics, 227(17), 8145-8164. Diamessis, P. J. & Redekopp, L. G. (2006). Numerical investigation of solitary internal wave-induced global instability in shallow water benthic boundary layers. Journal of physical oceanography, 36(5), 784-812. Escobar-Vargas, J., Diamessis, P., & Giraldo, F. (2012). High-order discontinuous element-based schemes for the inviscid shallow water equations: Spectral multidomain penalty and discontinuous Galerkin methods. Applied Mathematics and Computation, 218(9), 4825-4848. Fletcher, C. (1991). Computational techniques for fluid dynamics, volume 1 of Springer series in computational physics. Springer-Verlag, Berlin; New York, 2nd edition, 409 pp. Giraldo, F. & Warburton, T. (2008). A high-order triangular discontinuous Galerkin oceanic shallow water model. International Journal for Numerical Methods in Fluids, 56(7), 899-925. Gottlieb, D., & Hesthaven, J. S. (2001). Spectral methods for hyperbolic problems. Journal of Computational and Applied Mathematics, 128(1- 2), 83-131. Gottlieb, D., & Shu, C. W. (1997). On the Gibbs phenomenon and its resolution. SIAM review, 39(4), 644-668. Gottlieb, D., & Tadmor, E. (1985). Recovering pointwise values of discontinuous data within spectral accuracy. In Progress and supercomputing in computational fluid dynamics, 357-375. Gray, J., Wieland, M., & Hutter, K. (1999). Gravity-driven free surface flow of granular avalanches over complex basal topography. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 455(1985), 93-138. Gray, J. M. N. T. (2003). Rapid granular avalanches. In Dynamic Response of Granular and Porous Materials under Large and Catastrophic Deformations, 3-42. Springer. Hesthaven, J. S. & Gottlieb, A. D. (1996). A stable penalty method for the compressible navier-stokes equations: I. open boundary conditions. Journal on Scientific Computing, 17(3), 579–612. Hesthaven, J. S., Gottlieb, S., & Gottlieb, D. (2007). Spectral methods for time- dependent problems (Vol. 21). Cambridge University Press, 272 pp. Hesthaven, J. S. & Kirby, R. (2008). Filtering in Legendre spectral methods. Mathematics of Computation, 77(263), 1425-1452. Hesthaven, J. S. & Warburton, T. (2007). Nodal discontinuous Galerkin methods: algorithms, analysis, and applications. Springer Science & Business Media. 495 pp. Hutter, K. & Greve, R. (1993). Two-dimensional similarity solutions for finite- mass granular avalanches with coulomb-and viscous-type frictional resistance. Journal of Glaciology, 39(132), 357-372. Hutter, K. & Koch, T. (1991). Motion of a granular avalanche in an exponentially curved chute: experiments and theoretical predictions. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 334(1633), 93-138. Hutter, K.,Wang, Y., & Pudasaini, S. P. (2005). The Savage–Hutter avalanche model: how far can it be pushed? Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 363(1832), 1507–1528. Iverson, R. M. (2014). Debris flows: behaviour and hazard assessment. Geology Today, 30(1), 15-20. Iverson, R. M., & Denlinger, R. P. (2001). Flow of variably fluidized granular masses across three-dimensional terrain: 1. Coulomb mixture theory. Journal of Geophysical Research: Solid Earth (1978–2012), 106(B1), 537-552. Koch, T., Greve, R., & Hutter, K. (1994). Unconfined flow of granular avalanches along a partly curved surface. II. Experiments and numerical computations. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 445(1924), 415-435. Koschdon, K. & Schafer, M. (2003). A Lagrangian-Eulerian finite-volume method for simulating free surface flows of granular avalanches. In Dynamic Response of Granular and Porous Materials under Large and Catastrophic Deformations, 83-108. Springer. Kundu, P. & Cohen, L. (2008). Fluid Mechanics. Elsevier Acad. Press, London, England, 729 pp. 328 Mario Germán Trujillo-Vela, Jorge Alberto Escobar-Vargas, Alfonso Mariano Ramos-Cañón LeVeque, R. J. (2002). Finite volume methods for hyperbolic problems, volume 31. Cambridge university press, Cambridge, England, 529 pp. Maday, Y., Kaber, S. M. O., & Tadmor, E. (1993). Legendre pseudo spectral viscosity method for nonlinear conservation laws. SIAM Journal on Numerical Analysis, 30(2), 321-342. Navarra, A., Stern, W., & Miyakoda, K. (1994). Reduction of the Gibbs oscillation in spectral model simulations. Journal of climate, 7(8), 1169-1183. Ouyang, C., He, S., Xu, Q., Luo, Y., & Zhang, W. (2013). A MacCormack-TVD finite difference method to simulate the mass ow in mountainous terrain with variable computational domain. Computers & Geosciences, 52(0), 1-10. Pastor, M., Blanc, T., Haddad, B., Drempetic, V., Morles, M. S., Dutto, P., Stickle, M. M., Mira, P., & Merodo, J. F. (2015). Depth averaged models for fast landslide propagation: mathematical, rheological and numerical aspects. Archives of Computational Methods in Engineering, 22(1), 67-104. Pudasaini, S. P. (2012). A general two-phase debris ow model. Journal of Geophysical Research: Earth Surface, 117(F3), 799-819. Pudasaini, S. P. & Hutter, K. (2007). Avalanche dynamics: dynamics of rapid flows of dense granular avalanches. Springer, Berlin, 571 pp. Reed, W. H. & Hill, T. (1973). Triangular mesh methods for the neutron transport equation. Technical report, Los Alamos Scientific Lab., N. Mex. (USA). Savage, S. B. & Hutter, K. (1989). The motion of a finite mass of granular material down a rough incline. Journal of Fluid Mechanics, 199, 177- 215. Spiteri, R. J. & Ruuth, S. J. (2002). A new class of optimal high order strong stability-preserving time discretization methods. SIAM Journal on Numerical Analysis, 40(0), 469-491. Tai, Y. & Gray, J. (1998). Limiting stress states in granular avalanches. Annals of Glaciology, (26), 272-276. Tai, Y.-C., Hutter, K., & Gray, J. (2001). Dense granular avalanches: mathematical description and experimental validation. In Geomorphological fluid mechanics, 339-366. Springer. Tai, Y.-C., Noelle, S., Gray, J., & Hutter, K. (2002). Shock-capturing and front- tracking methods for granular avalanches. Journal of Computational Physics, 175(1), 269-301. Tai, Y.-C., Wang, Y., Gray, J., & Hutter, K. (1999). Methods of similitude in granular avalanche flows. In Advances in Cold-Region Thermal Engineering and Sciences, 415-428. Toro, E. F. (2013). Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Science & Business Media, Dordrecht; New York, 3rd edition. Trujillo-Vela, M. G. (2015). Simulación numérica de flujos de material desagregado. Método penalizado de multidominios espectrales. M.Sc. Thesis, Maestría en Hidrosistemas, Pontificia Universidad Javeriana, Bogotá D. C., Colombia Vandeven, H. (1991). Family of spectral filters for discontinuous problems. Journal of Scientific Computing, 6(2), 159-192. Vollmoller, P. (2004). A shock-capturing wave-propagation method for dry and saturated granular flows. Journal of Computational Physics, 199(1), 150-174. Wang, Y., Hutter, K., & Pudasaini, S. P. (2004). The Savage-Hutter theory: A system of partial differential equations for avalanche flows of snow, debris, and mud. Zeitschrift fur angewandte Mathematik und Mechanik ZAMM, 84(8), 507-527. Wieland, M., Gray, J., & Hutter, K. (1999). Channelized free-surface flow of cohesionless granular avalanches in a chute with shallow lateral curvature. Journal of Fluid Mechanics, 392, 73-100. Zhou, J. G. (2004). Lattice Boltzmann methods for shallow water flows. Springer, New York. 329A spectral multidomain penalty method solver for the numerical simulation of granular avalanches 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 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 (Hesthaven and 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). q qx y t x y t l x l y i N j N i j i j, , , ,( ) = ( ) ( ) ( ) = = ∑∑ 0 0 (19) where q(xi, yj, t) is the value of function at the discrete point (xi, yj), and li(x), lj(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 ∂ ( ) ∂ = ∂ ( ) ∂ ∂ ∂ = ∂ ∂ ( ) = ∑q q q x y t x x y t x x D x y ti i i i i N ik k j , , , , , ,    0 (20) where ∂ ∂ / x represents the mapping from the local coordinates system  ∈ − 1 1, , given by GLL points into the global coordinates system x ∈. Dik 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). ∂ ∂ + ∂ ( ) ∂ + ∂ ( ) ∂ = ( ) + −  ) = ( )∑q F G q S q e e e e l e e e l e l t q x y Q n 1 4 τ ⋅ , J J (21) where n e l,( ) is the unit vector at the subdomain boundary, l represents each one of the four edge of the domain, J represents the information flux between adjacent points, and e is the penalization coefficient (Escobar- Vargas, 2012). τ τ ω e s ≡ = 1 ∆ � (22) 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. Appendix B. Strong stability preserving Runge-Kutta In order to maintain the precision offered by the spatial discretization method, it is used a strong stability preserving Runge-Kutta method (SSP- RK) for the temporal discretization. By considering the following initial value problem d dt q R q= ( ) � (23) In this scheme, the prediction at the time t + 1 (i.e., q t+1) is based on both the existing solution at the time t (i.e., qt) and the evaluation of the forcing terms R(q) at different stages (Spiteri and Ruuth, 2002). q q0 = t (24) q q R qi k i i k k i k kt i= + ( ) = … = − ∑ 0 1 1 2α β, , , , , , ,∆ s (25) q qt s+ =1 (26) For consistency, we must have that k i i k = − ∑ = 0 1 1 , , i= 1, 2,…, s, where s is the number of stages of the SSP-RK approach, i k, and i k, are constant coefficients where all i k, ≥ 0 and i k, = 0 only if i k, = 0. Spiteri and Ruuth (2002) give these coefficients. t is the size of the time-step. Appendix C. Stabilization technique There are multiples methods for improving accuracy and stability of spectral and pseudo spectral techniques, such as physical space filters (Gottlieb and Tadmor, 1985, Gottlieb and Shu, 1997), spectral vanishing viscosity methods (Maday et al., 1993), dealiasing (Canuto et al., 2012) and over integration (Hesthaven and Kirby, 2008). A family of regularization techniques related to the first two above-mentioned is the spectral filters. Since there is not a dissipative term in the SMPM to remove the artificial oscillations, a spectral filter used by Vandeven, (1991), Gottlieb and Hesthaven, (2001), Diamessis et al., (2008) is necessary to reduce the high-frequency modes that appear close to the discontinuities (Hesthaven and Warburton, 2007). It was selected this spectral filter (exponential) because of the stability and robustness that it adds to the method (Hesthaven and Kirby, 2008, Hesthaven et al., 2007). A comparison of different spectral filtering methods and discussion of the advantages and weaknesses of using them is performed in Hesthaven and Kirby, (2008), Blackburn and Schmidt, (2003), Navarra et al., (1994). Based on the results of this work, it can be concluded that the spectral filter allows us to obtain stable solutions, whereas few solutions are obtained when a filter is not used. Thus, it is possible to avoid artificial and persistent oscillations around points where discontinuities appear. Also, the error in mass conservation has a smoother behavior over time when using a filter.