Multilayer shallow model for dry granular flows with a weakly non-hydrostatic pressure

Multilayer shallow model for dry granular flows with a weakly non-hydrostatic pressure C. Escalante ∗ , J. Garres-D´ıaz † , E.D. Fern´andez-Nieto ‡ , A. Mangeney § Abstract The multilayer model proposed in this paper is a generalization of the multilayer non-hydrostatic model for shallow granular flows Fern´andez-Nieto et al. [40], the multilayer model with µ(I) rheology Fern´andez-Nieto et al. [36], and the monolayer model with weakly non-hydrostatic pressure for dry granular flows Garres-D´ıaz et al.[42]. We show that the proposed model verifies a dissipative energy balance. A well-balanced numerical scheme is proposed to solve the equations based on a projection method and a hydrostatic reconstruction for the Coulomb friction terms. In order to reduce the computational cost associated with solving the linear system of the projection method, a precomputing of the initial guess for an iterative solver is proposed. This strategy allows us to reduce the computational time by around 70% when 20 layers are considered. In the numerical tests, we show that the proposed model can recover the in-depth velocity profiles typically observed in lab experiments and capture the flow/no-flow interface that appears in granular avalanches. During the initial stage of granular collapse simulations, the model is shown to improve the approximation of the mass profiles compared to other models and to predict the parabolic shape of the front velocity evolution with time, as observed in lab experiments. Interestingly, our numerical tests show that the ability of the granular flow to overcome obstacles strongly depends on the model used, which is of strong interest for landslide hazard assessment. Keywords: Multilayer models Non-hydrostatic pressure Finite Volume Granular flows. 1 Introduction Many efforts have been devoted in recent years to the understanding of granular flows, due to their importance in a wide variety of natural phenomena. Both dry and fluidized granular flows have been broadly studied from a theoretical point of view, and also many resources have been employed to perform experimental studies, in particular at the laboratory scale (see e.g. [43,28,48,1] for reviews). The understanding and quantification of the complex behavior of these flows is still a challenge, since their dynamics is strongly conditioned by many different physical processes (friction, dilatancy, fluid-grain interactions, particle segregation, etc.), see for instance [65,12,78,6]. In particular, the flow/no-flow local transition that controls erosion/deposition processes, involves many theoretical and numerical difficulties [55,52,12,63,36]. Here, we deal with dry and dense granular flows, which are nowadays described by the µ(I)-rheology, introduced in [53,54]. This rheology is based on a Drucker-Prager plasticity criterion, that defines the deviatoric tensor as      τ=µ(I)p ∥D∥Dif ∥D∥ = 0, ∥τ∥ ≤ µspif ∥D∥= 0, (1) where pis the pressure, µsis the tangent of the repose angle of the material, Dthe strain-rate tensor and µ(I) is a strain-rate and pressure dependent friction coefficient, which will be defined later. It is known that the incompressible µ(I)-rheology is ill-posed for low and high values of the inertial number I(see [7]). Many works have been devoted to ∗Dpto. Matem´atica Aplicada. U. M´alaga, 14071- M´alaga, Spain (escalan[email protected]) †Dpto. Matem´aticas. Edificio Einstein - U. C´ordoba. Campus de Rabanales, 14014-C´ordoba, Spain ([email protected] ‡Dpto. Matem´atica Aplicada I. ETS Arquitectura - U. Sevilla. Avda. Reina Mercedes S/N, 41012-Sevilla, Spain ([email protected]) §Institut de Physique du Globe de Paris, Seismology team, U. Paris-Diderot, Sorbonne Paris Cit´e, 75238, Paris, France. ANGE team, CEREMA, INRIA, Lab. J. Louis Lions, 75252, Paris, France. Institut Universitaire de France (IUF), France, ([email protected]) 1 improve the well-posed character of the µ(I)-rheology. A regularized incompressible µ(I) law, which improves the range of values for Iwhere it is well-posed, was proposed in [5]. This regularized expression was considered in [44,6] among others. A compressible version of the µ(I), called µ(I), ϕ(I)-rheology, has also been investigated (see e.g. [17,77]), considering that the solid volume fraction ϕis no more constant. Some authors suggested that this compressible version improves the well/ill-posed character of the µ(I) for steady problems (see [8,51]), although this is not necessarily true for transient problems (see [80]). As commented above, when the compressible rheological law is considered, the solid volume fraction ϕis also variable (depending on Iin the aerial case or on the viscous number Jin the submerged case), and dilation/contraction of the material occurs (see [77]). Dilatancy effects in a shallow depth-averaged model for dry avalanches was recently investigated in [11], where the model is obtained by removing the fluid phase in the two-phase model for debris flows with dilatancy [12]. In addition, in the last years, some authors have introduced non-local effects in the constitutive law (see e.g. [76,16]). However, the local µ(I)-rheology continues being very popular, and actually the most accepted law, since many questions about non-local models still need to be clarified. In the following, we focus on the incompressible µ(I)-rheology. This rheology has been used in models solving the full Navier-Stokes equations for free surface flows. The yield stress can be dealt with by defining a regularized µ(I)-viscosity as in the 2D solver of [55] or [56]. A different approach was followed in [52,63], where authors used a finite element solver combined with an Augmented Lagrangian method, also modelling sidewall friction effects. They showed a qualitative good agreement with granular collapse laboratory experiments. Solving the full Navier-Stokes system however requires huge computational efforts. For this reason, depth-averaged shallow models appear as an interesting alternative. Many authors have proposed such models for granular flows since the pioneering Savage-Hutter model [79] where Mohr-Coulomb friction terms with constant friction coefficients are added to the classical Saint-Venant system. A more sophisticated friction coefficient was proposed in [74,73,75], which depends on the thickness and velocity of the flow, e.g. on its Froude number. This variable friction coefficient has been widely used to model dry granular flows (see [62,60,72,19] among many others). Interestingly, when these depth-averaged models are considered in tilted (or local) coordinates over a reference plane or bottom, the flow could start flowing for slopes that could differ (being lower or higher) from the angle of repose of the material. This was shown in [10], where authors proposed a correction of the friction coefficient in terms of a second order correction of the pressure. This correction allows the model to verify the physical threshold to flow, defined by the angle of repose. A rigorous derivation of a shallow model with the µ(I)-rheology under a dimensional analysis is presented in [46], where authors also included a second order viscous term in the form of a horizontal diffusion term. To derive the model, they needed to assume a Bagnold profile for the downslope velocity even though velocity profiles are known to vary significantly from linear, Bagnold or S-shape profiles during the different phases of the flow (initial acceleration, developed flow, deceleration) (see [43]). The main drawback of such depth-averaged models is the fact that all the in-depth information is lost due to the shallow hypothesis and the depth-integration procedure. In particular, it makes it impossible to reproduce the lower static layer that exists in granular collapse or during the stopping phase of granular flows. This flow/no-flow (also called flowing/static) interface is shown to depend on the derivatives of the downslope velocity with respect to the direction normal to the topography, an information intrinsically missing in depth-averaged models [57]. As an intermediate step between full Navier-Stokes and depth-averaged models, multilayer (or layer-averaged) models [3,39] have been shown to be able to recover important in-depth variations of the flow properties such as velocity profiles. In [39] a multilayer model is deduced as the set of equations verified by a particular weak solution of the full Navier-Stokes system. There, the domain is subdivided in shallow layers and the velocity field is approximated by a layerwise constant function. Such a model is shown to reproduce the different shapes of the velocity profiles associated with different flow regimes. For example, in [36], simulations performed with a multilayer model for dry granular flows with the µ(I)-rheology were successfully compared to granular collapse experiments in [61], showing an excellent qualitative agreement concerning the shape of the deposits, of the velocity profiles, or the influence of an erodible bed on the travel distance (runout) of these flows. Moreover, this multilayer model was able to approximate the flow/no-flow interface separating the mobile and static layers of material during granular column collapse. The multilayer approach makes it also possible to properly account for sidewalls friction that has a strong influence on the flow dynamics [37]. All these depth-averaged models assume that the pressure is hydrostatic. However, during the initial destabilization, the mass geometry may significantly differ from a shallow layer, thus inducing strong non-hydrostatic effects arising from the non-negligible acceleration in the direction normal to the topography. This could explain part of the discrepancies observed when comparing the results of shallow depth-averaged models with granular collapse experiments, in particular during the initial instants e.g.[36]. Going beyond the hydrostatic assumption, dispersive models have become a very active topic of research in recent years, looking for an improvement of nonlinear dispersive properties. Two families of dispersive models can be considered: Boussinesq type systems, where high order derivatives of the variables are included in the model (see, e.g., [15,67] among others); and non-hydrostatic systems (see, e.g., [23,82]), where the pressure is split into hydrostatic and non-hydrostatic counterparts, involving new unknowns and restrictions in the model. Actually, most classical dispersive models for water flows can be reformulated as non-hydrostatic systems, as 2 shown in [32]. The advantage of the second approach is that only first-order derivatives appear in the model, which are easier to treat numerically. In addition, the specific structure of these models and their similarities to the Shallow Water Equations (SWE) enables the extension of numerous familiar numerical schemes for SWE to non-hydrostatic models, as demonstrated in [30]. One of the drawbacks of non-hydrostatic models is that they do not fall into the class of hyperbolic PDE systems. Furthermore, there is a need to address implicit problems to include the non-hydrostatic pressure terms to avoid the restrictive time step that would give an explicit scheme. Various numerical methods can be found in the literature for approximating hyperbolic-elliptic dispersive systems: the pressure-correction-based methods (e.g., [26,81] in the framework of Navier-Stokes equations), semi-implicit methods on staggered meshes (see, e.g., [4,24,50,82]), or more recent trends based on obtaining a hyperbolic relaxation system that can be discretized with explicit methods (see [20,27,32,34,47,49]). However, such non-hydrostatic models widely developed for water wave propagation, have been poorly applied to granular flows. The non-hydrostatic contribution associated to the acceleration in the direction normal to the slope described above is different from other non-hydrostatic effects such as those arising from the viscous stress tensor involved in the pressure [13], which are always ignored in granular flow models owing to their complexity. Recently, a shallow model for dry granular flows with a weakly non-hydrostatic pressure was introduced in [42]. The pressure is called weakly non-hydrostatic because only the contribution of the normal acceleration is considered in the pressure and not the ones coming from the stress tensor. This is a strong limitation since these last terms are of first order in the dimensional analysis based on the shallow parameter εwhile those associated to the normal velocity are of second order. This choice relies only on the difficulty to handle the stress contribution from a numerical point of view. Despite this limitation, very promising results were obtained. In particular, it was possible to recover the parabolic shape of the evolution of the front velocity with time while hydrostatic models failed to reproduce the initial front acceleration. In that work, the influence of the coordinate system (local or Cartesian) was studied, showing that the non-hydrostatic model in Cartesian coordinates produces roughly similar deposits than the hydrostatic model in local coordinates. This was partially suggested previously by [29], where they considered a Cartesian model with a correction related to the vertical acceleration. Concerning non-hydrostatic models for multilayer systems, a hierarchy of such models for Euler equations was introduced in [40], all of them verifying a dissipative energy balance. They are denoted as LDNHkmodels, where k is the degree of approximation of the vertical velocity. This work presents the extension of the LDNH0model to granular flows with the µ(I)-rheology, obtaining a multilayer model with weakly non-hydrostatic pressure. It can be seen as the extension of the shallow model in [42] to the multilayer framework, but also as the extension of the multilayer model [36] to the non-hydrostatic framework. In particular, after the promising results of the one-layer model with a weakly non-hydrostatic pressure [42], it is a natural question whether its multilayer version could improve the modelling of granular flows, adding the advantages of improving the in-depth description of the flow shown in [36]. Our goal is to present this weakly non-hydrostatic multilayer model. It is expected that this model recovers the good results of both previous models [36,42], and solves some of the drawbacks of each of them, being cheaper computationally than a fully non-hydrostatic model. The multilayer weakly non-hydrostatic µ(I)-model proposed in the paper verifies a dissipative energy balance. A finite volume discretization is also introduced, being well-balanced for steady solutions at rest with non-constant free surface. As expected, the prize to pay is the increase of the computational effort, mainly due to the need to solve the non-hydrostatic pressure in multilayer systems, which is made via an iterative linear system solver. Concerning that, we propose a simple strategy that allows us to notably reduce the computational cost. This paper is organized as follows. In Section 2the initial system is stated and the dimensional analysis is performed. Section 3is devoted to the multilayer approach of the system, presenting the final model. The numerical approximation is described in Section 4and some numerical tests are presented in Section 5. Concretely, we study the influence of the coordinate system in the multilayer case, and compare with laboratory-scale granular collapse experiments. We also perform convergence and well-balanced tests, and show the efficiency of the proposed strategy to speed up the iterative solver for the linear system associated to the pressure unknowns. Finally, some conclusions are presented in Section 6. 2 Derivation of a non-hydrostatic µ(I)-model In this section we derive the non-hydrostatic multilayer model that we propose. It is based on the dimensional analysis made in [42], and the non-hydrostatic extension for multilayer systems introduced in [40]. Firstly, we write the governing equations and the non-dimensional 2D Navier-Stokes system, under the assumption of a weakly non-hydrostatic pressure. 3 2.1 Governing equations Let us briefly summarize the governing equations describing these flows. We consider a 2D granular mass with constant density ρ∈Rand velocity u∈R2, whose dynamics is described by the 2D Navier-Stokes system (∇·u= 0, ρ(∂tu+ (u·∇)u) = ∇·σ+ρg, where σ=−pI+τis the total stress tensor, with p∈Rthe total (hydrostatic and non-hydrostatic) pressure and τthe deviatoric tensor. This tensor is defined by (1), which can be written in the regime when the material flows (∥D∥ = 0), as τ=ηD(u), where η∈Ris the viscosity coefficient and the stain-rate tensor D(u) = 1 2∇u+ (∇u)′. In the µ(I)-rheology, introduced in [53], the viscosity is defined as (see details in [55,36] among others) η=µ(I)p p||D(u)||2+δ2, with δ > 0 a regularization parameter (see e.g. [45,66]), ||D|| =√0.5D:Dand µ(I) the friction coefficient defined by µ(I) = µs+µ2−µs I0+II, where I0and µ2> µsare constant parameters depending on the material, and I=2ds||D(u)|| pp/ρs is the inertial number, with dsand ρsthe grain diameter and density, respectively. The grain density and the apparent flow density (ρ) are related by ρ=φsρswith φsthe solid volume fraction. As boundary conditions, similarly to what is done in [36,42], we consider: •At the free surface: the usual kinematic condition and normal stress balance Nt+u·nb+h= 0,and p= 0, with Nt,nb+hthe downward time-space normal vector to the free surface. •At the bottom: the non-penetration condition and a Coulomb friction law u·nb= 0,and σ nb−σ nb·nbnb=−µ(I)pu |u|,0′ , with nbthe downward space vector normal to the bottom and uthe first component of u. 2.2 Initial system and dimensional analysis Let us start by establishing the notation used here. We consider local (or tilted) coordinates, as usually done in granular flow models, that refer to a fixed inclined plane (a straight line in the 1D x−zcase). Concretely, we consider a plane given by eb(x)=(xend −x) tan θ, with (x, z) the Cartesian coordinates, xend the final point of the domain, and θ > 0 a positive slope angle. We adopt here the usual convention in geophysical flows, that considers positive angles θfor negative slopes. Let (X, Z) be the local coordinates, measured in the direction along and normal to the inclined plane eb, respectively, and u= (u, w) the velocity vector, whose components are the downslope and normal velocities. We shall also consider a local bottom, b(X), as the deviation from the reference plane, and h(X) the distance from b(X) to the free surface, both measured in the normal direction (see Figure 1). 4 Figure 1: Sketch of the Cartesian (blue) and local (red) reference system. The dimensional analysis performed here is the same as in [42] (see also [36,46] among others). Let us remind it for the purpose of completeness. Defining H, L and Uthe characteristic height, length and velocity, we define the usual shallowness parameter ε=H/L, and denoting with tildes (e·) the non-dimensional variables, we obtain (X, Z, t)=(Le X, H e Z, (L/U)et), (u, w)=(Ueu, εU ew), h=Heh, ρ =ρ0eρ, p =ρ0U2ep, (τXX , τXZ, τZZ ) = ρ0U2(εgτXX ,gτXZ , εgτZZ), with gτXX =eη∂ e Xeu, gτXZ =eη 2∂e Zeu+ε2∂e Xew,gτZZ =eη∂e Zew, and therefore e D(e u) = U H 1 2 2ε∂ e Xeu ∂e Zeu+ε2∂e Xew ∂e Zeu+ε2∂e Xew2ε∂e Zew . Note also that ∥e D(e u)∥=∂e Zeu/2 + O(ε). Finally, the Froude number is defined as Fr =U/√gH cos θ. Following the asymptotic expansion performed in [42] for the pressure, and in particular for the non-hydrostatic counterpart, we write ep=ρ0 Fr2eb+eh−e Z+εg pnh 1+ε2g pnh,(2) where the pressure is divided into a hydrostatic contribution, and the first and second order non-hydrostatic contributions g pnh 1,g pnh. In that work, looking at the non-dimensional vertical momentum equation, these contributions are identified to the stress tensor contributions and the normal acceleration, respectively. There, first order terms were neglected whereas second order terms were kept. The reason is first that dealing with the terms coming from the stress tensor is notably more difficult, especially from the numerical point of view. Furthermore, including the normal acceleration (in this sense, the obtained model is called weakly non-hydrostatic) in [42] allowed to obtain very promising results for the depth-averaged (e.g. one-layer) model. It is therefore interesting to extend this weakly non-hydrostatic model to its multilayer version. To this aim, we follow the same approach here. Taking into account (2), the non-dimensional 2D Navier-Stokes system reads (tildes are dropped for the sake of simplicity)                ∂Xu+∂Zw= 0, ρ0∂tu+u∂Xu+w∂Zu+∂Xp=1 ε ρ0 Fr2tan θ+ε∂XτXX +1 ε∂ZτXZ , ρ0ε2∂tw+u ∂Xw+w ∂Zw+ε∂Zpnh 1+ε2∂Zpnh =ε∂XτZX +ε∂ZτZZ . (3a) (3b) (3c) 5 This is the system that will be layer-averaged using the multilayer approach. Before, boundary conditions are written in non-dimensional form as ∂th+u|b+h∂X(b+h)−w|b+h= 0, pnh |b+h= 0, at the free surface, and u|b∂Xb−w|b= 0,η 2∂Zu|b =µ(I|b)p|b u|b u|b, at the bottom. In next section, we develop the multilayer discretization and layer-averaging of system (3). 3 Multilayer system with the µ(I)-rheology and a weakly non-hydrostatic pressure In this subsection we apply the multilayer discretization introduced in [39] to system (3). We first recall the main aspects of this approach and we focus later on the terms that differ from the previous work [36], where all the details of the hydrostatic case are described. Figure 2: Sketch of the multilayer discretization. We first consider (lα)1≤α≤Nstrictly positive coefficients verifying PN α=1 lα= 1. Then, for a positive height h(t, X) we consider Nshallow layers whose heights are hα(t, X) = lαh(t, X) (see Figure 2). Notice that then PN α=1 hα=hholds. Considering now the levels Zα+1/2=b+Pα β=1 hβ, we denote the layers by Ωα(t, X) = {(t, X, Z) : Zα−1/2≤Z≤Zα+1/2},for α= 1, . . . , N. Note that b=Z1/2<··· < Zα+1/2<··· < ZN+1/2=b+h is a partition of [b, b +h]. As notation, we define fα=1 hαZZα+1/2 Zα−1/2 f(t, X, Z)dZ, the averaged value in the layer Ωαof the function f. In order to deal with possible discontinuities of facross the interface defined by Z=Zα+1/2, we define f− α+1/2= lim Z→Zα+1/2 Z<Zα+1/2 f|Ωα(t), f+ α+1/2= lim Z→Zα+1/2 Z>Zα+1/2 f|Ωα+1(t), 6 and [[f]]α+1/2=f+ α+1/2−f− α+1/2denotes the jump of facross the interface. In case of a continuous function, then [[f]]α+1/2= 0 and therefore f+ α+1/2=f− α+1/2=fα+1/2. In order to obtain the multilayer discretization of system (3), some assumptions on the unknowns (u, w, q) are needed. Then, the mass and momentum equations are integrated on each layer Ωα, leading to 3 layer-integrated equations for each layer. In this paper, we present a generalization of the model denoted by LDNH0in [40]. In this model, the two components of the velocity uand ware approximated by constant functions in each layer: u(t, X, Z) = N X α=1 uα(t, X) 1 Ωα(Z), w(t, X, Z) = N X α=1 wα(t, X) 1 Ωα(Z), and the non-hydrostatic pressure pnh(t, X, Z), pnh |Ωα∈P1by a continuous linear function. Note that in this case: pnh α(t, X) = pnh(t, X, Zα) = pnh α+1/2(t, X) + pnh α−1/2(t, X) 2(4) where Zα=Zα−1/2+hα/2 is the midpoint of the layer Ωα. Actually, for each layer we have pnh(Z) = pnh α+pnh α+1/2−pnh α−1/2 hα (Z−Zα),for Z∈[Zα−1/2, Zα+1/2]. Under these assumptions and by appropriately defining the jump conditions across the interfaces, the LDNH0model for the Euler equation is obtained in [40]. Here, the main novelty is to deal with the viscous term ∂ZτXZ , which involves a complex viscosity, as detailed in the following sections. 3.1 Jump conditions at the interfaces From an asymptotic analysis, we recover the following non-dimensional conditions when imposing the normal flux jump condition for the mass and the momentum equations (all the details are presented in [36]): The jump condition associated with the mass conservation equation leads to the definition of the mass transference term at the interface Z=Zα+1/2, denoted by Gα+1/2:= G+ α+1/2=G− α+1/2(corresponding to −Γα+1/2in [40]), where G± α+1/2=∂tZα+1/2+u± α+1/2∂XZα+1/2−w± α+1/2. Morover, from that condition we recover [[u]]α+1/2·nα+1/2= 0, where nα+1/2=∂XZα+1/2,−1is the downward normal vector at the interface. From the jump condition for the momentum equation, it follows that ε2τ± XX,α+1/2∂XZα+1/2−τ± XZ,α+1/2=ε2eτXX,α+1/2∂XZα+1/2−eτXZ,α+1/2±1 2ρ0εGα+1/2u+ α+1/2−u− α+1/2, τ± ZX,α+1/2∂XZα+1/2−τ± ZZ,α+1/2=eτZX,α+1/2∂XZα+1/2−eτZZ,α+1/2±1 2ρ0εGα+1/2w+ α+1/2−w− α+1/2, where we recall that τXX =η∂Xu,τXZ =1 2η∂Zu+ε2∂Xw,τZZ =η∂Zw, and eτis a consistent approximation of τat the Z=Zα+1/2(see [36]). 3.2 Layer-averaging procedure In this subsection we detail the layer-averaging of system (3) in the framework of the multilayer approach. We tackle the mass and the momentum conservation equation separately. 7 3.2.1 Mass conservation equation The mass conservation equation (3a) is integrated along the normal direction in each layer, leading to 0 = ZZα+1/2 Zα−1/2 (∂Xu+∂Zw)dZ =∂X(hαuα)−u− α+1/2∂XZα+1/2−w− α+1/2+u+ α−1/2∂XZα−1/2−w+ α−1/2, and therefore lα(∂th+∂X(huα)) = Gα+1/2−Gα−1/2,for α= 1, . . . , N. Notice that combining previous equations, the explicit formula for Gα+1/2is found (see [2]) Gα+1/2= α X β=1 lβ∂X(h(uβ−u)) ,with u= N X γ=1 lγuγ.(5) In particular, by summing up all the equations we have the mass conservation equation ∂th+∂X(hu)=0,(6) where G1/2=GN+1/2= 0 have been assumed as boundary conditions. 3.2.2 Horizontal momentum conservation equation We split in this case equation (3b) in convective, pressure, and viscous terms. The usual integration procedure gives ZZα+1/2 Zα−1/2 (∂tu+u∂Xu+w∂Zu)dZ =∂t(hαuα) + ∂Xhαu2 α−u− α+1/2Gα+1/2+u+ α−1/2Gα−1/2. Note that in the present approach (LDNH0model) it holds that u− α+1/2=u+ α−1/2=uα, since u|Ωα∈P0. Pressure terms are now integrated within the layer in the normal direction, leading to ZZα+1/2 Zα−1/2 ∂Xρ0 Fr2eb+b+h−z+εpnh 1+ε2pnhdz =hα ρ0 Fr2∂Xeb+b+h +∂Xεhαpnh 1,α +ε2hαpnh α−εpnh 1+ε2pnh|Zα+1/2 ∂XZα+1/2+εpnh 1+ε2pnh|Zα−1/2 ∂XZα−1/2, where we have used eb(X) = (Xend −X) sin θ. Finally, viscous terms are integrated, which yields ZZα+1/2 Zα−1/2ε∂XτXX +1 ε∂ZτXZ dz =ε∂X ZZα+1/2 Zα−1/2 η∂Xu dz!+1 εε2τ+ XX,α−1/2∂XZα−1/2−τ+ XZ,α−1/2−1 εε2τ− XX,α+1/2−τ− XZ,α+1/2 =ε∂X ZZα+1/2 Zα−1/2 η∂Xu dz!+εeτXX,α−1/2∂XZα+1/2−1 εeτXZ,α−1/2−εeτXX,α+1/2∂XZα+1/2−1 εeτXZ,α+1/2 +1 2ρ0Gα−1/2u+ α−1/2−u− α−1/2+1 2ρ0Gα+1/2u+ α+1/2−u− α+1/2. Defining now ˜ Kα+1/2=εeτXX,α+1/2∂XZα+1/2−1 εeτXZ,α+1/2and UH Zα+1/2=uα+1 −uα hα+1/2 , we obtain that ˜ Kα+1/2=Kα+1/2+O(ε),where Kα+1/2=−1 ε ηα+1/2 2UH Zα+1/2. 8 We can collect finally all the terms, keep terms up to first order and terms of order ε2, go to the dimensional variables, and obtain the horizontal momentum equation ρ∂t(hαuα) + ∂Xhαu2 α+gcos θhα∂Xeb+b+h+∂Xhαpnh α=Kα−1/2−Kα+1/2 +pnh α+1/2∂XZα+1/2−pnh α−1/2∂XZα−1/2+uα+1/2ρGα+1/2−uα−1/2ρGα−1/2,(7) for α= 1, . . . , N and uα+1/2=uα+uα+1 2. 3.2.3 Vertical momentum conservation equation Following the same procedure for equation (3c), we obtain the momentum conservation equation in the normal direction for each layer. Let us just detail the viscous terms. We have ZZα+1/2 Zα−1/2 ε(∂XτZX +∂ZτZZ )dz =ε∂X ZZα+1/2 Zα−1/2 η∂Zu+ε2∂Xwdz!+ετ+ ZX,α−1/2∂XZα−1/2−τ+ ZZ,α−1/2−ετ− ZX,α+1/2∂XZα+1/2−τ− ZZ,α+1/2 =ε∂X ZZα+1/2 Zα−1/2 η∂Zu+ε2∂Xwdz!+εeτZX,α−1/2∂XZα−1/2−eτZZ,α−1/2−εeτZX,α+1/2∂XZα+1/2−eτZZ,α+1/2 +1 2ρ0ε2Gα−1/2w+ α−1/2−w− α−1/2+1 2ρ0ε2Gα+1/2w+ α+1/2−w− α+1/2. It is observed that ε∂X ZZα+1/2 Zα−1/2 η∂Zu+ε2∂Xwdz!+εeτZX,α−1/2∂XZα−1/2−eτZZ,α−1/2−εeτZX,α+1/2∂XZα+1/2−eτZZ,α+1/2 is defined only by terms of order O(ε) and O(ε3). Then, collecting all the terms and keeping also terms of order ε2, the vertical momentum equation in dimensional variables is ρ(∂t(hαwα) + ∂X(hαuαwα)) = pnh α−1/2−pnh α+1/2+wα+1/2ρGα+1/2−wα−1/2ρGα−1/2,(8) for α= 1, . . . , N, where wα+1/2= (w+ α+1/2+w− α+1/2)/2. Notice that here first order term in εhave been neglected whereas second order terms are kept, similarly to the hypothesis made for the pressure (2). On the one hand, the discretization of the neglected first order terms by non-hydrostatic multilayer models has not been tackled yet in the literature. It entails important difficulties even in the case of a hydrostatic pressure (see [18]). On the other hand, including second order terms gave promising results in [42]. 3.3 Final multilayer non-hydrostatic model Notice that we have 3N+1 unknowns (h, uα, wα, pnh α) and 2N+1 equations till now. Then, the incompressibility condition is used in each layer. Concretely, we integrate the incompressibility equation (3a) for Z∈[Zα−1/2, Zα] and, taking into 9 Step 3: Non-hydrostatic pressure correction In the last step, the non-hydrostatic effects are added using the momentum equations (9b)-(9c) together with the incompressibility constraints (9d). Thus, we implicitly solve the following system: (∂tU=T(U,P, ∂XU, ∂XP, ∂XZ), C(U, ∂XU, ∂XZ) = 0,(26) where Tis evaluated at time tn+1 corresponding to the last term in (17). As usual, we will use an implicit projection correction method that will lead us to solve a linear system. Let us recall the constraint at the layer αby Eα= 0,∀α= 1, . . . , N, where Eα=wα−uα∂XZα+ α−1 X β=1 ∂X(lβhuβ) + ∂X(lαhuα) 2. As it can be seen, the imposition of the constraint for α= 2, . . . , N involves variables from the layer β= 1 to β=α−1. To simplify the set of constraints (9d), and to improve the diagonal dominance of the linear operator to be inverted, we equivalently consider Nα= 0,∀α= 1, . . . , N, where N1=E1,N2=E2−E1, . . . , Nα=EN−EN−1. For clarity, we will assume in what follows that lα= 1/N. Then, we define the vector Cin (13b) by C(U, ∂XU, ∂XZ)=2hN (N1, . . . , NN)′,(27) that reads C=       2Nhw1−2Nhu1∂xZ1+h∂x(hu1) 2N(hw2−hw1)−2N(hu2∂xZ2−hu1∂xZ1) + h∂xhu2+h∂xhu1 . . . 2N(hwN−hwN−1)−2N(huN∂xZN−huN−1∂xZN−1) + h∂xhuN+h∂xhuN−1        ,(28) where we recall that Zα=Zα+1/2−1 2Nh. We consider now a projection method, taking into account the fluctuation e Pn+1 defined in (15). To do that, we will impose the incompressibility constraints on the staggered mesh xi+1/2 CUn+1 i+1/2,(∂XUn+1)i+1/2,(∂XZ)i+1/2=0,(29) where Un+1 i+1/2=Un+2/3 i+1/2+ ∆tTUn+1 i+1/2,e Pn+1 i+1/2,(∂XUn+1)i+1/2,(∂Xe Pn+1)i+1/2,(∂XZ)i+1/2,(30) and the conserved variables and their derivatives on the staggered mesh are approximated from the averaged values Un+2/3 i+1/2=Un+2/3 i+Un+2/3 i+1 2,Un+1 i+1/2=Un+1 i+Un+1 i+1 2,(∂XUn+1)i+1/2=Un+1 i+1 −Un+1 i ∆x. Similarly, the component (∂XZα−1/2) of (∂XZ) is given by (∂XZα−1/2)i+1/2=h+ i+1/2−h− i+1/2−(hi+1 −hi) ∆x+α−1 N·hi+1 −hi ∆x, where the first addend is an approximation of ∂xbusing the hydrostatic reconstruction (20). Finally, the derivatives of the fluctuation of the non-hydrostatic pressures are approximated by (∂Xe Pn+1)i+1/2=e Pn+1 i+3/2−e Pn+1 i−1/2 2∆x. 16 Remark 3 Notice that the factor hin the definition of Cin (27)leads to a rewriting of the incompressibility constraints in terms of solely conserved variables (see (28)). That gives a more efficient and stable numerical implementation since no division by hmust be performed when computing C. Now, putting Eqs. (30) into (29) yields a linear system A·P=1 ∆tR,(31) where Pis a vector that contains a rearrangement of the non-hydrostatic pressure fluctuations unknowns P=   P1 . . . PN   ,Pα=     epnh,n+1 α,1/2 . . . epnh,n+1 α,Nx+1/2      , α = 1, . . . , N, Ais a block tridiagonal matrix of dimension N·(Nx+ 1) ×N·(Nx+ 1) given by A=       M1,1M1,2 M2,1M2,2M2,3 M3,2 ...... ......MN−1,N MN,N−1MN,N        ,(32) being Mi,j tridiagonal matrices of dimension (Nx+1)×(Nx+1) (see details in Appendix A). Finally, Ris the right-hand side of the linear system that is given by R=   R1 . . . RN   ,Rα=       CαUn+2/3 1/2,(∂XUn+2/3)1/2,(∂XZ)1/2 . . . CαUn+2/3 Nx+1/2,(∂XUn+2/3)Nx+1/2,(∂XZ)Nx+1/2        , α = 1, . . . , N, (33) where the sub-index αin Cαdenotes the α−th component of the vector (28). To numerically approximate the solution of the linear system (31), we consider the following block Gauss-Seidel solver: Algorithm 1: P(0) ←I.C. 2: s←1 3: Er ←0 4: do 5: for α= 1,2, . . . , N do 6: Solve the tridiagonal system: Mα,αP(s) α=1 ∆tRα−Mα,α−1P(s−1) α−1−Mα,α+1P(s−1) α+1 , 7: Er ←Er + ∆x||P(s−1) α−P(s) α||p 8: P(s−1) α←P(s) α 9: end for 10: while Er ≥ϵp tol. In line 1 of the algorithm, I.C. denotes the initial condition for the iterative solver. In practice, we take as initial condition the fluctuation at the previous time step e Pn=Pn−Pn−1,although other choices are possible. As stated in line 6 of 17 the algorithm, a tridiagonal linear system must be solved for each layer and at each loop of the iterative method, which is carried out with the efficient Thomas’ algorithm. The Error variable Er is used to compute the lpnorms for p= 1,or p= 2,used as the stopping criteria and given by ||P(s−1) α−P(s) α||p=1 N Nx X i=0 hi+1/2epnh,(s−1) α,i+1/2−epnh,(s) α,i+1/2 p,for p= 1,2. Alternative, the l∞norm can be considered, and therefore Er in line 7 is updated as Er ←max Er, max i=0,...,Nxnepnh,(s−1) α,i+1/2−epnh,(s) α,i+1/2o‘ whereas the stopping condition in line 10 is given by Er ≥ϵp tol.The influence of the choice of that different norms will be studied in Subsection 5.1.1. Finally, ϵp tol will be set in the numerical experiments. Remark 4 The proposed linear solver has shown the ability to converge for a few iterations (less than 62 for 5layers, as shown in Table 2) considering the strategy described in Subsection 4.1, that can be improved. One could use other linear solvers, such as LU factorizations or Krylov gradient-based methods. The first one is discarded since the matrix to be inverted (32)depends on the time step, increasing the computational effort. Krylov-based methods converge rapidly, especially if a convenient preconditioner is adopted and a matrix-free version of the algorithm to maintain low storage of the final implementation is considered (see, for instance [64]). In general, matrix-free linear solvers make the algorithm more convenient for future GPU implementations for 3D domains where memory is an important issue. The methodology employed here follows similar ideas presented in [31,59], where the resulting linear systems are solved using iterative Jacobi methods combined with a scheduled relaxation parallelized for 3D domains that show excellent computational performance. Some other authors have employed similar methods. For instance, in ([31,33]), different ad-hoc implementations of the Gauss-Seidel iterative procedures were employed for 2D domains. Once the non-hydrostatic pressure fluctuations e Pn+1 i+1/2have been computed, the discharges Un+1 ican be updated from Un+1 i=Un+2/3 i+ ∆tTUn+1 i,e Pn+1 i,(∂XUn+1)i,(∂Xe Pn+1)i,(∂XZ)i, where similarly, a second order point value approximation in the center of the cell will be used using (22). Finally, the non-hydrostatic unknowns are updated to Pn+1 i+1/2by (16). Remark 5 Since non-hydrostatic terms appear only in the momentum equations, the final numerical scheme is well-balanced for the water at rest solutions thanks to the hydrostatic reconstruction performed in step 1 and the definition (25), as we prove numerically in Subsection 5.3. It is positive preserving for the water height. It is a consequence of the nonnegativity of the HLL scheme for a flat bottom ([9]) and the numerical treatment for wet/dry areas [21] (see Remark 2). Moreover, the final numerical scheme is linearly L∞−stable under the usual CFL condition given in (14)(see [22]). Remark 6 Concerning boundary conditions, here we use standard methods for wall/free-outflow boundary conditions as in [31]. Concretely, we use a ghost-cell technique that duplicates the height and the pressures at the boundaries, and either duplicates or takes the opposite value for the discharges, depending on whether free-outflow or wall boundary conditions are considered. That is, a discretization of homogeneous Neumann boundary conditions are used for the height and pressures, and homogeneous Neumann (free-outflow case) or Dirichlet (wall case) boundary conditions are used for the discharges. In addition, the matrix A(32)for the pressure resolution is modified according to (37). 4.1 Strategy to speed up non-hydrostatic multilayer simulations As commented at the beginning of Section 4, the computational cost associated with non-hydrostatic multilayer systems becomes high, especially when increasing the number of vertical layers. Most of this effort is due to the third step in the scheme, that is, the elliptic problem solving the non-hydrostatic contributions. For this reason, any strategy to speed up this step is appreciated. We propose a straightforward strategy that allows us to significantly reduce the computational cost when using many layers. This strategy will be analyzed from the numerical point of view in Subsection 5.1. 18 The goal is to reduce the number of iterations in the non-hydrostatic iterative solver, or at least when considering a high number of vertical layers. We achieve it by improving the initial guess P(0), taking the solution provided by the linear system at the previous time step. The main idea is to solve the elliptic problem related to dispersive effects in a coarse vertical mesh, with a lower number of layers Nini guess < N, so that its solution is a better initial guess for the original system with Nlayers. In this way, we would solve the dispersive linear system twice: firstly with Nini guess layers, whose resolution is much faster than the original problem; secondly the full system with an improved initial guess, so that less iterations are needed to converge to the solution Pn+1. In order to be able to build the first linear system in the coarse vertical discretization, we assume a conformal mesh as follows. We define Ngr =N/Nini guess ∈Nrepresenting the number of layers that merge into a single layer. Then, the domain is divided along the normal direction into layers Ωβ,red, where Ωβ,red = Ngr [ j=1 ΩNgr(β−1)+jfor β= 1,...Nini guess. Now, we need to define auxiliary vector variables with size Nini guess for this domain subdivision. Those are given from the original vector variables by simply computing appropriate averaged values with the data of the merged layers. For instance, for the horizontal velocities, we define uβ,red =1 Ngr Ngr X j=1 uNgr(β−1)+jfor β= 1,...Nini guess, and analogously for all the variables involved in resolving the elliptic problem (30). Concretely, it is done to define Ured, Pred and e Pred, which are the input variables for the first stage where the linear system for the reduced elliptic operator is solved. Once we have the solution for the non-hydrostatic variable e Pred, it is used to generate the initial guess for the second stage, where the full linear system with Nlayers is solved. Thus, we define, for α= 1, . . . , N, epnh,(0) α+1/2=epnh β+1/2,red with 1 ≤β≤Nini guess such that Zβ−1/2≤Zα−1/2< Zβ+1/2, and epnh,(0) N+1/2= 0. These variables define the iterative algorithm’s initial condition P(0). It is expected that the number of iterations needed in this case for the full system is lower than the original one, where the solution at the previous time step is taken. Therefore, the elliptic problem’s computational time is also reduced, despite solving two linear systems. In Subsection 5.1.2 we will see that this strategy is effective for a granular collapse problem and different vertical discretizations. Let us remark that, although this strategy is effective for the numerical experiments performed here, other strategies could be used to speed up the iterative solver. 5 Numerical tests We perform here some numerical tests with the proposed multilayer model and compare the results to those obtained with precedent models, namely the multilayer hydrostatic [36] and the one-layer weakly non-hydrostatic [42] models. The goal is to show that the multilayer model with weakly non-hydrostatic pressure combines some of the strengths of these previous models. In particular, we will quantify the influence of the multilayer and of the non-hydrostatic approaches in several practical cases. It will be shown in Subsection 5.5. In particular, we recover important results of both models, although the prize is a higher computational cost. Before, concerning this computational effort, it will be evaluated in Subsection 5.1, together with the efficiency of the proposed numerical strategy to speed up the non-hydrostatic step of the numerical scheme presented in Section 4. Next, in Subsection 5.4, we will quantify the influence of the coordinate system (Cartesian or local) in the multilayer model in order confirm the results obtained with the one-layer model (see [42]). Note that this cannot be done by simulating granular collapse experiments since the associated initial conditions cannot be defined in a Cartesian reference system. We then move to the comparison with laboratory granular collapses experiments. In all the tests, we consider the material properties given in Table 1. In particular, the µ(I) coefficients are increased by 0.1 to account for sidewall friction, following the approach of [52]. For the numerical parameters, unless otherwise is specified, we consider ϵp tol = 10−6as tolerance for the iterative solver finding the non-hydrostatic pressure (see line 10 of the algorithm in Section 4), and δ= 10−5as regularization parameter in the viscosity coefficient (equation (10)). Note that this choice of δis in agreement with other works for viscoplastic fluids, as [41,25,58]. 19 µsµ2ds(mm) I0φs tan(25.5◦)≃0.48 tan(36◦)≃0.73 0.7 0.279 0.62 Table 1: Values of the rheological parameters used in the numerical tests. 5.1 Analysis of the computational efforts for non-hydrostatic simulations We investigate here the computational efforts associated to the non-hydrostatic multilayer model. The computation of the non-hydrostatic pressure described in Section 4is based on an iterative solver to solve a linear system whose size is N×(Nx+ 1), where Nis the number of layers in the normal direction and Nxthe number of cells in the horizontal discretization. Then, the computational time could grow dramatically when increasing the number of normal layers. This fact motivates the following attempt to reduce the computational effort. We first study the influence of changing the norm used in the stopping criterion of the iterative solver for the non-hydrostatic pressure, and secondly we evaluate the proposed strategy in Subsection 4.1 to speed up non-hydrostatic simulations. We consider here a computational domain [−0.2,3] mand 640 cells for the horizontal discretization, and the dam break problem given by the initial conditions eb(X) = (3 −X) sin θ, b(X)=0, h0(X) = 0.14 + hiif x≤0, hiotherwise,(34) with hi= 1.82 mm and θ= 22◦. We choose a high slope case since the multilayer framework is more appropriate for large values of the slopes, as it will be concluded in 5.5. The material parameters are in Table 1. A wall boundary condition is employed upstream while a free-outflow condition is used downstream. Finally, we use CFL = 0.7 and tf= 3 s. 5.1.1 Influence of the choice of the stopping criterion in the iterative solver As explained in Section 4, the computation of the non-hydrostatic pressure is based on an iterative solver to find the solution of a linear system of size N×(Nx+ 1). Therefore, the computational time is sensitive to the norm used in the stopping criterion for this solver. Layers l2- NºIter l2- tcomp l∞- NºIter l∞- tcomp 20 431 779.7 693 1096.7 10 145 84.8 211 197.7 8 101 58.8 142 114.3 5 46 18.1 62 36.4 2 10 4.2 14 5.2 Table 2: Computational times tcomp (s) and maximum number of iterations made in a single time step by the iterative solver for different number of normal layers Nand norms l2,l∞. We have measured the computational time and the number of iterations made by the solver for different number of vertical layers. Table 2shows the computational times and the maximum number of iterations needed in a time step for N= 20,10,8,5,2 normal layers, and Figure 3shows the iterations needed at each time step for N= 5,10,20 layers (Figure 3a) and the computational times for both norms (Figure 3b). When using the l2norm the computational effort is lower than with the l∞norm for all cases, as expected. In particular, for N= 20, the simulation needs 46% less time. It is due to the fact that the number of iterations needed with the l∞norm is always greater than with the l2norm (see Figure 3a). In particular, we see the difference in the number of iterations needed when using the different norms in the stopping criterion, after a small time. Note that when using the l1norm, the results are similar to the l2norm. In addition, we have checked that there are no significant differences in the results when using the different norms in the iterative solver. Then, the l2-norm is used from now on in the simulations, also in subsections 5.4 and 5.5. Most of the computational cost is spent in the Step 3, that is, the elliptic problem related to the non-hydrostatic pressure. In Table 3we show the percentage of the computational time that the algorithm spends in this step and also for Step 2 (viscous terms) for different number of layers. We see that it is 28% for the one-layer model, but it increases quickly with the number of layers. For 5 layers it is 81%, and it reaches 98% for 20 layers. In the non-hydrostatic case, 20 0123 0 50 100 150 200 250 0123 0 200 400 600 800 0 2 5 8 10 15 20 1 2 3 5 13 18 Figure 3: (a): Number of iterations needed for N= 5,10,20 layers with the l2and l∞norms. (b): Computational time in minutes for different number of layers with norms l2, l∞. we see that Step 2 does not spend a significant amount of resources. We also see that Step 2 spends approximately 10% of the computational time for N= 1,2, ..., 20 in the hydrostatic case. Layers 1 2 5 10 20 tcomp Step 3/tcomp (NH) 28% 58% 81% 92% 98% tcomp Step 2/tcomp (NH) 5% 4% 2% 0.6% 0.2% tcomp Step 2/tcomp (H) 9% 11% 11% 10% 9% Table 3: Percentage of computational time spent by steps 2 (viscous terms resolution) and 3 (elliptic problem associated to the pressure) of the scheme in Section 4in the non-hydrostatic case (NH). We also show this percentage for Step 2 in the hydrostatic case (H). 5.1.2 Speed-up strategy for non-hydrostatic simulations We investigate in this subsection the efficiency of the algorithm to speed up the numerical scheme for non-hydrostatic simulations that has been proposed in Subsection 4.1. (a) N= 8 layers (b) N= 10 layers Nini guess - 4 2 tcomp (s) 47.3 44.6 77.46 Speed-up - 1.06 0.61 Nini guess - 5 2 tcomp (s) 84.8 56.5 128.2 Speed-up - 1.5 0.66 (c) N= 20 layers Nini guess - 10 5 4 2 tcomp (s) 594.2 183.2 166.8 215.2 564.6 Speed-up - 3.24 3.56 2.76 1.05 Table 4: Computational times and speed-ups for non-hydrostatic simulations with (a) 8 (b) 10, (c) 20 vertical layers, reducing to Nini guess layers for the first step of the algorithm to solve the non-hydrostatic pressure. 21 0 5 10 15 20 0 2 4 6 8 10 0 0.1 0.2 0.3 0.4 0.5 0 100 200 300 400 0.6 1 1.5 2 2.5 0 50 100 Figure 4: (a): Computational time in minutes for different number of layers used to the precomputing of the initial guess for the pressure solver for N= 8,10,20 layers. (b): Number of iterations needed by the iterative solver in the case N= 20 for the simulation without using a precomputed guess (FS 20 layers, dashed black line) and when an initial guess is precomputed with 10 layers (first step with reduced system with 10 layers - RS-20-10 - solid cyan line and second step with the full system - FS-20-10 - dashed red line). Dotted gray line is the total of iterations (step 1 and step 2) in the case of precomputing the initial guess. In Table 4we show the computational times when reducing to Nini guess layers in the first step of the algorithm to compute a good initial guess for the pressure iterative solver for simulations with 8, 10 and 20 layers. We see that the proposed algorithm is effective, especially when using a large number of normal layers. Actually, there is no significant gain in computational effort when using 8 layers. However, it is useful for 10 and 20 layers. When the vertical resolution is increased to 20 layers, the proposed strategy allows us to significantly reduce the computational cost with an appropriate choice of Nini guess. In this case, when reducing to Nini disp = 10 layers for the first step of the pressure computation, the speed-up is 3.24, which means that with the proposed strategy the scheme is faster by 69% than the scheme without using the precomputed initial guess. For Nini guess = 5 the speed-up is 3.56 (72% faster) and for Nini guess = 5 the speed-up is 2.76 (64% faster). When using Nini guess = 2 there is no significant gain in computational time (it is only 5% faster). However, this strategy is not effective for all the choices. For instance, when using Nini disp = 2 with 10 or 8 layers, the computational time increases. It is clear that in the case of 20 layers the optimal choice is to use Nini disp = 5. In Figure 4we study the case of 20 layers and Nini disp = 10 in more detail. Figure 4a shows the computational times needed for different choices of Nini disp, while in Figure 4b we see an analysis when Nini disp = 10 in terms of the number of iterations needed. The more expensive part of the scheme is the resolution of the non-hydrostatic pressure as it is shown in Table 3, where the iterative solver is used. We see that the fact of solving a first system with a reduced number of layers (10 layers in this case) gives us a good initial guess for the resolution of the full system with 20 layers. It makes that the scheme with the proposed strategy needs much less iterations to converge than the scheme without the precomputing of the initial guess. Actually, by solving directly the full system, the maximum number of iterations is 431, whereas with the proposed strategy, the maximum number of iterations needed for the full system is 71. In summary, the proposed strategy is effective, especially when using a large number of layers, and it allows us to significantly reduce the computational time. The key is to reduce the number of iterations needed to solve the full linear system, even making a previous resolution with a smaller linear system at any time step. 5.2 Convergence test We perform here a convergence test showing the first order accuracy of the proposed numerical method. We consider the same test configuration as in [42]. Concretely, we consider a granular mass at rest with a smooth free surface, which flows over a slope, and measure the order at an intermediate time t= 0.15 s. The computational domain is [−0.5,1.5] m, CFL = 0.7 and free-outflow boundary conditions are assumed. The bottom and initial conditions are given by eb(X) = (1.5−X) sin 20◦, b(X)=0, h0(X) = 0.02 + 0.4e−5(x−0.25)2. 22 The same material properties as in previous test are used (see Table 1). We measure here the order of accuracy for the height and the velocity field in both directions by increasing Nand Nxat the same time. That is, we refine the mesh in both directions (x, z). Here, L1absolute errors are computed, which are given by ∥h−href ∥1= ∆x M X i=1 hi−ehref i,∥v−vref ∥1= ∆x Nx X i=1 hi N N X α=1 vα,i −evref α,i ,(35) where href ,evref must be understood as the projections of the reference solutions onto the spaces of the local solutions. In previous equation, v(z) denotes the horizontal (u(z)) or vertical velocity (w(z)). Here the solution of the multilayer model with Nx= 6400 and N= 128 will be the reference solution. (Nx, N) Error hOrder hError u(z) Order u(z) Error w(z) Order w(z) (100,2) 1.59×10−2- 4.14×10−2- 1.59×10−2- (200,4) 1.22×10−20.39 2.62×10−20.66 1.22×10−20.68 (400,8) 6.40×10−30.93 1.16×10−21.18 6.40×10−31.29 (800,16) 3.06×10−31.06 5.58×10−31.05 3.06×10−31.22 (1600,32) 1.35×10−31.18 2.65×10−31.07 1.35×10−31.18 Table 5: L1errors and related orders for the height (h) and the velocity fields (u(z), w(z)) for the solution at time t= 0.15 s. Table 5shows the L1errors and numerical convergence rates for the height (h), horizontal velocity (u(z)) and vertical velocity (w(z)) at time t= 0.15. We measure the errors at a short time by two main reasons: first, solutions with shocks and overshooting may appear for longer times; second, the computational cost associated to the reference solution considered here (Nx, N) = (6400,128) is huge. We see that the proposed scheme is first order accurate for all the variables, as expected. 5.3 Well-balanced test We perform now a test to check the well-balancing of the scheme for non flat steady states defined by (12). Let us consider the domain [−3,1] mwith 300 cells and free-outflow boundary conditions. The initial conditions (see figure 5) are eb(X) = 0, b(X) = 2 X+ 4 −1 2+e−10(X+3/2)2+ 0.15 RAND(X), h0(X)=0.8−µs(X−2) −b(X), being RAND(X)∈[0,1] a random value generated at run time, and the flow is assumed at rest (hu)α= (hw)α= 0 for α= 1, ..., N. We take the rheological parameter in Table 1except for µs= tan(20◦). CFL = 0.7 is used here and the final time simulation is tfin = 1 s. -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 0 1 2 3 Figure 5: Bottom topography and initial free surface the well-balanced test. The L1errors for all the variables are measured, as in (35), where in this case the reference solution evref is the initial state. Table 6shows the maximum along the simulated time of these errors for different values of the regularization 23 parameter δ. We show these errors for the scheme described in Section 4with the definition of bhn+2/3 igiven by (25) (Table 6a) and by considering bhn+2/3 i=hn+2/3 i(Table 6b), in order to show its influence. We see here that the definition (25) is essential to get the well-balanced property of the scheme. As commented before, it is not possible to obtain zero velocity or null variation of the variables because of the regularization parameter. In fact, there is a relation between δ and the order of magnitude of the variations (and the residual velocities). Actually, they are of the same order. A similar behaviour was observed in [38] for a hydrostatic multilayer model for viscoplatic (Herschel-Bulkley) fluids. (a) WB: bhn+2/3 i(25) (b) No WB: bhn+2/3 i=hn+2/3 i Error hError u(z) Error w(z) δ= 10−55.34×10−55.38×10−54.21×10−5 δ= 10−10 5.34×10−10 5.39×10−10 4.21×10−10 δ= 10−15 2.99×10−16 5.39×10−15 4.21×10−15 Error hError u(z) Error w(z) δ= 10−55.22×10−41.20×10−37.42×10−4 δ= 10−10 1.76×10−41.93×10−43.48×10−4 δ= 10−15 4.62×10−56.06×10−51.87×10−4 Table 6: Maximum along the time of L1errors for the height (h) and the velocity fields (u(z), w(z)) with respect to the initial state by considering the definition of b hn+2/3 igiven by (a) eq. (25)(well-balanced scheme); (b) hn+2/3 i(no well-balanced scheme). Let us remark that, since in Step 3 of the scheme we solve the non-hydrostatic pressure, which in this test is of the order of δ, we must use ϵp tol < δ. Then, we use here ϵp tol =δ/100. Notice that we recover the expected well-balanced property of the scheme even in the limit case ∂Xeb+b+h=µs. In addition, if we change the slope of the initial condition by tan(20.01◦) (i.e., 0.01◦greater than the angle of repose), then we recover that the error between the initial and final height is of order 10−5for any value of δ= 10−5,10−12,10−15. That is, the material is not at rest in this case. We must also remark that we are able to use very small values of δin this test because it is a steady at rest solution. In general transient problem we use a large value (δ= 10−5) for stability and convergence reasons (see [58]). In the following, we move to tests evaluating the ability of the model to recover some important physical behaviours. We first study the influence of the coordinate system and then compare the our results with dam break laboratory experiments. 5.4 Influence of the coordinate system Let us recall the differences between the local and Cartesian models. In thin-layer landslide models, the downslope velocity should be computed in the direction parallel to the topography and curvature effects should be added (see [14,60,70]). A reasonable approach, with the advantage of being much simpler conceptually and from the practical point of view, is to take a reference plane with a slope θcorresponding to the mean slope of the topography, and then consider variations of the bottom with respect to this reference plane (see e.g. Figure 1). It is actually the local coordinate system considered in this work. In this approach, the downslope velocity is tangent to the reference plane eband the equations are depth-averaged in the direction normal to the plane. In this model, the flow is considered thin in the normal direction compared to its extension parallel to the reference plane. On the contrary, in Cartesian models, the velocity is computed in the horizontal x-direction, the flow is assumed thin in the vertical direction compared to the horizontal one, and depth-averaging is performed in the vertical direction. When the slope is not small, it is well known that these models are less accurate than local models. Actually, Cartesian models overestimate the velocity. Another difference between these reference systems in the particular case of non-hydrostatic models, is the influence of the bottom in the terms pnh α±1/2Zα±1/2in the horizontal momentum equations (9b). Concretely, Zα±1/2is defined in terms of the local bottom bottom band not ebin the local coordinates case, contrary to the Cartesian case. 5.4.1 Numerical set-up To compare the results in Cartesian and local coordinate systems, we consider a test with an initial granular mass having a smooth shape so that its thickness can be defined in both systems. We take the initial condition as in [42]. It consists of a granular mass flowing in a complex topography with a bump located in the middle of the domain. We consider a Cartesian domain x∈[−3,3] m and 640 nodes. The material is initially at rest, and we set a final time tf= 10 s and CFL = 0.5. In this case we use free-outflow conditions at both boundaries. Figure 6shows the initial conditions for the bottom and the mass thickness in both coordinates system. In the Cartesian case, they are bCart(x)=1−tanh(x) + κe−10(x−1)2+ 0.5e−10(x−3)2, 24 and hCart(x) = max(η(x)−bCart(x),0) + hi,with η(x) = y0−1 + e−0.5(x−x0)8,if x≤0, 0,otherwise, with κ= 0.3, hi= 10−4m, x0=−1.53 m and y0= 1.91 m. In that (Cartesian) case, we fix θ= 0◦and ebCart(x) = 0. Note also that a small thickness hihas been added in order to avoid dealing with wet/dry areas and remove the influence of the treatment for such cases. In order to define these functions in the local coordinate system, the reference plane ebloc =−0.7−tan θ(x−3), with θ= 25◦, is considered. Now the local bottom bloc(X) is the distance from ebloc(x) to the bottom measured in the direction normal to the plane ebloc(x). Analogously, hloc(X) is measured as the distance from bloc(X) to the free surface. Figure 6: Sketch of the initial conditions in both Cartesian and local coordinates. κ= 0.3. The variables defined in the Cartesian system are in blue and those in the Local system are in red. 5.4.2 Local/Cartesian and hydrostatic/non-hydrostatic model comparison The goal here is to study two main aspects. We will investigate (i) the influence of the Coordinate system, which should be smaller for non-hydrostatic (NH) models than for hydrostatic (H) ones and (ii) the similarity between the non-hydrostatic Cartesian and the hydrostatic local models. Concerning (ii), [29] suggested that a Cartesian model with a correction associated to the normal acceleration produces similar results than a hydrostatic local model. In [42] we found that it is only true for the final deposit but not at intermediate times. Figure 7shows the final deposits for the local and Cartesian models with a hydrostatic and weakly non-hydrostatic pressure, and 10 vertical layers. The slope angle of the bottom bCart is also represented and varies from −15◦to 50◦(we recall that positive angles correspond to negative slopes). We see that, qualitatively, the hydrostatic local model gives a deposit similar to the non-hydrostatic Cartesian model, except near the initial location of the mass. In Figure 8a we show the time evolution of the relative errors between the mass thicknesses simulated with these two models. We see that, for the multilayer cases (5, 10 and 20 layers), the error is around 20% at final times whereas higher errors are obtained at intermediate times. These results are in agreement with [42] and contradicts what was assumed in [29]. In addition, we see that the differences are smaller in the multilayer case than in the one-layer case. In Figure 8b we show the influence of the coordinate system (Cartesian versus Local) on the simulations for both hydrostatic and non-hydrostatic models. This comparison is performed for the one-layer and the multilayer (5, 10 and 20 layers) cases. We see that, as expected, non-hydrostatic models are less dependent on the coordinate system. Interestingly, the influence of the coordinate system is smaller in the multilayer case than in the one-layer case. Moreover, in the hydrostatic case, the differences strongly decrease when increasing with the number of layers, whereas for the non-hydrostatic case the results for 5, 10 and 20 layers are almost identical. In order to show how the differences in the models could affect practical applications, we consider a bottom with a higher bump in the middle of the domain (at x= 1). Namely, we use κ= 0.45 here. Figure 9shows the deposits for all models. Interestingly, the simulated mass does not overcome the bump for the NH-Local model, whereas the other models predict that a significant amount of material is flowing further down to the valley. Such situation is usual in real applications where hazard assessment requires to estimate if a given landslide will overcome a topography relief such as this bump and eventually reach a town located in a valley [68,71,69]. The NH-Local model would predict that such 25 W(cm) model µsµ2µw 10 NH-µwtan(20.9◦)≃0.38 tan(32.76◦)≃0.64 tan(10.5◦) 20 NH tan(23.4◦)≃0.43 tan(34.74◦)≃0.69 0 20 NH-µwtan(20.9◦)≃0.38 tan(32.76◦)≃0.64 tan(10.5◦) Table 8: Modified rheological parameters considered when adding sidewalls friction in the simulations. Data for NH model with W= 10 m are in Table 1. to right-hand side of equation (24). Notice that, in practice, just the diagonal element of the linear system defined by (24) is modified. The models without the correction in the coefficients (but adding the term (36)) are denoted by NH −µw, and the rheological parameters used are summarized in Table 8for the widths W= 10,20 cm. Figure 18 shows the mass thickness at an intermediate time and when the material is at rest (deposit)computed with the non-hydrostatic models with and without the correction of the friction coefficients, for slopes θ= 0◦,16◦,22◦and W= 10 cm. We see that for θ= 0◦the results are similar in both cases. However, for large values of the slope, in particular for θ= 22◦, the final deposits are dramatically overestimated without the extra 0.1 in the friction coefficients. Actually, for this slope the flow is not at rest yet. Note that in this case the slope is higher than the friction angle arctan(µs) = 20.9◦. As a consequence, for large values of the slope, it is not possible to remove the correction in the friction coefficients. Interestingly, for the intermediate time the results are similar in all cases. The models NH-µw slightly improve the results in the initial part of the horizontal domain. Figure 19 shows the front velocity in these cases in addition to the results for a channel width W= 20 cm. We see that for high slopes, friction including sidewall effects is not enough to stop the material, and for θ= 0◦there is no significant difference between the results with or without the sidewall friction term. Figure 15: Flow/no-flow interface at different times at θ= 9.78◦(left-hand side) and θ= 22◦(right-hand side) and hi= 0 mm for the multilayer models with non-hydrostatic (dashed green lines) and hydrostatic (dashed magenta lines) pressure. Orange dash-dotted lines are results in [63]. 32 Figure 16: Normal profiles of the downslope normalized velocity (10uα,i/(maxβ,j |uβ,j|),β= 1, ..., N, j = 1, ..., Nx) obtained with the NH-20 l. model for θ= 22◦and hi= 1.82 mm during granular collapse at different positions, at (a) t= 0.60 s and (b) t= 1.60 s. Dashed green lines are flow/no-flow interfaces. (a) θ= 0◦,hi= 1 mm (b) θ= 22◦,hi= 1.82 mm Model Err(%) Err(%) tcomp (tini) (tfin) (s) NH-20 l. - - 1096.7 (19.08 min) NH-10 l. 0.30 0.33 368.2 (6.1 min) NH-8 l. 0.45 0.47 197.7 (3.3 min) NH-5 l. 0.79 0.88 60.1 (1.0 min) NH-2 l. 1.68 1.93 14.2 NH-1 l. 3.05 2.30 6.4 H-20 l. 10.01 8.32 56.1 H-10 l. 10.41 8.58 28.8 H-1 l. 11.65 11.25 4.7 Model Err(%) Err(%) tcomp (tini) (tfin) (s) NH-20 l. - - 1517.4 (25.3 min) NH-10 l. 0.92 1.92 544.6 (9.1 min) NH-8 l. 1.31 2.87 290.6 (4.8 min) NH-5 l. 2.19 5.58 89.2 (1.5 min) NH-2 l. 4.70 14.18 21.0 NH-1 l. 6.47 22.80 9.7 H-20 l. 7.36 11.21 79.4 H-10 l. 8.23 12.90 41.0 H-1 l. 11.78 31.93 7.1 Table 9: l2-errors made in the height by the different models compared to the reference solution corresponding to the non-hydrostatic model with 20 layers at intermediate times tini = 0.24 s(θ= 0◦),0.32 s(θ= 22◦), and at final time tfin 1.5s(θ= 0◦),3s(θ= 22◦).tcomp is the computational time needed to reach tfin. 5.5.3 Quantification of model differences By comparing the NH-10 l., NH-1 l. and H-10 l. models, we have shown that the proposed NH-10 l. model could be seen as a good compromise between the NH-1 l. and H-10 l. models. However, the price to pay is a higher computational effort. We have also concluded that multilayer effects have more influence for large slopes (at least for current models), while non-hydrostatic effects are more important for small slopes. Let us now quantify the improvement made using the multilayer model with a different number of layers. In Table 9we give the difference between the mass thickness 33 0 5 10 15 0 0.5 1 1.5 0 5 10 15 0 0.5 1 1.5 0 5 10 15 0 0.5 1 1.5 0 5 10 15 0 0.5 1 1.5 0 5 10 15 0 0.5 1 1.5 Figure 17: Time evolution of the normalized velocity of the front computed with hydrostatic (H) and non-hydrostatic (NH) models, and experimental data (solid-circle blue lines) for hi= 0 mm. Here h0= 0.14 m, v0=√gcos θh0m/s and τc=ph0/(gcos θ)s. calculated with the hydrostatic and non-hydrostatic models at an intermediate time and for the mass deposit for granular collapse on θ= 0◦and θ= 22◦, taking as reference the solution of the non-hydrostatic model with 20 vertical layers. We see that for θ= 0◦, if one uses 2 layers in the non-hydrostatic model the error in the thickness is below 2%, whereas we need at least 10 layers for θ= 22◦to have the same result. It is clear then that the multilayer approach is needed for large slopes, but it can be simplified for short values of θif one does not need additional information on in-depth variations of the flow. One should note that this fact (the influence of multilayer effects) could vary if the approximation of the flow/no-flow interface is improved, since its size is currently overestimated (see Figure 15). Improving the approximation of the flow/no-flow interface would make the multilayer approach more decisive for small slope values, where a wide static layer of material is expected to be recovered. Let us remark that we have also performed this quantitative analysis in the case of a rigid bed (hi= 0 mm) and with a smooth initial condition instead of the abrupt granular collapse, obtaining similar results to Table 9. Notice also that these differences shows the layer-independence of the results for N= 10 for large and small slopes, that is the number of layers used in the comparisons with experimental data. Concerning the computational times shown in Table 9, the non-hydrostatic model is more expensive than the hydrostatic one, as expected. Namely, we pay the prize of solving the elliptic operator associated to the non-hydrostatic pressure, which is coupled for all cells and layers, as seen in Section 4and analyzed in Subsection 5.1. 6 Conclusions In this work, we have introduced a multilayer model including the µ(I)-rheology and a weakly non-hydrostatic pressure to simulate dry granular flows. It is called a weakly non-hydrostatic model, in the sense that after a dimensional analysis, non-hydrostatic first-order terms related to the rheology are neglected, whereas second order terms related to normal 34 Figure 18: Granular mass over a plane of slope θ= 0◦(hi= 0 mm), 16◦(hi= 1.4mm), 22◦(hi= 1.82 mm), at an intermediate time and when the mass has stopped (deposit) for the laboratory experiments (solid-circle blue line), and non-hydrostatic models with and without the sidewall friction term of equation (36)(see Tables 1and 8). Note that at θ= 22◦, the granular mass it is not at rest yet. acceleration are kept. A well-balanced numerical discretization is also proposed, based on a three-step splitting procedure to deal with viscosity effects and non-hydrostatic pressure. The scheme is well-balanced for those steady states at rest where the slope of the free surface is lower than the angle of repose of the material. This is achieved thanks to a hydrostatic reconstruction technique involving the Coulomb friction term at the bottom. Thus, the proposed model can be seen as an extension of previous works. Concretely, it generalizes the hydrostatic multilayer model [36] to the non-hydrostatic framework and the one-layer non-hydrostatic model [42] to the multilayer framework. Then, all the good results associated with the discretization in the direction normal to the slope in [36] are present in this model. In particular, we show the ability of the model to reproduce changes in the shape of the velocity profiles and to approximate the flow/no-flow interface, which is not possible with one-layer models. In addition to the multilayer approach, the main achievements of the one-layer non-hydrostatic model of [42] are also present in this model. In particular, the simulation of the mass thickness at short times is improved, and the acceleration/deceleration of the front velocity is recovered, contrary to hydrostatic models where the first acceleration phase is not captured. Concerning the influence of the coordinate system, we have reinforced two of the results given in [42]: (i) non-hydrostatic models are less dependent on the coordinate system than hydrostatic models. Moreover, this dependence is even smaller when a normal discretization is considered; (ii) the non-hydrostatic model in Cartesian coordinates produces deposits similar to the hydrostatic model in local coordinates, but the simulation differ at intermediate times. Furthermore, the differences between the deposits are smaller when the multilayer approach is used. Thus, on the one hand, the novel model combines the strengths of both previous models, solving some of their drawbacks. On the other hand, the computational cost of this model is higher than for the previous ones, mainly due to the resolution of the linear system associated with the non-hydrostatic pressure, whose size is (M+ 1) ×N, being Mthe number of cells and Nthe number of normal layers. In order to speed up the scheme, we have proposed a first precomputing of the initial guess for the iterative solver. This strategy allows us to reduce around 70% the computational time when N= 20 layers. 35 In summary, the proposed model could be a reasonable compromise between efficiency and accuracy, recovering the good results of previous works at the same time. However, the results here suggest that an important improvement of these models would be to consider the viscous contributions of the momentum equation in the direction normal to the slope, which are neglected here. It corresponds to solving the full Navier-Stokes system in the layer-averaged (multilayer) framework. That is a difficult task that will be tackled in the future. Acknowledgements This work is partially supported by grants RTI2018-096064-B-C21 and RTI2018-096064-B-C22 funded by MCIN/AEI/10.13039/501100011033 and “ERDF A way of making Europe”, and by projects PID2020-114688RB-I00, PID2022-137637NB-C21 and PID2022-137637NB-C22, the ERC contract ERC-CG-2013-PE10-617472 SLIDEQUAKES, the DT-GEO Digital Twin project and the EnvSeis Doctoral Network. C. Escalante and J. Garres-D´ıaz have been also partially supported by project PID2020-114688RB-I00. J. Garres-D´ıaz has been also partially supported by the European Union - NextGenerationEU program. 0 5 10 15 0 0.5 1 1.5 0 5 10 15 0 0.5 1 1.5 0 5 10 15 0 0.5 1 1.5 0 5 10 15 0 0.5 1 1.5 0 5 10 15 0 0.5 1 1.5 0 5 10 15 0 0.5 1 1.5 Figure 19: Time evolution of the normalized velocity of the front computed with non-hydrostatic models with and without the correction in the µ(I)coefficients accounting for the sidewall friction, and experimental data (solid-circle blue lines), for W= 10,20 cm. Here h0= 0.14 m, v0=√gcos θh0m/s and τc=ph0/(gcos θ)s. 36 A Coefficients of the linear system associated to dispersive effects Here we describe the coefficients appearing in (32) where we will assume free-outflow boundary conditions for simplicity. Therefore, the right-hand side vector in (33) is modified as follows R=   R1 . . . RN   ,Rα=                    0 0 CαUn+2/3 5/2,(∂XUn+2/3)5/2,(∂XZ)5/2 . . . CαUn+2/3 Nx−5/2,(∂XUn+2/3)Nx−5/2,(∂XZ)Nx−5/2 0 0                    , α = 1, . . . , N. The tridiagonal matrices Mα,α−1, Mα,α,and Mα,α+1 for α= 2, . . . , N in (32) are described, for i, j = 3, . . . , Nx−2, by (Mα,α−1)i,j =                          h2 i+1/2 ∆x2 −hi+1/2 ξα,i+1/2+ξα+1/2,i+1/2 ∆x+N2,for j=i−1, −2 h2 i+1/2 ∆x2+ 4ξ2 α,i+1/2+ 2hi+1/2∂xξα,i+1/2+ 2N2,for j=i, h2 i+1/2 ∆x2+hi+1/2 ξα,i+1/2+ξα+1/2,i+1/2 ∆x+N2,for j=i+ 1, (Mα,α)i,j =                          2 h2 i+1/2 ∆x2 −hi+1/2 ξα,i+1/2+ξα+3/2,i+1/2−ξα+1,i+1/2−ξα−1/2,i+1/2 ∆x −2N2,for j=i−1, −4 h2 i+1/2 ∆x2 −4ξ2 α,i+1/2+ξ2 α+1,i+1/2−2hi+1/2∂xξα,i+1/2+∂xξα+1,i+1/2−4N2,for j=i, 2 h2 i+1/2 ∆x2+hi+1/2 ξα,i+1/2+ξα+3/2,i+1/2−ξα+1,i+1/2−ξα−1/2,i+1/2 ∆x −2N2,for j=i+ 1 (Mα,α+1)i,j =                          h2 i+1/2 ∆x2+hi+1/2 ξα+1,i+1/2+ξα+1/2,i+1/2 ∆x+N2,for j=i−1, −2 h2 i+1/2 ∆x2+ 4ξ2 α+1,i+1/2 −2hi+1/2∂xξα+1,i+1/2+ 2N2,for j=i, h2 i+1/2 ∆x2 −hi+1/2 ξα+1,i+1/2+ξα+1/2,i+1/2 ∆x+N2,for j=i+ 1, with MN,N+1 =0,and ξα,i+1/2= (α−1/2)hi+1 −hi ∆x+Nh+ i+1/2−h− i+1/2−(hi+1 −hi) ∆x, that is proportional to the derivative of the midpoint of the layer Zα,and we approximate their second order derivatives as follows ∂xξα,1+1/2=1 ∆xminmod ξα,i+3/2−ξα,i+1/2, ξα,i+1/2−ξα,i−1/2,1 2ξα,i+3/2−ξα,i−1/2, minmod being the usual slope limiter function. 37 The special cases of M1,1,and M1,2,corresponds to the discretization at the lowest layer, where the constraint (27) is simply given by 2hNN1= 2hNE1,and reads (M1,1)i,j =                          h2 i+1/2 ∆x2 −hi+1/2 ξα+3/2,i+1/2−ξα+1,i+1/2 ∆x −N2,for j=i−1, −2 h2 i+1/2 ∆x2 −4ξ2 α+1,i+1/2+ 2hi+1/2∂xξα+1,i+1/2−2N2,for j=i, h2 i+1/2 ∆x2+hi+1/2 ξα+3/2,i+1/2−ξα+1,i+1/2 ∆x −N2,for j=i+ 1 (M1,2)i,j =                          h2 i+1/2 ∆x2+hi+1/2 ξα+1,i+1/2+ξα+1/2,i+1/2 ∆x+N2,for j=i−1, −2 h2 i+1/2 ∆x2+ 2N2,for j=i, h2 i+1/2 ∆x2 −hi+1/2 ξα+1,i+1/2+ξα+1/2,i+1/2 ∆x+N2,for j=i+ 1. Regarding the horizontal wall/free-outflow boundary conditions along the independent variable x, we consider, for α= 1,2,...N, (Mα,α)1,1= (Mα,α)2,2= (Mα,α)Nx−1,Nx−1= (Mα,α)Nx,Nx= 1, (Mα,α)1,2= (Mα,α)2,3= (Mα,α)Nx−1,Nx−2= (Mα,α)Nx,Nx−1=−1,(37) and for the rest of matrices, Mα,α±1,they are given as usually when applying ghost cells techniques for all the variables. B Proof of Theorem 1 This appendix gives detailed proof of Theorem 1concerning the energy balance for model (9). Some of the steps of this proof are similar to the proof of the energy inequality verified by LDNH0model in [40]. However, we give here all the details for the sake of completeness. Proof: Using the horizontal momentum equation (9b) as well as the layer height evolution equation (9a), we obtain hα∂tuα+∂xu2 α 2+uα(Gα+1/2−Gα−1/2) + ∂xhαpnh α ρ+gcos θhα∂x(eb+b+h) =1 ρKα−1/2−Kα+1/2+pnh α+1/2 ρ∂xzα+1/2−pnh α−1/2 ρ∂xzα−1/2+uα+1/2Gα+1/2−uα−1/2Gα−1/2. (38) Now, we sum up (38) multiplied by uαwith the equation obtained by multiplying the mass conservation equation by u2 α/2 + gcos θ(eb+b+h), getting ∂thα u2 α 2+gcos θ(eb+b+h)∂thα+∂xhαuαu2 α 2+gcos θ(eb+b+h) −gcos θ(eb+b+h)Gα+1/2−Gα−1/2+1 ρuα∂xhαpnh α−1 ρuαpnh α+1/2∂xzα+1/2−pnh α−1/2∂xzα−1/2 −Gα+1/2u2 α 2−uαuα+1/2+Gα−1/2u2 α 2−uαuα−1/2+1 ρuαKα−1/2−Kα+1/2 (39) Similar steps are followed for the equations of the vertical velocity (9c), and using the mass equation, we get ∂thα w2 α 2+∂xhαuα w2 α 2+1 ρwαpnh α+1/2−pnh α−1/2= −Gα+1/2w2 α 2−wαewα+1/2+Gα−1/2w2 α 2−wαewα−1/2.(40) 38 By defining Eα=hαu2 α+w2 α 2+gcos θeb+b+h 2, and summing up equations (39) and (40), after some straightforward computations we obtain ∂tEα+∂xEαuα+gcos θhα h 2uα+1 ρPNH,α =MTα+1 ρuαKα−1/2−Kα+1/2+gcos θ 2(hα∂th−h∂thα),(41) where we have used that gcos θ(eb+b+h)∂thα=∂thαgcos θeb+b+h 2+gcos θ 2(h∂thα−hα∂th). In equation (41), PNH,α and MTαcollect the terms related to non-hydrostatic pressure and mass transference, respectively, in the layer Ωα. The non-hydrostatic contributions are PNH,α =uα∂xhαpnh α+pnh α+1/2−uα∂xzα+1/2+wα−pnh α−1/2−uα∂xzα−1/2+wα =∂xhαuαpnh α−pnh α+1/2 α X β=1 ∂x(hβuβ) + pnh α−1/2 α−1 X β=1 ∂x(hβuβ) where we have used restriction (4) and (9d). Notice that in the particular cases α= 1 and α=N, using the non-penetration condition and pnh |b+h= 0, we have PNH,1=∂xh1u1pnh 1−pnh 3/2∂x(h1u1), PNH,N =∂xhNuNpnh N+pnh N−1/2 N−1 X β=1 ∂x(hβuβ). It leads to a conservative contribution since the terms on pnh α±1/2vanish when summing up from α= 1, . . . , N: N X β=1 PNH,β = N X β=1 ∂xhβpnh βuβ. Concerning the mass transfer terms, they are MTα=gcos θeb+b+hGα+1/2−Gα−1/2+Gα+1/2 2uαuα+1 +wαwα+1−Gα−1/2 2uαuα−1+wαwα−1. It is easy to see that, since we have set G1/2=GN+1/2= 0, it holds N X β=1 MTβ= 0,and g 2 N X β=1 (hα∂tH−H∂thα)=0, and thus, summing up (41) in all the layers, the proof is finished. Finally, by using the definitions of Kα+1/2for α= 0, . . . , N, the viscous terms satisfy N X β=1 uαKβ−1/2−Kβ+1/2=u1K1/2−uNKN+1/2+ N−1 X β=1 Kβ+1/2(uβ+1 −uβ) =−µ(I1/2)p1/2 u2 1 |u1|− N−1 X β=1 ηβ+1/2 2 (uβ+1 −uβ)2 hβ+1/2 , that is then a dissipative contribution, as expected. □ 39 References [1] B. Andreotti, Y. Forterre, and O. Pouliquen. Granular Media: Between Fluid and Solid. Cambridge University Press, 2013. 1 [2] E. Audusse, M-O. Bristeau, and A. Decoene. Numerical simulations of 3D free surface flows by a multilayer Saint-Venant model. International Journal for Numerical Methods in Fluids, 56(3):331–350, 2008. 8 [3] E. Audusse, M.-O. Bristeau, B. Perthame, and J. Sainte-Marie. A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation. ESAIM: Mathematical Modelling and Numerical Analysis, 45(1):169–200, jun 2010. 2,11 [4] N. A¨ıssiouene, M.-O. Bristeau, E. Godlewski, A. Mangeney, C. Par´es Madro˜nal, and J. Sainte-Marie. A two-dimensional method for a family of dispersive shallow water models. The SMAI journal of computational mathematics, 6:187–226, September 2020. 3 [5] T. Barker and J. M. N. T. Gray. Partial regularisation of the incompressible µ(I)-rheology for granular flow. Journal of Fluid Mechanics, 828:5–32, aug 2017. 2 [6] T. Barker, M. Rauter, E. S. F. Maguire, C. G. Johnson, and J. M. N. T. Gray. Coupling rheology and segregation in granular flows. Journal of Fluid Mechanics, 909, dec 2020. 1,2 [7] T. Barker, D. G. Schaeffer, P. Bohorquez, and J. M. N. T. Gray. Well-posed and ill-posed behaviour of the µ(I)-rheology for granular flow. Journal of Fluid Mechanics, 779:794–818, 9 2015. 1 [8] T. Barker, D. G. Schaeffer, M. Shearer, and J. M. N. T. Gray. Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2201):20160846, may 2017. 2 [9] F. Bouchut. Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws. Birkh¨auser Basel, 2004. 14,18 [10] F. Bouchut, J. M. Delgado-S´anchez, E. D. Fern´andez-Nieto, A. Mangeney, and G. Narbona-Reina. A bed pressure correction of the friction term for depth-averaged granular flow models. Applied Mathematical Modelling, 106:627–658, jun 2022. 2,28 [11] F. Bouchut, E. D. Fern´andez-Nieto, E. H. Kon´e, A. Mangeney, and G. Narbona-Reina. Dilatancy in dry granular flows with a compressible µ(i) rheology. Journal of Computational Physics, 429:110013, mar 2021. 2 [12] F. Bouchut, E. D. Fern´andez-Nieto, A. Mangeney, and G. Narbona-Reina. A two-phase two-layer model for fluidized granular flows with dilatancy effects. Journal of Fluid Mechanics, 801:166–221, August 2016. 1,2 [13] F. Bouchut, I. Ionescu, and A. Mangeney. An analytic approach for the evolution of the static/flowing interface in viscoplastic granular flows. Communications in Mathematical Sciences, 14(8):2101–2126, 2016. 3 [14] F. Bouchut and M. Westdickenberg. Gravity driven shallow water models for arbitrary topography. Communications in Mathematical Sciences, 2(3):359–389, 09 2004. 24,26 [15] J. Boussinesq. Th´eorie des ondes et des remous qui se propagent le long d’un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond. J. Math. Pures Appl., pages 55–108, 1872. 2 [16] M. Bouzid, M. Trulsson, P. Claudin, E. Cl´ement, and B. Andreotti. Nonlocal rheology of granular flows across yield conditions. Physical Review Letters, 111(23), dec 2013. 2 [17] F. Boyer, ´ E. Guazzelli, and O. Pouliquen. Unifying suspension and granular rheology. Physical Review Letters, 107(18), oct 2011. 2 [18] M.-O. Bristeau, C. Guichard, B. di Martino, and J. Sainte-Marie. Layer-averaged euler and Navier-Stokes equations. Communications in Mathematical Sciences, 15(5):1221–1246, 2017. 9 40 [19] M. Brunet, L. Moretti, A. Le Friant, A. Mangeney, E. D. Fern´andez Nieto, and F. Bouchut. Numerical simulation of the 30–45 ka debris avalanche flow of Montagne Pel´ee volcano, Martinique: from volcano flank collapse to submarine emplacement. Natural Hazards, 87(2):1189–1222, apr 2017. 2 [20] S. Busto, M. Dumbser, C. Escalante, N. Favrie, and S. Gavrilyuk. On high order ADER Discontinuous Galerkin schemes for first order hyperbolic reformulations of nonlinear dispersive systems. Journal of Scientific Computing, 87(2), March 2021. 3 [21] M. J. Castro, J. M. Gonz´alez-Vida, and C. Par´es. Numerical treatment of wet/dry fronts in shallow flows with a modified Roe scheme. Mathematical Models and Methods in Applied Sciences, 16(06):897–931, jun 2006. 15,18 [22] M. J. Castro D´ıaz and E. D. Fern´andez-Nieto. A class of computationally fast first order finite volume solvers: PVM methods. SIAM Journal on Scientific Computing, 34(4):A2173–A2196, jan 2012. 14,18 [23] V. Casulli. A semi-implicit finite difference method for non-hydrostatic free-surface flows. International Journal for Numerical Methods in Fluids, 30(4):425–440, jun 1999. 2 [24] V. Casulli and P. Zanolli. Semi-implicit numerical modeling of nonhydrostatic free-surface flows for environmental problems. Mathematical and Computer Modelling, 36(9-10):1131–1149, dec 2002. 3 [25] J. Chauchat and M. M´edale. A three-dimensional numerical model for dense granular flows based on the µ(I)-rheology. Journal of Computational Physics, 256(0):696–712, 2014. 19 [26] A. J. Chorin. Numerical solution of the Navier-Stokes equations. Mathematics of Computation, 22(104):745–762, 1968. 3 [27] A. Dedner, F. Kemm, D. Kr¨oner, C.-D. Munz, T. Schnitzer, and M. Wesenberg. Hyperbolic divergence cleaning for the MHD equations. Journal of Computational Physics, 175(2):645–673, January 2002. 3 [28] R Delannay, A Valance, A Mangeney, O Roche, and P Richard. Granular and particle-laden flows: from laboratory experiments to field observations. Journal of Physics D: Applied Physics, 50(5):053001, 2017. 1 [29] R. P. Denlinger and R. M. Iverson. Granular avalanches across irregular three-dimensional terrain: 1. Theory and computation. Journal of Geophysical Research: Earth Surface, 109(F1), mar 2004. 3,25 [30] C. Escalante, T. Morales de Luna de Luna, and M. J. Castro. Non-hydrostatic pressure shallow flows: GPU implementation using finite volume and finite difference scheme. Applied Mathematics and Computation, 338(338):631–659, dec 2018. 3 [31] C. Escalante, E. D. Fern´andez-Nieto, T. Morales de Luna, and M. J. Castro. An efficient two-layer non-hydrostatic approach for dispersive water waves. Journal of Scientific Computing, 79(1):273–320, oct 2018. 18 [32] C. Escalante and T. Morales de Luna. A general non-hydrostatic hyperbolic formulation for Boussinesq dispersive shallow flows and its numerical approximation. Journal of Scientific Computing, 83(3), jun 2020. 3 [33] C. Escalante S´anchez, E. D. Fern´andez-Nieto, T. Morales de Luna, Y. Penel, and J. Sainte-Marie. Numerical simulations of a dispersive model approximating free-surface Euler equations. Journal of Scientific Computing, 89(3), oct 2021. 11,18 [34] N. Favrie and S. Gavrilyuk. A rapid numerical method for solving serre–green–naghdi equations describing long free surface gravity waves. Nonlinearity, 30(7):2718–2736, May 2017. 3 [35] E. D. Fern´andez-Nieto, F. Bouchut, D. Bresch, M. J. Castro D´ıaz, and A. Mangeney. A new Savage-Hutter type model for submarine avalanches and generated tsunami. Journal of Computational Physics, 227(16):7720–7754, 2008. 26 [36] E. D. Fern´andez-Nieto, J. Garres-D´ıaz, A. Mangeney, and G. Narbona-Reina. A multilayer shallow model for dry granular flows with the µ(I)-rheology: application to granular collapse on erodible beds. Journal of Fluid Mechanics, 798:643–681, jun 2016. 1,2,3,4,5,6,7,11,19,28,29,35 [37] E. D. Fern´andez-Nieto, J. Garres-D´ıaz, A. Mangeney, and G. Narbona-Reina. 2D granular flows with the µ(I) rheology and side walls friction: a well-balanced multilayer discretization. Journal of Computational Physics, 356:192–219, 2018. 2,14,15,31 41