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

Cipriano Escalante, Enrique Domingo Fernández-Nieto, José Garres-Díaz, Anne Mangeney

► To cite this version:

Cipriano Escalante, Enrique Domingo Fernández-Nieto, José Garres-Díaz, Anne Mangeney. Multilayer shallow model for dry granular flows with a weakly non-hydrostatic pressure. Journal of Scientific Computing, 2023, 96 (3), pp.88. 10.1007/s10915-023-02299-y. hal-03941686

HAL Id: hal-03941686

https://hal.science/hal-03941686v1

Submitted on 21 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multilayer shallow model for dry granular flows with a weakly non-hydrostatic pressure

C. Escalante,* J. Garres-Díaz,† E.D. Fernández-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ández-Nieto et al. [40], the multilayer model with $\mu(I)$ rheology Fernández-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 $\mu(I)$-rheology, introduced in [53, 54]. This rheology is based on a Drucker-Prager plasticity criterion, that defines the deviatoric tensor as

$$ \left\{ \begin{array}{ll} \tau = \frac{\mu(I)p}{\|D\|} D & \text{if } \|D\| \neq 0, \\ \|\tau\| \le \mu_s p & \text{if } \|D\| = 0, \end{array} \right. \qquad (1) $$

where $p$ is the pressure, $\mu_s$ is the tangent of the repose angle of the material, $D$ the strain-rate tensor and $\mu(I)$ is a strain-rate and pressure dependent friction coefficient, which will be defined later. It is known that the incompressible $\mu(I)$-rheology is ill-posed for low and high values of the inertial number $I$ (see [7]). Many works have been devoted to

<footnote>*Dpto. Matemática Aplicada. U. Málaga, 14071- Málaga, Spain (escalante@uma.es)</footnote>

<footnote>†Dpto. Matemáticas. Edificio Einstein - U. Córdoba. Campus de Rabanales, 14014-Córdoba, Spain (jgarres@uco.es)</footnote>

<footnote>‡Dpto. Matemática Aplicada I. ETS Arquitectura - U. Sevilla. Avda. Reina Mercedes S/N, 41012-Sevilla, Spain (edofer@us.es)</footnote>

<footnote>§Institut de Physique du Globe de Paris, Seismology team, U. Paris-Diderot, Sorbonne Paris Cité, 75238, Paris, France. ANGE team, CEREMA, INRIA, Lab. J. Louis Lions, 75252, Paris, France. Institut Universitaire de France (IUF), France, (mangeney@ipgp.fr)</footnote>
improve the well-posed character of the $\mu(I)$-rheology. A regularized incompressible $\mu(I)$ law, which improves the range of values for $I$ where it is well-posed, was proposed in [5]. This regularized expression was considered in [44, 6] among others. A compressible version of the $\mu(I)$, called $\mu(I)$, $\phi(I)$-rheology, has also been investigated (see e.g. [17, 77]), considering that the solid volume fraction $\phi$ is no more constant. Some authors suggested that this compressible version improves the well/ill-posed character of the $\mu(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 $\phi$ is also variable (depending on $I$ in the aerial case or on the viscous number $J$ in 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 $\mu(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 $\mu(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 $\mu(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 $\mu(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 $\mu(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
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 $\epsilon$ 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 LDNH$_k$ models, where k is the degree of approximation of the vertical velocity. This work presents the extension of the LDNH$_0$ model to granular flows with the $\mu(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 $\mu(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 2 the initial system is stated and the dimensional analysis is performed. Section 3 is devoted to the multilayer approach of the system, presenting the final model. The numerical approximation is described in Section 4 and 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 $\mu(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.
## 2.1 Governing equations

Let us briefly summarize the governing equations describing these flows. We consider a 2D granular mass with constant density $\rho \in \mathbb{R}$ and velocity $\mathbf{u} \in \mathbb{R}^2$, whose dynamics is described by the 2D Navier-Stokes system

$$ \left\{ \begin{array}{l} \nabla \cdot \mathbf{u} = 0, \\ \rho (\partial_t u + (\mathbf{u} \cdot \nabla) \mathbf{u}) = \nabla \cdot \boldsymbol{\sigma} + \rho g, \end{array} \right. $$

where $\boldsymbol{\sigma} = -p\mathbf{I} + \boldsymbol{\tau}$ is the total stress tensor, with $p \in \mathbb{R}$ the total (hydrostatic and non-hydrostatic) pressure and $\boldsymbol{\tau}$ the deviatoric tensor. This tensor is defined by (1), which can be written in the regime when the material flows ($\|\mathbf{D}\| \neq 0$), as

$$ \boldsymbol{\tau} = \eta \mathbf{D}(\mathbf{u}), $$

where $\eta \in \mathbb{R}$ is the viscosity coefficient and the stain-rate tensor $\mathbf{D}(\mathbf{u}) = \frac{1}{2}(\nabla \mathbf{u} + (\nabla \mathbf{u})')$. In the $\mu(I)$-rheology, introduced in [53], the viscosity is defined as (see details in [55, 36] among others)

$$ \eta = \frac{\mu(I)p}{\sqrt{\|\mathbf{D}(\mathbf{u})\|^2 + \delta^2}}, $$

with $\delta > 0$ a regularization parameter (see e.g. [45, 66]), $\|\mathbf{D}\| = \sqrt{0.5\mathbf{D}:\mathbf{D}}$ and $\mu(I)$ the friction coefficient defined by

$$ \mu(I) = \mu_s + \frac{\mu_2 - \mu_s}{I_0 + I} I, $$

where $I_0$ and $\mu_2 > \mu_s$ are constant parameters depending on the material, and

$$ I = \frac{2d_s \|\mathbf{D}(\mathbf{u})\|}{\sqrt{p/\rho_s}} $$

is the inertial number, with $d_s$ and $\rho_s$ the grain diameter and density, respectively. The grain density and the apparent flow density ($\rho$) are related by $\rho = \varphi_s\rho_s$ with $\varphi_s$ the 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

$$ N_t + \mathbf{u} \cdot \tilde{\mathbf{n}}^{b+h} = 0, \quad \text{and} \quad p = 0, $$

with $(N_t, \tilde{\mathbf{n}}^{b+h})$ the downward time-space normal vector to the free surface.

* At the bottom: the non-penetration condition and a Coulomb friction law

$$ \mathbf{u} \cdot \mathbf{n}^b = 0, \quad \text{and} \quad \boldsymbol{\sigma} \tilde{\mathbf{n}}^b - ((\boldsymbol{\sigma} \tilde{\mathbf{n}}^b) \cdot \tilde{\mathbf{n}}^b) \tilde{\mathbf{n}}^b = \left(-\mu(I) p \frac{u}{|u|}, 0\right)' $$

with $\tilde{\mathbf{n}}^b$ the downward space vector normal to the bottom and $u$ the first component of $\mathbf{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-z$ case). Concretely, we consider a plane given by $\tilde{b}(x) = (x_{end} - x) \tan \theta$, with $(x, z)$ the Cartesian coordinates, $x_{end}$ the final point of the domain, and $\theta > 0$ a positive slope angle. We adopt here the usual convention in geophysical flows, that considers positive angles $\theta$ for negative slopes. Let $(X, Z)$ be the local coordinates, measured in the direction along and normal to the inclined plane $\tilde{b}$, respectively, and $\mathbf{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).
<figure><img src="image_2.png" /><figcaption>Figure 1: Sketch of the Cartesian (blue) and local (red) reference system.</figcaption></figure>

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 U the characteristic height, length and velocity, we define the usual shallowness parameter ε = H/L, and denoting with tildes (·) the non-dimensional variables, we obtain

$$ (X, Z, t) = (L\tilde{X}, H\tilde{Z}, (L/U)\tilde{t}), $$

$$ (u, w) = (U\tilde{u}, \varepsilon U\tilde{w}), $$

$$ h = H\tilde{h}, \quad \rho = \rho_0\tilde{\rho}, \quad p = \rho_0 U^2 \tilde{p}, $$

$$ (\tau_{XX}, \tau_{XZ}, \tau_{ZZ}) = \rho_0 U^2 (\varepsilon \tilde{\tau}_{XX}, \tilde{\tau}_{XZ}, \varepsilon \tilde{\tau}_{ZZ}), $$

with

$$ \tilde{\tau}_{XX} = \tilde{\eta} \partial_{\tilde{X}} \tilde{u}, \quad \tilde{\tau}_{XZ} = \frac{\tilde{\eta}}{2} (\partial_{\tilde{Z}} \tilde{u} + \varepsilon^2 \partial_{\tilde{X}} \tilde{w}), \quad \tilde{\tau}_{ZZ} = \tilde{\eta} \partial_{\tilde{Z}} \tilde{w}, $$

and therefore

$$ \tilde{D}(\tilde{u}) = \frac{U}{H} \frac{1}{2} \begin{pmatrix} 2\varepsilon \partial_{\tilde{X}} \tilde{u} & \partial_{\tilde{Z}} \tilde{u} + \varepsilon^2 \partial_{\tilde{X}} \tilde{w} \\ \partial_{\tilde{Z}} \tilde{u} + \varepsilon^2 \partial_{\tilde{X}} \tilde{w} & 2\varepsilon \partial_{\tilde{Z}} \tilde{w} \end{pmatrix}. $$

Note also that $||\tilde{D}(\tilde{u})|| = |\partial_{\tilde{Z}}\tilde{u}|/2 + \mathcal{O}(\varepsilon)$. Finally, the Froude number is defined as $Fr = U/\sqrt{gH \cos\theta}$.

Following the asymptotic expansion performed in [42] for the pressure, and in particular for the non-hydrostatic counterpart, we write

$$ \tilde{p} = \frac{\rho_0}{Fr^2} (\tilde{b} + \tilde{h} - \tilde{Z}) + \varepsilon \tilde{p}_1^{nh} + \varepsilon^2 \tilde{p}^{nh}, \quad (2) $$

where the pressure is divided into a hydrostatic contribution, and the first and second order non-hydrostatic contributions $\tilde{p}_1^{nh}, \tilde{p}^{nh}$. 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)

$$ \partial_X u + \partial_Z w = 0, \quad (3a) $$

$$ \rho_0 (\partial_t u + u \partial_X u + w \partial_Z u) + \partial_X p = \frac{1}{\varepsilon} \frac{\rho_0}{Fr^2} \tan \theta + \varepsilon \partial_X \tau_{XX} + \frac{1}{\varepsilon} \partial_Z \tau_{XZ}, \quad (3b) $$

$$ \rho_0 \varepsilon^2 (\partial_t w + u \partial_X w + w \partial_Z w) + \varepsilon \partial_Z p_1^{nh} + \varepsilon^2 \partial_Z p^{nh} = \varepsilon \partial_X \tau_{XZ} + \varepsilon \partial_Z \tau_{ZZ}. \quad (3c) $$
This is the system that will be layer-averaged using the multilayer approach. Before, boundary conditions are written in non-dimensional form as

$$ \partial_t h + u_{|_{b+h}} \partial_X (b+h) - w_{|_{b+h}} = 0, \quad p_{|_{b+h}}^{nh} = 0, $$

at the free surface, and

$$ u_{|_b} \partial_X b - w_{|_b} = 0, \quad \left( \frac{\eta}{2} \partial_Z u \right)_{|_b} = \mu(I_{|_b}) p_{|_b} \frac{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 $\mu(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><img src="image_3.png" /><figcaption>Figure 2: Sketch of the multilayer discretization.</figcaption></figure>

We first consider $(l_\alpha)_{1\le\alpha\le N}$ strictly positive coefficients verifying $\sum_{\alpha=1}^N l_\alpha = 1$. Then, for a positive height $h(t, X)$ we consider $N$ shallow layers whose heights are $h_\alpha(t, X) = l_\alpha h(t, X)$ (see Figure 2). Notice that then $\sum_{\alpha=1}^N h_\alpha = h$ holds. Considering now the levels $Z_{\alpha+1/2} = b + \sum_{\beta=1}^\alpha h_\beta$, we denote the layers by

$$ \Omega_\alpha(t, X) = \{(t, X, Z) : Z_{\alpha-1/2} \le Z \le Z_{\alpha+1/2}\}, \quad \text{for } \alpha = 1, \dots, N. $$

Note that

$$ b = Z_{1/2} < \dots < Z_{\alpha+1/2} < \dots < Z_{N+1/2} = b+h $$

is a partition of $[b, b+h]$. As notation, we define

$$ f_{\alpha} = \frac{1}{h_{\alpha}} \int_{Z_{\alpha-1/2}}^{Z_{\alpha+1/2}} f(t, X, Z) \, dZ, $$

the averaged value in the layer $\Omega_\alpha$ of the function $f$. In order to deal with possible discontinuities of $f$ across the interface defined by $Z = Z_{\alpha+1/2}$, we define

$$ f_{\alpha+1/2}^{-} = \lim_{\substack{Z \to Z_{\alpha+1/2} \\ Z < Z_{\alpha+1/2}}} f_{|_{\Omega_{\alpha}(t)}}^{+}, \quad f_{\alpha+1/2}^{+} = \lim_{\substack{Z \to Z_{\alpha+1/2} \\ Z > Z_{\alpha+1/2}}} f_{|_{\Omega_{\alpha+1}(t)}}^{-} $$
and [[f]]$_{α+1/2}$ = $f_{α+1/2}^+$ − $f_{α+1/2}^-$ denotes the jump of f across 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 LDNH$_{0}$ in [40]. In this model, the two components of the velocity u and w are approximated by constant functions in each layer:

$$
u(t, X, Z) = \sum_{\alpha=1}^{N} u_{\alpha}(t, X) \mathbb{1}_{\Omega_{\alpha}}(Z), \quad w(t, X, Z) = \sum_{\alpha=1}^{N} w_{\alpha}(t, X) \mathbb{1}_{\Omega_{\alpha}}(Z),
$$

and the non-hydrostatic pressure $p^{nh}(t, X, Z)$, $p_{|\Omega_\alpha}^{nh} \in \mathbb{P}_1$ by a continuous linear function. Note that in this case:

$$
p_{\alpha}^{nh}(t, X) = p^{nh}(t, X, Z_{\alpha}) = \frac{p_{\alpha+1/2}^{nh}(t, X) + p_{\alpha-1/2}^{nh}(t, X)}{2} \quad (4)
$$

where $Z_α = Z_{α-1/2} + h_α/2$ is the midpoint of the layer $\Omega_α$. Actually, for each layer we have

$$
p^{nh}(Z) = p_{\alpha}^{nh} + \frac{p_{\alpha+1/2}^{nh} - p_{\alpha-1/2}^{nh}}{h_{\alpha}} (Z - Z_{\alpha}), \quad \text{for } Z \in [Z_{\alpha-1/2}, Z_{\alpha+1/2}].
$$

Under these assumptions and by appropriately defining the jump conditions across the interfaces, the LDNH$_{0}$ model for the Euler equation is obtained in [40]. Here, the main novelty is to deal with the viscous term $\partial_Z \tau_{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_{\alpha+1/2}$, denoted by $G_{\alpha+1/2} := G_{\alpha+1/2}^+ = G_{\alpha+1/2}^-$ (corresponding to $-\Gamma_{\alpha+1/2}$ in [40]), where

$$
G_{\alpha+1/2}^{\pm} = \partial_t Z_{\alpha+1/2} + u_{\alpha+1/2}^{\pm} \partial_X Z_{\alpha+1/2} - w_{\alpha+1/2}^{\pm}.
$$

Morover, from that condition we recover

$$
[[u]]_{\alpha+1/2} \cdot \vec{n}_{\alpha+1/2} = 0,
$$

where $\vec{n}_{\alpha+1/2} = (\partial_X Z_{\alpha+1/2}, -1)$ is the downward normal vector at the interface.

From the jump condition for the momentum equation, it follows that

$$
\epsilon^2 \tau_{XX, \alpha+1/2}^{\pm} \partial_X Z_{\alpha+1/2} - \tau_{XZ, \alpha+1/2}^{\pm} = \epsilon^2 \tilde{\tau}_{XX, \alpha+1/2} \partial_X Z_{\alpha+1/2} - \tilde{\tau}_{XZ, \alpha+1/2} \pm \frac{1}{2} \rho_0 \epsilon G_{\alpha+1/2} (u_{\alpha+1/2}^{+} - u_{\alpha+1/2}^{-}),
$$

$$
\tau_{ZX,\alpha+1/2}^{\pm} \partial_X Z_{\alpha+1/2} - \tau_{ZZ,\alpha+1/2}^{\pm} = \tilde{\tau}_{ZX,\alpha+1/2} \partial_X Z_{\alpha+1/2} - \tilde{\tau}_{ZZ,\alpha+1/2} \pm \frac{1}{2} \rho_0 \epsilon G_{\alpha+1/2} (w_{\alpha+1/2}^{+} - w_{\alpha+1/2}^{-}),
$$

where we recall that $\tau_{XX} = \eta \partial_X u$, $\tau_{XZ} = \frac{1}{2} \eta (\partial_Z u + \epsilon^2 \partial_X w)$, $\tau_{ZZ} = \eta \partial_Z w$, and $\tilde{\tau}$ is a consistent approximation of $\tau$ at the $Z = Z_{\alpha+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.
### 3.2.1 Mass conservation equation

The mass conservation equation (3a) is integrated along the normal direction in each layer, leading to

$$0 = \int_{Z_{\alpha-1/2}}^{Z_{\alpha+1/2}} (\partial_X u + \partial_Z w) dZ = \partial_X (h_\alpha u_\alpha) - (u_{\alpha+1/2}^{-} \partial_X Z_{\alpha+1/2} - w_{\alpha+1/2}^{-}) + (u_{\alpha-1/2}^{+} \partial_X Z_{\alpha-1/2} - w_{\alpha-1/2}^{+}),$$

and therefore

$$l_{\alpha} (\partial_t h + \partial_X (hu_{\alpha})) = G_{\alpha+1/2} - G_{\alpha-1/2}, \quad \text{for } \alpha = 1, \dots, N.$$

Notice that combining previous equations, the explicit formula for $G_{\alpha+1/2}$ is found (see [2])

$$G_{\alpha+1/2} = \sum_{\beta=1}^{\alpha} l_{\beta} \partial_X (h(u_{\beta} - \bar{u})), \quad \text{with} \quad \bar{u} = \sum_{\gamma=1}^{N} l_{\gamma} u_{\gamma}. \qquad (5)$$

In particular, by summing up all the equations we have the mass conservation equation

$$\partial_t h + \partial_X (h \bar{u}) = 0, \qquad (6)$$

where $G_{1/2} = G_{N+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

$$\int_{Z_{\alpha-1/2}}^{Z_{\alpha+1/2}} (\partial_t u + u \partial_X u + w \partial_Z u) dZ = \partial_t (h_\alpha u_\alpha) + \partial_X (h_\alpha u_\alpha^2) - u_{\alpha+1/2}^{-} G_{\alpha+1/2} + u_{\alpha-1/2}^{+} G_{\alpha-1/2}.$$

Note that in the present approach (LDNH₀ model) it holds that $u_{\alpha+1/2}^{-} = u_{\alpha-1/2}^{+} = u_{\alpha}$, since $u_{|\Omega_{\alpha}} \in \mathbb{P}_0$. Pressure terms are now integrated within the layer in the normal direction, leading to

$$\int_{Z_{\alpha-1/2}}^{Z_{\alpha+1/2}} \partial_X \left( \frac{\rho_0}{Fr^2} (\tilde{b} + b + h - z) + \varepsilon p_1^{nh} + \varepsilon^2 p^{nh} \right) dz = h_\alpha \frac{\rho_0}{Fr^2} \partial_X (\tilde{b} + b + h) \\ + \partial_X \left( \varepsilon h_\alpha p_{1,\alpha}^{nh} + \varepsilon^2 h_\alpha p_{\alpha}^{nh} \right) - \left( \varepsilon p_1^{nh} + \varepsilon^2 p^{nh} \right)_{|_{Z_{\alpha+1/2}}} \partial_X Z_{\alpha+1/2} + \left( \varepsilon p_1^{nh} + \varepsilon^2 p^{nh} \right)_{|_{Z_{\alpha-1/2}}} \partial_X Z_{\alpha-1/2},$$

where we have used $\tilde{b}(X) = (X_{end} - X) \sin \theta$. Finally, viscous terms are integrated, which yields

$$\int_{Z_{\alpha-1/2}}^{Z_{\alpha+1/2}} \left( \varepsilon \partial_X \tau_{XX} + \frac{1}{\varepsilon} \partial_Z \tau_{XZ} \right) dz \\ = \varepsilon \partial_X \left( \int_{Z_{\alpha-1/2}}^{Z_{\alpha+1/2}} \eta \partial_X u \, dz \right) + \frac{1}{\varepsilon} \left( \varepsilon^2 \tau_{XX, \alpha-1/2}^{+} \partial_X Z_{\alpha-1/2} - \tau_{XZ, \alpha-1/2}^{+} \right) - \frac{1}{\varepsilon} \left( \varepsilon^2 \tau_{XX, \alpha+1/2}^{-} - \tau_{XZ, \alpha+1/2}^{-} \right) \\ = \varepsilon \partial_X \left( \int_{Z_{\alpha-1/2}}^{Z_{\alpha+1/2}} \eta \partial_X u \, dz \right) + \left( \varepsilon \tilde{\tau}_{XX, \alpha-1/2} \partial_X Z_{\alpha+1/2} - \frac{1}{\varepsilon} \tilde{\tau}_{XZ, \alpha-1/2} \right) - \left( \varepsilon \tilde{\tau}_{XX, \alpha+1/2} \partial_X Z_{\alpha+1/2} - \frac{1}{\varepsilon} \tilde{\tau}_{XZ, \alpha+1/2} \right) \\ + \frac{1}{2} \rho_0 G_{\alpha-1/2} (u_{\alpha-1/2}^{+} - u_{\alpha-1/2}^{-}) + \frac{1}{2} \rho_0 G_{\alpha+1/2} (u_{\alpha+1/2}^{+} - u_{\alpha+1/2}^{-}).$$

Defining now

$$\tilde{K}_{\alpha+1/2} = \varepsilon \tilde{\tau}_{XX, \alpha+1/2} \partial_X Z_{\alpha+1/2} - \frac{1}{\varepsilon} \tilde{\tau}_{XZ, \alpha+1/2} \quad \text{and} \quad U_{\mathcal{Z}}^{H}_{\alpha+1/2} = \frac{u_{\alpha+1} - u_{\alpha}}{h_{\alpha+1/2}},$$

we obtain that

$$\tilde{K}_{\alpha+1/2} = K_{\alpha+1/2} + \mathcal{O}(\varepsilon), \quad \text{where} \quad K_{\alpha+1/2} = -\frac{1}{\varepsilon} \frac{\eta_{\alpha+1/2}}{2} U_{\mathcal{Z}}^{H}_{\alpha+1/2}.$$
We can collect finally all the terms, keep terms up to first order and terms of order $\epsilon^2$, go to the dimensional variables, and obtain the horizontal momentum equation

$$ \rho \left( \partial_t (h_\alpha u_\alpha) + \partial_X (h_\alpha u_\alpha^2) + g \cos \theta h_\alpha \partial_X (\tilde{b} + b + h) \right) + \partial_X (h_\alpha p_\alpha^{nh}) = K_{\alpha-1/2} - K_{\alpha+1/2} \\ + p_{\alpha+1/2}^{nh} \partial_X Z_{\alpha+1/2} - p_{\alpha-1/2}^{nh} \partial_X Z_{\alpha-1/2} + u_{\alpha+1/2} \rho G_{\alpha+1/2} - u_{\alpha-1/2} \rho G_{\alpha-1/2}, \quad (7) $$

for $\alpha = 1, \dots, N$ and

$$ u_{\alpha+1/2} = \frac{u_{\alpha} + u_{\alpha+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

$$ \int_{Z_{\alpha-1/2}}^{Z_{\alpha+1/2}} \varepsilon (\partial_X \tau_{ZX} + \partial_Z \tau_{ZZ}) dz \\ = \varepsilon \partial_X \left( \int_{Z_{\alpha-1/2}}^{Z_{\alpha+1/2}} \eta (\partial_Z u + \varepsilon^2 \partial_X w) dz \right) + \varepsilon \left( \tau_{ZX,\alpha-1/2}^+ \partial_X Z_{\alpha-1/2} - \tau_{ZZ,\alpha-1/2}^+ \right) - \varepsilon \left( \tau_{ZX,\alpha+1/2}^- \partial_X Z_{\alpha+1/2} - \tau_{ZZ,\alpha+1/2}^- \right) \\ = \varepsilon \partial_X \left( \int_{Z_{\alpha-1/2}}^{Z_{\alpha+1/2}} \eta (\partial_Z u + \varepsilon^2 \partial_X w) dz \right) + \varepsilon \left( \tilde{\tau}_{ZX,\alpha-1/2} \partial_X Z_{\alpha-1/2} - \tilde{\tau}_{ZZ,\alpha-1/2} \right) - \varepsilon \left( \tilde{\tau}_{ZX,\alpha+1/2} \partial_X Z_{\alpha+1/2} - \tilde{\tau}_{ZZ,\alpha+1/2} \right) \\ + \frac{1}{2} \rho_0 \varepsilon^2 G_{\alpha-1/2} (w_{\alpha-1/2}^+ - w_{\alpha-1/2}^-) + \frac{1}{2} \rho_0 \varepsilon^2 G_{\alpha+1/2} (w_{\alpha+1/2}^+ - w_{\alpha+1/2}^-). $$

It is observed that

$$ \varepsilon \partial_X \left( \int_{Z_{\alpha-1/2}}^{Z_{\alpha+1/2}} \eta (\partial_Z u + \varepsilon^2 \partial_X w) dz \right) + \varepsilon \left( \tilde{\tau}_{ZX,\alpha-1/2} \partial_X Z_{\alpha-1/2} - \tilde{\tau}_{ZZ,\alpha-1/2} \right) - \varepsilon \left( \tilde{\tau}_{ZX,\alpha+1/2} \partial_X Z_{\alpha+1/2} - \tilde{\tau}_{ZZ,\alpha+1/2} \right) $$

is defined only by terms of order $O(\varepsilon)$ and $O(\varepsilon^3)$.

Then, collecting all the terms and keeping also terms of order $\varepsilon^2$, the vertical momentum equation in dimensional variables is

$$ \rho (\partial_t (h_\alpha w_\alpha) + \partial_X (h_\alpha u_\alpha w_\alpha)) = p_{\alpha-1/2}^{nh} - p_{\alpha+1/2}^{nh} + w_{\alpha+1/2} \rho G_{\alpha+1/2} - w_{\alpha-1/2} \rho G_{\alpha-1/2}, \quad (8) $$

for $\alpha = 1, \dots, N$, where $w_{\alpha+1/2} = (w_{\alpha+1/2}^+ + w_{\alpha+1/2}^-)/2$. Notice that here first order term in $\varepsilon$ 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_\alpha, w_\alpha, p_\alpha^{nh}$) and $2N+1$ equations till now. Then, the incompressibility condition is used in each layer. Concretely, we integrate the incompressibility equation (3a) for $Z \in [Z_{\alpha-1/2}, Z_\alpha]$ and, taking into
account equations (6), (7) and (8), the final weakly non-hydrostatic multilayer system reads

$$ \partial_t h + \partial_X (h \bar{u}) = 0, \qquad (9a) $$

$$ -K_{\alpha+1/2} + p_{\alpha+1/2}^{nh} \partial_X Z_{\alpha+1/2} - p_{\alpha-1/2}^{nh} \partial_X Z_{\alpha-1/2} + u_{\alpha+1/2} \rho G_{\alpha+1/2} - u_{\alpha-1/2} \rho G_{\alpha-1/2}, \qquad (9b) $$

$$ \rho (\partial_t (h_\alpha w_\alpha) + \partial_X (h_\alpha u_\alpha w_\alpha)) = p_{\alpha-1/2}^{nh} - p_{\alpha+1/2}^{nh} + w_{\alpha+1/2} \rho G_{\alpha+1/2} - w_{\alpha-1/2} \rho G_{\alpha-1/2}, \qquad (9c) $$

$$ w_{\alpha} - u_{\alpha} \partial_X Z_{\alpha} + \sum_{\beta=1}^{\alpha-1} \partial_X (l_{\beta} h u_{\beta}) + \frac{\partial_X (l_{\alpha} h u_{\alpha})}{2} = 0, \quad \text{for } \alpha = 1, \dots, N \qquad (9d) $$

where $G_{\alpha+1/2}$ and $\bar{u}$ are given by (5), $\tilde{b}(X) = (X_{end} - X) \sin \theta$ is the reference plane, and the velocities at the interfaces are

$$ u_{\alpha+1/2} = \frac{u_{\alpha+1} + u_{\alpha}}{2}, \qquad w_{\alpha+1/2} = \frac{w_{\alpha+1} + w_{\alpha}}{2}, $$

since $u|_{\Omega_\alpha}, w|_{\Omega_\alpha} \in \mathbb{P}_0$. The viscous terms are

$$ K_{\alpha+1/2} = -\frac{1}{2} \eta_{\alpha+1/2} \mathcal{U}_{\mathcal{Z} \alpha+1/2}^{H}, \quad \text{with} \quad \mathcal{U}_{\mathcal{Z} \alpha+1/2}^{H} = \frac{u_{\alpha+1} - u_{\alpha}}{h_{\alpha+1/2}} $$

and $h_{\alpha+1/2} = (h_{\alpha+1} + h_{\alpha})/2$, for $\alpha = 1, \dots, N-1$. In the previous equation, the viscosity coefficients are

$$ \eta_{\alpha+1/2} = \frac{\mu(I_{\alpha+1/2})p_{\alpha+1/2}}{\sqrt{(\mathcal{U}_{\mathcal{Z} \alpha+1/2}^{H})^2 / 4 + \delta^2}}, \qquad (10) $$

with

$$ p_{\alpha+1/2} = \rho g \cos \theta h \sum_{\beta=\alpha+1}^{N} l_{\beta} + p_{\alpha+1/2}^{nh} \quad \text{and} \quad I_{\alpha+1/2} = \frac{d_s |\mathcal{U}_{\mathcal{Z} \alpha+1/2}^{H}|}{\sqrt{p_{\alpha+1/2}/\rho_s}}. $$

At the free surface we have $K_{N+1/2} = 0$, and at the bottom $K_{1/2}$ is defined by the Coulomb friction law

$$ K_{1/2} = -\mu(I_{1/2})p_{1/2} \frac{u_1}{|u_1|}. \qquad (11) $$

It can be easily proved that the previous model satisfies a dissipative energy balance as a consequence of the energy balance associated to model LDNH$_0$ [40]. Actually, the following result is obtained.

**Theorem 1** Let $E_\alpha$ be the energy of layer $\Omega_\alpha$ that is defined by

$$ E_{\alpha} = h_{\alpha} \left( \frac{u_{\alpha}^{2} + w_{\alpha}^{2}}{2} + g \cos \theta \left( \tilde{b} + b + \frac{h}{2} \right) \right). $$

System (9) satisfies the dissipative energy inequality

$$ \rho \partial_t \left( \sum_{\alpha=1}^{N} E_{\alpha} \right) + \rho \partial_X \left[ \sum_{\alpha=1}^{N} u_{\alpha} \left( E_{\alpha} + g \cos \theta h_{\alpha} \frac{h}{2} + \frac{p_{\alpha}^{nh}}{\rho} \right) \right] \leq -\rho g \cos \theta h |u_1| \mu(I) - \sum_{\alpha=1}^{N-1} \frac{\eta_{\alpha+1/2}}{2} \frac{(u_{\alpha+1} - u_{\alpha})^2}{h_{\alpha+1/2}}. $$
The proof of the theorem can be deduced following [40] and [36]. For the sake of completeness we give the details in Appendix B.

Concerning the steady states of system (9), we are interested in preserving steady states at rest where the slope of the free surface is lower than the angle of repose of the granular material. Concretely, the solutions verifying

$$u_{\alpha} = w_{\alpha} = p_{\alpha}^{nh} = 0, \quad \text{for } \alpha = 1, \dots, N, \quad \text{and} \quad \left| \partial_{x} (h + b + \tilde{b}) \right| \le \mu_{s}. \qquad (12)$$

Notice that these solutions are analogous to the steady states preserved by previous models: hydrostatic multilayer model [36] and also its one-layer version, non-hydrostatic one-layer model [42], and even the classical Savage-Hutter model [79].

**Remark 1** From system (9) one can recover previous models in [36, 42]. In particular, the non-hydrostatic monolayer model [42] is obtained when fixing the number of layers $N = 1$ in (9), and the hydrostatic model [36] is recovered using equations (9a), (9b) and vanishing non-hydrostatic pressure $p_{\alpha}^{nh} = 0$.

# 4 Numerical approximation

In the literature, multilayer hydrostatic systems have been discretized by combining finite volume methods with a splitting procedure (see, for instance, [3, 39]) where viscous terms are treated (semi-)implicitly. On the other hand, some other numerical methods have been proposed for non-hydrostatic pressure models. For instance, in [42] a one-layer non-hydrostatic pressure system for granular flows is discretized. In [33], a multilayer non-hydrostatic pressure without viscous effects is discretized. In the latter cases, the presence of non-hydrostatic terms on the model demands solving an elliptic problem implicitly. To our knowledge, it is the first time that a numerical scheme for a multilayer system with non-constant viscosity and non-hydrostatic pressure has been exposed, and we will describe it in the following.

We perform a three-step numerical approximation following [42]. Firstly, we deal with the hyperbolic counterpart (18) in an explicit way, including also some non-hydrostatic contributions. It is worth mentioning that these non-hydrostatic terms result from considering the pressure written in an incremental form (see (15)), as explained later. This step is based on an explicit path-conservative finite volume method, with a hydrostatic reconstruction, taking care of the well-balanced property. Secondly, the viscous effects are implicitly discretized by solving system (23). In this step, the viscous terms are linearized in order to not solve an implicit nonlinear problem. Finally, system (26) is implicitly solved to compute the non-hydrostatic pressure unknown, and then we appropriately correct the velocity and pressure fields.

The computational cost of this multilayer non-hydrostatic pressure system notably increases with the number of vertical layers, mainly because of the elliptic problem associated with the non-hydrostatic pressure. Then, designing some strategies to make these methods more efficient when a fine vertical resolution is needed is convenient. In this sense, we propose a new and very efficient procedure to solve the elliptic problem. It will be analyzed in the numerical tests (see Subsection 5.1). In what follows, we shall describe each step in detail.

System (9) can be written in compact form as

$$\left\{
\begin{array}{l}
\partial_t U + \partial_x F(U) + B(U)\partial_x U + S(U)\partial_x (\tilde{b} + b + h) = E(U, P) + T(U, P, \partial_x U, \partial_x P, \partial_x Z), \quad (13a) \\
C(U, \partial_x U, \partial_x Z) = 0,
\end{array}
\right.
$$

where $U = (h, hu_1, \dots, hu_N, hw_1, \dots, hw_N)' \in \mathbb{R}^{2N+1}$ and $P = (p_{1/2}^{nh}, \dots, p_{N+1/2}^{nh})' \in \mathbb{R}^{N+1}$ denote the conserved variables and non-hydrostatic unknowns (at the interfaces), and

$$F(U) = \begin{pmatrix} h\bar{u} \\ hu_1^2 \\ \vdots \\ hu_N^2 \\ hu_1w_1 \\ \vdots \\ hu_Nw_N \end{pmatrix}, \quad S(U) = \begin{pmatrix} 0 \\ g \cos \theta h \\ \vdots \\ g \cos \theta h \\ 0 \\ \vdots \\ 0 \end{pmatrix}, \quad E(U, P) = \frac{1}{\rho} \begin{pmatrix} 0 \\ (K_{1-1/2} - K_{1+1/2})/l_1 \\ \vdots \\ (K_{N-1/2} - K_{N+1/2})/l_N \\ 0 \\ \vdots \\ 0 \end{pmatrix}$$
are the convective flux, hydrostatic pressure, and viscous terms, respectively. $B(U)\partial_X U$ contains the non-conservative products involving the momentum transfer across the interfaces

$$B(U)\partial_X U = \begin{pmatrix} 0 \\ (u_{1+1/2}G_{1+1/2} - u_{1-1/2}G_{1-1/2})/l_1 \\ \vdots \\ (u_{N+1/2}G_{N+1/2} - u_{N-1/2}G_{N-1/2})/l_N \\ (w_{1+1/2}G_{1+1/2} - w_{1-1/2}G_{1-1/2})/l_1 \\ \vdots \\ (w_{N+1/2}G_{N+1/2} - w_{N-1/2}G_{N-1/2})/l_N \end{pmatrix}$$

Concerning non-hydrostatic terms, the vector

$$\mathcal{T} = \frac{1}{\rho} \begin{pmatrix} 0 \\ \left(-\partial_X(hp_1^{nh}) + p_{1+1/2}^{nh}\partial_X Z_{1+1/2} - p_{1-1/2}^{nh}\partial_X Z_{1-1/2}\right)/l_1 \\ \vdots \\ \left(-\partial_X(hp_N^{nh}) + p_{N+1/2}^{nh}\partial_X Z_{N+1/2} - p_{N-1/2}^{nh}\partial_X Z_{N-1/2}\right)/l_N \\ \left(p_{1-1/2}^{nh} - p_{1+1/2}^{nh}\right)/l_1 \\ \vdots \\ \left(p_{N-1/2}^{nh} - p_{N+1/2}^{nh}\right)/l_N \end{pmatrix},$$

contains the non-hydrostatic corrections for each momentum equation, where $Z = (Z_{1/2}, \dots, Z_{N+1/2})' \in \mathbb{R}^{N+1}$ is the vector containing the internal interfaces. Finally, $C(U, \partial_X U, \partial_X Z)$ describes the incompressibility constraints (9d), which will be rewritten and better described in Step 3 of the description of this numerical algorithm.

We consider a set of $N_x$ control volumes $\Omega_i = [x_{i-1/2}, x_{i+1/2}]$, with constant length $\Delta x$ as usual in the framework of finite volume methods. We denote by $x_i = (x_{i-1/2} + x_{i+1/2})/2$ the center of each volume cell. For any time $t^n \in [0, t_f]$ we consider the cell averages

$$U_i^n = \frac{1}{\Delta x} \int_{\Omega_i} U(x, t^n) dx.$$

For convenience, the non-hydrostatic unknowns are located at the half-integer locations $x_{i+1/2}$, $i = 0, \dots, N_x$. Then, we denote the non-hydrostatic pressure unknowns at time $t^n$ by

$$P^n = (P_{1/2}^n, \dots, P_{N_x+1/2}^n)' \in \mathbb{R}^{(N_x+1) \cdot (N+1)},$$

where

$$P_{i+1/2}^n = (p_{1/2}^{nh}(x_{i+1/2}, t^n), \dots, p_{N+1/2}^{nh}(x_{i+1/2}, t^n))' \in \mathbb{R}^{N+1},$$

that corresponds with a second-order approximation of the cell averages of the non-hydrostatic pressure on the staggered grid $[x_i, x_{i+1}]$.

Concerning the time discretization, the interval $[0, t_f]$ is discretized in time steps $\Delta t^n = t^{n+1} - t^n$, which are computed satisfying the usual stability condition

$$\Delta t^n = \min_{i=1,\dots,N_x} \Delta t_i^n, \quad \text{with} \quad \Delta t_i^n = \frac{CFL \Delta x}{\max_{i=1,\dots,N_x} (\lambda_i^n)}, \qquad (14)$$

where $CFL \in (0, 1)$, and $\lambda_i^n = |\bar{u}_i^n| + \sqrt{g \cos \theta h_i^n}$ is an upper bound of the system's eigenvalues.
Before describing each step of the scheme, let us mention that we consider an incremental decomposition of the non-hydrostatic unknowns to increase the order of accuracy of the projection method for the elliptic problem (see [50]). Thus, we define the fluctuation of the non-hydrostatic pressure from time $t^n$ to $t^{n+1}$ as follows:

$$ \tilde{P}_{i+1/2}^{n+1} = P_{i+1/2}^{n+1} - P_{i+1/2}^{n}, \qquad (15) $$

that is $\tilde{P}^{n+1} = P^{n+1} - P^n$. Note that this approach makes non-hydrostatic contribution $P_{i+1/2}^{n+1}$, at time $t^{n+1}$, to be written as sum of two terms:

$$ P^{n+1} = P^n + \tilde{P}^{n+1}, \qquad (16) $$

and therefore, we can split the non-hydrostatic contribution at the horizontal momentum equation as follows

$$ \mathcal{T}(U, P^{n+1}, \partial_X U, \partial_X P^{n+1}, \partial_X Z) = \mathcal{T}(U, P^n, \partial_X U, \partial_X P^n, \partial_X Z) + \mathcal{T}(U, \tilde{P}^{n+1}, \partial_X U, \partial_X \tilde{P}^{n+1}, \partial_X Z). \qquad (17) $$

We will discretize the first terms explicitly in time at step 1 of the scheme, whereas the second terms will be considered implicitly in step 3.

We then propose a numerical scheme based on three steps described in what follows. Let us previously remark that, as usually in splitting procedures, the outputs of first and second steps are denoted by $U^{n+1/3}$ and $U^{n+2/3}$ respectively. The results of step 3 is $U^{n+1}$, $P^{n+1}$.

### Step 1: hydrostatic problem

The first step focuses on solving explicitly the underlying system obtained when viscous and non-hydrostatic terms related to $\tilde{P}^{n+1}$ are removed from system (13):

$$ \partial_t U + \partial_X F_C(U) + B(U)\partial_X U + S(U)\partial_X (\tilde{b} + b + h) = \mathcal{T}(U, P, \partial_X U, \partial_X P, \partial_X Z), \qquad (18) $$

where $\mathcal{T}$ will be evaluated explicitly at time $t^n$ as suggested in (17).

Here, we consider the extension to the multilayer case of the scheme presented in [42] for the one-layer case, where the main difference is the presence of non-conservative terms coming from the mass transference between the vertical layers. Thus, the finite volume method is described as

$$ U_i^{n+1/3} = U_i^n + \frac{\Delta t}{\Delta x} \left( \mathcal{F}_{i-1/2}^n - \mathcal{F}_{i+1/2}^n + \frac{1}{2} \left( \mathcal{B}_{i-1/2}^n + \mathcal{B}_{i+1/2}^n + S_{i-1/2}^n + S_{i+1/2}^n \right) \right) + \mathcal{T} (U_i^n, P_i^n, (\partial_X U_i^n), (\partial_X P_i^n), (\partial_X Z_i^n)), \qquad (19) $$

with

$$ \mathcal{B}_{i+1/2} = \frac{1}{2} (\mathcal{B}(U_i) + \mathcal{B}(U_{i+1})) (U_{i+1} - U_i), $$

and

$$ S_{i+1/2} = \frac{1}{2} (S(U_i) + S(U_{i+1})) (h_{i+1/2}^{+} - h_{i+1/2}^{-}), $$

where $h_{i+1/2}^{\pm}$ are the reconstructed states

$$ \begin{aligned} h_{i+1/2}^{-} &= \max(0, h_i - (\Delta b_{i+1/2})_{+}), & \text{with } (\Delta b_{i+1/2})_{+} &= \max(0, b_{i+1} - b_i). \qquad (20) \\ h_{i+1/2}^{+} &= \max(0, h_{i+1} - (-\Delta b_{i+1/2})_{+}), & & \end{aligned} $$

The numerical flux $\mathcal{F}_{i+1/2}^n$ is defined as

$$ \mathcal{F}_{i+1/2}^n = \frac{1}{2} (\mathcal{F}_C(U_{i+1}^n) + \mathcal{F}_C(U_i^n)) - \frac{1}{2} \mathcal{D}_{i+1/2}^n, $$
being $D_{i+1/2}^n$ the numerical diffusion of the scheme, which will be defined within the framework of Polynomial Viscosity Methods (PVM) introduced in [22] for non-conservative hyperbolic systems. Concretely, we use the HLL method for non-conservative systems, which reads

$$ \mathcal{D}_{i+1/2} = \alpha_0 (\tilde{U}_{i+1/2}^+ - \tilde{U}_{i+1/2}^-) + \alpha_1 (\mathbf{F}_C(\mathbf{U}_{i+1}^n) - \mathbf{F}_C(\mathbf{U}_i^n) + \mathcal{B}_{i+1/2}^n + \tilde{\mathcal{S}}_{i+1/2}^n), \quad (21) $$

with

$$ \tilde{\mathcal{S}}_{i+1/2}^n = \frac{1}{2} (\mathbf{S}(\mathbf{U}_{i+1}^n) + \mathbf{S}(\mathbf{U}_i^n)) (\tilde{h}_{i+1/2}^{+,n} - \tilde{h}_{i+1/2}^{-,n}), $$

and $\tilde{U}_{i+1/2}^{\pm}$ the reconstructed states defined in (20). These reconstructed states are essential to satisfy the well-balanced property for the set of steady solutions defined by (12). In (21), the parameter $\alpha_i$ are

$$ \alpha_0 = \frac{S_R |S_L| - S_L |S_R|}{S_R - S_L}, \quad \alpha_1 = \frac{|S_R| - |S_L|}{S_R - S_L}, $$

defined in terms of the minimum ($S_L$) and maximum ($S_R$) wave speeds. These values are

$$ \begin{aligned} S_L &= \min \left( \bar{u}_i - \sqrt{g \cos \theta h_i}, \bar{u}_{i+1/2} - \sqrt{g \cos \theta h_{i+1/2}} \right), \\ S_R &= \max \left( \bar{u}_{i+1} + \sqrt{g \cos \theta h_{i+1}}, \bar{u}_{i+1/2} + \sqrt{g \cos \theta h_{i+1/2}} \right), \end{aligned} $$

with $h_{i+1/2}$, $\bar{u}_{i+1/2}$ the usual Roe's averaged states.

The main difficulty in this first step consists of preserving steady states defined by (12), where not only classical lake-at-rest solutions ($u_\alpha = 0$ and $\partial_X (\tilde{b} + b + h) = 0$) are involved, but also others with non-constant free surface. The latter case corresponds with those steady states where friction is greater than pressure forces. In that case, the usual hydrostatic reconstruction must be modified, accounting for the friction term, in order to preserve those steady states. To this aim, following [9] (see also [37, 42]), we define $\tilde{h}_{i+1/2}^{\pm}$ as in (20) by replacing $(\Delta b_{i+1/2})_+$ by

$$ (\tilde{\Delta} b_{i+1/2})_+ = \max (0, b_{i+1} - b_i + \Delta C_{i+1/2}), $$

accounting for the Coulomb friction. Concretely, the term $\Delta C_{i+1/2}$ is responsible for balancing friction and pressure forces and vanishing the numerical diffusion when it is necessary to obtain $\partial_t h = 0$. This term is defined as $\Delta C_{i+1/2} = -f_{i+1/2} \Delta x/(g \cos \theta)$ with

$$ f_{i+1/2} = - \operatorname{proj}_{g \cos \theta \mu(I_{1/2})} \left( \frac{\bar{u}_{i+1/2}}{\Delta t} - \frac{g \cos \theta (h_{i+1} + b_{i+1} - h_i - b_i)}{\Delta x} \right), $$

being

$$ \operatorname{proj}_{g \cos \theta \mu(I_{1/2})} (X) = \begin{cases} X, & \text{if } |X| \le g \cos \theta \mu(I_{1/2}), \\ g \cos \theta \mu(I_{1/2}) \frac{X}{|X|}, & \text{if } |X| > g \cos \theta \mu(I_{1/2}). \end{cases} $$

Finally, the reconstructed states $\tilde{U}_{i+1/2}^{\pm}$ in (21) are

$$ \begin{aligned} \tilde{U}_{i+1/2}^{-} &= (\tilde{h}_{i+1/2}^{-}, \tilde{h}_{i+1/2}^{-} u_{1,i}, \dots, \tilde{h}_{i+1/2}^{-} u_{N,i}, \tilde{h}_{i+1/2}^{-} w_{1,i}, \dots, \tilde{h}_{i+1/2}^{-} w_{N,i})', \\ \tilde{U}_{i+1/2}^{+} &= (\tilde{h}_{i+1/2}^{+}, \tilde{h}_{i+1/2}^{+} u_{1,i+1}, \dots, \tilde{h}_{i+1/2}^{+} u_{N,i+1}, \tilde{h}_{i+1/2}^{+} w_{1,i+1}, \dots, \tilde{h}_{i+1/2}^{+} w_{N,i+1})'. \end{aligned} $$

Once we have discretized the hydrostatic system, the explicit contributions of the non-hydrostatic pressure $\mathcal{T}$ in (19) are added. To this aim, let us recall that $p_\alpha^{nh}$ is given by (4). Moreover we write $\partial_X (hp_\alpha^{nh})$ in terms of $h$, $p_\alpha^{nh}$ and its derivatives by using the product rule, and

$$ P_i = \frac{P_{i+1/2} + P_{i-1/2}}{2}, \quad (\partial_X P)_i = \frac{P_{i+1/2} - P_{i-1/2}}{\Delta x}, \quad (\partial_X U)_i = \frac{U_{i+1} - U_{i-1}}{2\Delta x}. \quad (22) $$
as second-order point-values approximations of the averaged values. Finally, note that $Z_i$ can be written as

$$Z_i = (b_i, b_i + l_1 h_i, b_i + (l_1 + l_2)h_i, \dots, b_i + h_i)'$$

Note that the previous vector components describe the vertical interfaces $\Gamma_{\alpha+1/2}$, $\alpha = 0, \dots, N$, and the derivative of a given component $Z_{\alpha-1/2}$ is discretized as follows:

$$(\partial_x Z_{\alpha-1/2})_i = \frac{b_{i+1} - b_{i-1}}{2\Delta x} + \frac{h_{i+1} - h_{i-1}}{2\Delta x} \sum_{\beta=1}^{\alpha-1} l_{\beta}.$$

## Step 2: Viscous terms

In this step the viscous terms are accounted, by solving the system

$$\partial_t U = E(U, P). \qquad (23)$$

First, the heights are no longer modified along the scheme due to the form of vectors $E$ and $T$. Then, we trivially have $h_i^{n+1} = h_i^{n+2/3} = h_i^{n+1/3}$. To introduce the viscous terms, we consider an implicit method, where the viscous terms are linearized by freezing the viscosity coefficients at time $t^n$, as in [37], in order to avoid the resolution of a non-linear system. It leads to a $N \times N$ tridiagonal linear system at each cell $\Omega_i$. Concretely, the new values at this second step for the horizontal discharges (the vertical ones are not involved in $E(U, P)$) are given by

$$U_i^{n+2/3} = U_i^{n+1/3} + \Delta t E(U_i^n, U_i^{n+1/3}, U_i^{n+2/3}, P_i^n).$$

The previous equation is written by components, as

$$\hat{h}_i^{n+2/3} u_{\alpha,i}^{n+2/3} = (hu)_i^{n+1/3} + \frac{\Delta t}{2\rho l_\alpha} \left( \eta_{\alpha+\frac{1}{2}} (U_i^n, P_i^n) \frac{u_{\alpha+1,i}^{n+2/3} - u_{\alpha,i}^{n+2/3}}{l_{\alpha+\frac{1}{2}} h_i^{n+1}} - \eta_{\alpha-\frac{1}{2}} (U_i^n, P_i^n) \frac{u_{\alpha,i}^{n+2/3} - u_{\alpha-1,i}^{n+2/3}}{l_{\alpha-\frac{1}{2}} h_i^{n+1}} \right), \quad (24)$$

with $\eta_{\alpha+1/2}$ given by (10), and

$$\hat{h}_i^{n+2/3} = \begin{cases} h_i^{n+2/3} + \frac{1}{2}(h_{i-1}^{n+2/3} + h_{i+1}^{n+2/3}) & \text{if } |h_{i+1} + b_{i+1} + \tilde{b}_{i+1} - (h_{i-1} + b_{i-1} + \tilde{b}_{i-1})| \le 2\Delta x \mu_s \\ h_i^{n+2/3} & \text{otherwise.} \end{cases} \quad (25)$$

a reconstruction of $h^{n+2/3}$ in the cell $\Omega_i$. This definition, that is similar to the one used in [42], is essential to the well-balanced property of the scheme, as it is shown in the numerical tests (see Subsection 5.3). Notice that first condition in (25) holds in the steady state case $|\partial_x (h + b + \tilde{b})| \le \mu_s$. Let us remark that equation (24) must be accordingly modified for according with the boundary conditions $K_{N+1/2} = 0$, and $K_{1/2}$ in (11). Eq. (24) leads to $N_x$ tridiagonal linear systems, whose dimension is $N \times N$, which are solved using the Thomas' algorithm. Notice that nonlinear systems would be obtained in case of considering a fully implicit scheme in this step.

**Remark 2** In these steps, two numerical aspects should be noticed. First, in order to be able to deal with wet/dry areas: (i) in Step 1, a standard wet/dry treatment (see [21]) is used to avoid the appearance of negative height values in case of emerging bottom situations. It consists of detecting the emerging bottom (for instance $h_i > 0$, $h_{i+1} = 0$ and $h_i + b_i < b_{i+1}$) and using the modified state $\tilde{U}_i = (h_i, 0, \dots, 0)'$ to compute the numerical fluxes (the symmetric case is analogous); (ii) in Step 2, when the total height is lower than a given tolerance, then the viscous terms are discretized as a Coulomb friction term for a one-layer model (see [37] for details). In practice, we set this tolerance as $10^{-6}$. Thus, a one-layer system with Coulomb friction is solved in those cells satisfying $h_i < 10^{-6}$ for $i = 1, \dots, N_x$, instead of the linear system (24), which removes the vertical (multilayer) structure in these cells. Second, we cannot obtain $u_\alpha = 0$ because of the regularization method, but minimal velocities. In practice, we reach velocities $u_\alpha \sim 10^{-6}, 10^{-7}$ for $\delta = 10^{-5}$ (see Subsection 5.3).
### 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:

$$ \left\{ \begin{array}{l} \partial_t U = \mathcal{T}(U, P, \partial_X U, \partial_X P, \partial_X Z), \\ \mathcal{C}(U, \partial_X U, \partial_X Z) = 0, \end{array} \right. \qquad (26) $$

where $\mathcal{T}$ is evaluated at time $t^{n+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 $\alpha$ by

$$ \mathcal{E}_{\alpha} = 0, \quad \forall \alpha = 1, \dots, N, \text{ where } \mathcal{E}_{\alpha} = w_{\alpha} - u_{\alpha} \partial_X Z_{\alpha} + \sum_{\beta=1}^{\alpha-1} \partial_X (l_{\beta} h u_{\beta}) + \frac{\partial_X (l_{\alpha} h u_{\alpha})}{2}. $$

As it can be seen, the imposition of the constraint for $\alpha = 2, \dots, N$ involves variables from the layer $\beta = 1$ to $\beta = \alpha - 1$. To simplify the set of constraints (9d), and to improve the diagonal dominance of the linear operator to be inverted, we equivalently consider

$$ \mathcal{N}_{\alpha} = 0, \quad \forall \alpha = 1, \dots, N, \text{ where } \mathcal{N}_1 = \mathcal{E}_1, \quad \mathcal{N}_2 = \mathcal{E}_2 - \mathcal{E}_1, \quad \dots, \quad \mathcal{N}_{\alpha} = \mathcal{E}_N - \mathcal{E}_{N-1}. $$

For clarity, we will assume in what follows that $l_{\alpha} = 1/N$. Then, we define the vector $\mathcal{C}$ in (13b) by

$$ \mathcal{C}(U, \partial_X U, \partial_X Z) = 2hN(\mathcal{N}_1, \dots, \mathcal{N}_N)', \qquad (27) $$

that reads

$$ \mathcal{C} = \begin{pmatrix} 2Nhw_1 - 2Nhu_1\partial_X Z_1 + h\partial_X(hu_1) & & \\ & 2N(hw_2 - hw_1) - 2N(hu_2\partial_X Z_2 - hu_1\partial_X Z_1) + h\partial_X hu_2 + h\partial_X hu_1 & \\ & & \vdots \\ & & 2N(hw_N - hw_{N-1}) - 2N(hu_N\partial_X Z_N - hu_{N-1}\partial_X Z_{N-1}) + h\partial_X hu_N + h\partial_X hu_{N-1} \end{pmatrix}, \qquad (28) $$

where we recall that $Z_{\alpha} = Z_{\alpha+1/2} - \frac{1}{2N}h$. We consider now a projection method, taking into account the fluctuation $\tilde{P}^{n+1}$ defined in (15). To do that, we will impose the incompressibility constraints on the staggered mesh $x_{i+1/2}$

$$ \mathcal{C}(\tilde{U}_{i+1/2}^{n+1}, (\partial_X \tilde{U}^{n+1})_{i+1/2}, (\partial_X Z)_{i+1/2}) = 0, \qquad (29) $$

where

$$ \tilde{U}_{i+1/2}^{n+1} = \tilde{U}_{i+1/2}^{n+2/3} + \Delta t \mathcal{T} \left( \tilde{U}_{i+1/2}^{n+1}, \tilde{P}_{i+1/2}^{n+1}, (\partial_X \tilde{U}^{n+1})_{i+1/2}, (\partial_X \tilde{P}^{n+1})_{i+1/2}, (\partial_X Z)_{i+1/2} \right), \qquad (30) $$

and the conserved variables and their derivatives on the staggered mesh are approximated from the averaged values

$$ \tilde{U}_{i+1/2}^{n+2/3} = \frac{\tilde{U}_i^{n+2/3} + \tilde{U}_{i+1}^{n+2/3}}{2}, \quad \tilde{U}_{i+1/2}^{n+1} = \frac{\tilde{U}_i^{n+1} + \tilde{U}_{i+1}^{n+1}}{2}, \quad (\partial_X \tilde{U}^{n+1})_{i+1/2} = \frac{\tilde{U}_{i+1}^{n+1} - \tilde{U}_i^{n+1}}{\Delta x}. $$

Similarly, the component $(\partial_X Z_{\alpha-1/2})$ of $(\partial_X Z)$ is given by

$$ (\partial_X Z_{\alpha-1/2})_{i+1/2} = \frac{h_{i+1/2}^{+} - h_{i+1/2}^{-} - (h_{i+1} - h_i)}{\Delta x} + \frac{\alpha - 1}{N} \cdot \frac{h_{i+1} - h_i}{\Delta x}, $$

where the first addend is an approximation of $\partial_x b$ using the hydrostatic reconstruction (20). Finally, the derivatives of the fluctuation of the non-hydrostatic pressures are approximated by

$$ (\partial_X \tilde{P}^{n+1})_{i+1/2} = \frac{\tilde{P}_{i+3/2}^{n+1} - \tilde{P}_{i-1/2}^{n+1}}{2\Delta x}. $$
**Remark 3** Notice that the factor h in the definition of C in (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 h must be performed when computing C.

Now, putting Eqs. (30) into (29) yields a linear system

$$ \mathcal{A} \cdot \mathcal{P} = \frac{1}{\Delta t} \mathcal{R}, \qquad (31) $$

where P is a vector that contains a rearrangement of the non-hydrostatic pressure fluctuations unknowns

$$ \mathcal{P} = \begin{pmatrix} \mathcal{P}_1 \\ \vdots \\ \mathcal{P}_N \end{pmatrix}, \quad \mathcal{P}_\alpha = \begin{pmatrix} \tilde{p}_{\alpha,1/2}^{nh,n+1} \\ \vdots \\ \tilde{p}_{\alpha,N_x+1/2}^{nh,n+1} \end{pmatrix}, \quad \alpha = 1, \dots, N, $$

A is a block tridiagonal matrix of dimension $N \cdot (N_x + 1) \times N \cdot (N_x + 1)$ given by

$$ \mathcal{A} = \begin{pmatrix} M_{1,1} & M_{1,2} & & & \\ M_{2,1} & M_{2,2} & M_{2,3} & & \\ & & M_{3,2} & \ddots & \\ & & & \ddots & M_{N-1,N} \\ & & & & M_{N,N-1} \\ & & & & M_{N,N} \end{pmatrix}, \qquad (32) $$

being $M_{i,j}$ tridiagonal matrices of dimension $(N_x+1) \times (N_x+1)$ (see details in Appendix A). Finally, $\mathcal{R}$ is the right-hand side of the linear system that is given by

$$ \mathcal{R} = \begin{pmatrix} \mathcal{R}_1 \\ \vdots \\ \mathcal{R}_N \end{pmatrix}, \quad \mathcal{R}_\alpha = \begin{pmatrix} \mathcal{C}_\alpha (U_{1/2}^{n+2/3}, (\partial_X U^{n+2/3})_{1/2}, (\partial_X Z)_{1/2}) \\ \vdots \\ \mathcal{C}_\alpha (U_{N_x+1/2}^{n+2/3}, (\partial_X U^{n+2/3})_{N_x+1/2}, (\partial_X Z)_{N_x+1/2}) \end{pmatrix}, \quad \alpha = 1, \dots, N, \qquad (33) $$

where the sub-index $\alpha$ in $\mathcal{C}_\alpha$ denotes the $\alpha$-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: $\mathcal{P}^{(0)} \leftarrow I.C.$

2: $s \leftarrow 1$

3: $Er \leftarrow 0$

4: **do**

5: **for** $\alpha = 1, 2, \dots, N$ **do**

6: Solve the tridiagonal system: $M_{\alpha,\alpha} \mathcal{P}_\alpha^{(s)} = \frac{1}{\Delta t} \mathcal{R}_\alpha - M_{\alpha,\alpha-1} \mathcal{P}_{\alpha-1}^{(s-1)} - M_{\alpha,\alpha+1} \mathcal{P}_{\alpha+1}^{(s-1)}$

7: $Er \leftarrow Er + \Delta x \|\mathcal{P}_\alpha^{(s-1)} - \mathcal{P}_\alpha^{(s)}\|^p$

8: $\mathcal{P}_\alpha^{(s-1)} \leftarrow \mathcal{P}_\alpha^{(s)}$

9: **end for**

10: **while** $Er \ge \epsilon_{tol}^p$

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 $\tilde{\mathcal{P}}^n = \mathcal{P}^n - \mathcal{P}^{n-1}$, although other choices are possible. As stated in line 6 of
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 $l_p$ norms for $p = 1$, or $p = 2$, used as the stopping criteria and given by

$$ ||\mathcal{P}_{\alpha}^{(s-1)} - \mathcal{P}_{\alpha}^{(s)}||^p = \frac{1}{N} \sum_{i=0}^{N_x} h_{i+1/2} \left| \tilde{p}_{\alpha,i+1/2}^{nh,(s-1)} - \tilde{p}_{\alpha,i+1/2}^{nh,(s)} \right|^p, \quad \text{for } p = 1, 2. $$

Alternative, the $l_\infty$ norm can be considered, and therefore $Er$ in line 7 is updated as

$$ Er \leftarrow \max \left\{ Er, \max_{i=0, \dots, N_x} \left\{ \left| \tilde{p}_{\alpha, i+1/2}^{nh,(s-1)} - \tilde{p}_{\alpha, i+1/2}^{nh,(s)} \right| \right\} \right\}^i $$

whereas the stopping condition in line 10 is given by $Er \ge \epsilon_{tol}^p$. The influence of the choice of that different norms will be studied in Subsection 5.1.1. Finally, $\epsilon_{tol}^p$ 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 5 layers, 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 $\tilde{P}_{i+1/2}^{n+1}$ have been computed, the discharges $U_i^{n+1}$ can be updated from

$$ U_i^{n+1} = U_i^{n+2/3} + \Delta t \mathcal{T} \left( U_i^{n+1}, \tilde{P}_i^{n+1}, (\partial_x U_i^{n+1})_i, (\partial_x \tilde{P}_i^{n+1})_i, (\partial_x Z)_i \right), $$

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 $P_{i+1/2}^{n+1}$ by (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^\infty$-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.
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 $N_{ini guess} < N$, so that its solution is a better initial guess for the original system with N layers. In this way, we would solve the dispersive linear system twice: firstly with $N_{ini 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 $P^{n+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 $N_{gr} = N/N_{ini guess} \in \mathbb{N}$ representing the number of layers that merge into a single layer. Then, the domain is divided along the normal direction into layers $\Omega_{\beta,red}$, where

$$ \Omega_{\beta,red} = \bigcup_{j=1}^{N_{gr}} \Omega_{N_{gr}(\beta-1)+j} \quad \text{for } \beta = 1, \dots, N_{ini guess}. $$

Now, we need to define auxiliary vector variables with size $N_{ini 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_{\beta,red} = \frac{1}{N_{gr}} \sum_{j=1}^{N_{gr}} u_{N_{gr}(\beta-1)+j} \quad \text{for } \beta = 1, \dots, N_{ini guess}, $$

and analogously for all the variables involved in resolving the elliptic problem (30). Concretely, it is done to define $U_{red}$, $P_{red}$ and $\tilde{P}_{red}$, 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 $\tilde{P}_{red}$, it is used to generate the initial guess for the second stage, where the full linear system with N layers is solved. Thus, we define, for $\alpha = 1, \dots, N$,

$$ \tilde{p}_{\alpha+1/2}^{nh,(0)} = \tilde{p}_{\beta+1/2,red}^{nh} $$

with $1 \le \beta \le N_{ini guess}$ such that $Z_{\beta-1/2} \le Z_{\alpha-1/2} < Z_{\beta+1/2}$, and $\tilde{p}_{N+1/2}^{nh,(0)} = 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 $\mu(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 $\epsilon_{tol}^p = 10^{-6}$ as tolerance for the iterative solver finding the non-hydrostatic pressure (see line 10 of the algorithm in Section 4), and $\delta = 10^{-5}$ as regularization parameter in the viscosity coefficient (equation (10)). Note that this choice of $\delta$ is in agreement with other works for viscoplastic fluids, as [41, 25, 58].
<table><caption>Table 1: Values of the rheological parameters used in the numerical tests.</caption><thead><tr><th>μ<sub>s</sub></th><th>μ<sub>2</sub></th><th>d<sub>s</sub> (mm)</th><th>I<sub>0</sub></th><th>φ<sub>s</sub></th></tr></thead><tbody><tr><td>tan(25.5°) ≈ 0.48</td><td>tan(36°) ≈ 0.73</td><td>0.7</td><td>0.279</td><td>0.62</td></tr></tbody></table>

## 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 4 is based on an iterative solver to solve a linear system whose size is $N \times (N_x + 1)$, where N is the number of layers in the normal direction and $N_x$ the 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] m$ and 640 cells for the horizontal discretization, and the dam break problem given by the initial conditions

$$ \tilde{b}(X) = (3 - X) \sin \theta, \quad b(X) = 0, \quad h_0(X) = \begin{cases} 0.14 + h_i & \text{if } x \le 0, \\ h_i & \text{otherwise,} \end{cases} \quad (34) $$

with $h_i = 1.82 mm$ and $\theta = 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 $t_f = 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 \times (N_x + 1)$. Therefore, the computational time is sensitive to the norm used in the stopping criterion for this solver.

<table><caption>Table 2: Computational times $t_{comp}$ (s) and maximum number of iterations made in a single time step by the iterative solver for different number of normal layers $N$ and norms $l_2, l_\infty$.</caption><thead><tr><th>Layers</th><th>l<sub>2</sub> - N<sup>0</sup> Iter</th><th>l<sub>2</sub> - t<sub>comp</sub></th><th>l<sub>∞</sub> - N<sup>0</sup> Iter</th><th>l<sub>∞</sub> - t<sub>comp</sub></th></tr></thead><tbody><tr><td>20</td><td>431</td><td>779.7</td><td>693</td><td>1096.7</td></tr><tr><td>10</td><td>145</td><td>84.8</td><td>211</td><td>197.7</td></tr><tr><td>8</td><td>101</td><td>58.8</td><td>142</td><td>114.3</td></tr><tr><td>5</td><td>46</td><td>18.1</td><td>62</td><td>36.4</td></tr><tr><td>2</td><td>10</td><td>4.2</td><td>14</td><td>5.2</td></tr></tbody></table>

We have measured the computational time and the number of iterations made by the solver for different number of vertical layers. Table 2 shows 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 3 shows 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 $l_2$ norm the computational effort is lower than with the $l_\infty$ 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_\infty$ norm is always greater than with the $l_2$ norm (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 $l_1$ norm, the results are similar to the $l_2$ norm.

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 $l_2$-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 3 we 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,
<figure><img src="image_4.png" /><figcaption>Figure 3: (a): Number of iterations needed for $N = 5, 10, 20$ layers with the $l_2$ and $l_\infty$ norms. (b): Computational time in minutes for different number of layers with norms $l_2, l_\infty$.</figcaption></figure>

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.

<table><caption>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 4 in the non-hydrostatic case (NH). We also show this percentage for Step 2 in the hydrostatic case (H).</caption><thead><tr><th>Layers</th><th>1</th><th>2</th><th>5</th><th>10</th><th>20</th></tr></thead><tbody><tr><td>t<sub>comp</sub> Step 3/t<sub>comp</sub> (NH)</td><td>28%</td><td>58%</td><td>81%</td><td>92%</td><td>98%</td></tr><tr><td>t<sub>comp</sub> Step 2/t<sub>comp</sub> (NH)</td><td>5%</td><td>4%</td><td>2%</td><td>0.6%</td><td>0.2%</td></tr><tr><td>t<sub>comp</sub> Step 2/t<sub>comp</sub> (H)</td><td>9%</td><td>11%</td><td>11%</td><td>10%</td><td>9%</td></tr></tbody></table>

### 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.

<table><caption>Table 4: Computational times and speed-ups for non-hydrostatic simulations with (a) 8 (b) 10, (c) 20 vertical layers, reducing to $N_{ini guess}$ layers for the first step of the algorithm to solve the non-hydrostatic pressure.</caption><thead><tr><th colspan="3">(a) N = 8 layers</th><th colspan="3">(b) N = 10 layers</th></tr><tr><th>N<sub>ini guess</sub></th><th>-</th><th>4</th><th>-</th><th>5</th><th>2</th></tr></thead><tbody><tr><td>t<sub>comp</sub> (s)</td><td>47.3</td><td>44.6</td><td>77.46</td><td>84.8</td><td>128.2</td></tr><tr><td>Speed-up</td><td>-</td><td>1.06</td><td>0.61</td><td>-</td><td>0.66</td></tr></tbody></table><table><thead><tr><th colspan="3">(c) N = 20 layers</th></tr><tr><th>N<sub>ini guess</sub></th><th>-</th><th>10</th><th>5</th><th>4</th><th>2</th></tr></thead><tbody><tr><td>t<sub>comp</sub> (s)</td><td>594.2</td><td>183.2</td><td>166.8</td><td>215.2</td><td>564.6</td></tr><tr><td>Speed-up</td><td>-</td><td>3.24</td><td>3.56</td><td>2.76</td><td>1.05</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>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.</figcaption></figure>

In Table 4 we show the computational times when reducing to $N_{\text{ini 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 $N_{\text{ini guess}}$. In this case, when reducing to $N_{\text{ini 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 $N_{\text{ini guess}} = 5$ the speed-up is 3.56 (72% faster) and for $N_{\text{ini guess}} = 5$ the speed-up is 2.76 (64% faster). When using $N_{\text{ini 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 $N_{\text{ini 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 $N_{\text{ini disp}} = 5$.

In Figure 4 we study the case of 20 layers and $N_{\text{ini disp}} = 10$ in more detail. Figure 4a shows the computational times needed for different choices of $N_{\text{ini disp}}$, while in Figure 4b we see an analysis when $N_{\text{ini 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

$$ \tilde{b}(X) = (1.5 - X) \sin 20^{\circ}, \quad b(X) = 0, \quad h_0(X) = 0.02 + 0.4e^{-5(x-0.25)^2}. $$
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 N and $N_x$ at the same time. That is, we refine the mesh in both directions $(x, z)$. Here, $L^1$ absolute errors are computed, which are given by

$$ \|h - h^{ref}\|_1 = \Delta x \sum_{i=1}^{M} |h_i - \tilde{h}_i^{ref}|, \quad \|v - v^{ref}\|_1 = \Delta x \sum_{i=1}^{N_x} \frac{h_i}{N} \sum_{\alpha=1}^{N} |v_{\alpha,i} - \tilde{v}_{\alpha,i}^{ref}|, \quad (35) $$

where $h^{ref}$, $\tilde{v}^{ref}$ 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 $N_x = 6400$ and $N = 128$ will be the reference solution.

<table><caption>Table 5: $L^1$ errors and related orders for the height ($h$) and the velocity fields ($u(z)$, $w(z)$) for the solution at time $t = 0.15$ s.</caption><thead><tr><th>(N<sub>x</sub>, N)</th><th>Error h</th><th>Order h</th><th>Error u(z)</th><th>Order u(z)</th><th>Error w(z)</th><th>Order w(z)</th></tr></thead><tbody><tr><td>(100,2)</td><td>1.59×10<sup>-2</sup></td><td>-</td><td>4.14×10<sup>-2</sup></td><td>-</td><td>1.59×10<sup>-2</sup></td><td>-</td></tr><tr><td>(200,4)</td><td>1.22×10<sup>-2</sup></td><td>0.39</td><td>2.62×10<sup>-2</sup></td><td>0.66</td><td>1.22×10<sup>-2</sup></td><td>0.68</td></tr><tr><td>(400,8)</td><td>6.40×10<sup>-3</sup></td><td>0.93</td><td>1.16×10<sup>-2</sup></td><td>1.18</td><td>6.40×10<sup>-3</sup></td><td>1.29</td></tr><tr><td>(800,16)</td><td>3.06×10<sup>-3</sup></td><td>1.06</td><td>5.58×10<sup>-3</sup></td><td>1.05</td><td>3.06×10<sup>-3</sup></td><td>1.22</td></tr><tr><td>(1600,32)</td><td>1.35×10<sup>-3</sup></td><td>1.18</td><td>2.65×10<sup>-3</sup></td><td>1.07</td><td>1.35×10<sup>-3</sup></td><td>1.18</td></tr></tbody></table>

Table 5 shows the $L^1$ errors 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 ($N_x, 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]$ $m$ with 300 cells and free-outflow boundary conditions. The initial conditions (see figure 5) are $\tilde{b}(X) = 0$,

$$ b(X) = \frac{2}{X+4} - \frac{1}{2} + e^{-10(X+3/2)^2} + 0.15 \text{RAND}(X), \quad h_0(X) = 0.8 - \mu_s(X-2) - b(X), $$

being $\text{RAND}(X) \in [0, 1]$ a random value generated at run time, and the flow is assumed at rest $(hu)_\alpha = (hw)_\alpha = 0$ for $\alpha = 1, \dots, N$. We take the rheological parameter in Table 1 except for $\mu_s = \tan(20^\circ)$. $CFL = 0.7$ is used here and the final time simulation is $t_{fin} = 1$ s.

<figure><img src="image_6.png" /><figcaption>Figure 5: Bottom topography and initial free surface the well-balanced test.</figcaption></figure>

The $L^1$ errors for all the variables are measured, as in (35), where in this case the reference solution $\tilde{v}^{ref}$ is the initial state. Table 6 shows the maximum along the simulated time of these errors for different values of the regularization
parameter $\delta$. We show these errors for the scheme described in Section 4 with the definition of $\hat{h}_i^{n+2/3}$ given by (25) (Table 6a) and by considering $\hat{h}_i^{n+2/3} = h_i^{n+2/3}$ (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 $\delta$ 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 viscoplastic (Herschel-Bulkley) fluids.

<table><caption>Table 6: Maximum along the time of L¹ errors for the height (h) and the velocity fields (u(z), w(z)) with respect to the initial state by considering the definition of $\hat{h}_i^{n+2/3}$ given by (a) eq. (25) (well-balanced scheme); (b) $h_i^{n+2/3}$ (no well-balanced scheme).</caption><thead><tr><th colspan="4">(a) WB: ĥ<sub>i</sub><sup>n+2/3</sup> (25)</th><th colspan="4">(b) No WB: ĥ<sub>i</sub><sup>n+2/3</sup> = h<sub>i</sub><sup>n+2/3</sup></th></tr><tr><th></th><th>Error h</th><th>Error u(z)</th><th>Error w(z)</th><th></th><th>Error h</th><th>Error u(z)</th><th>Error w(z)</th></tr></thead><tbody><tr><td>δ = 10<sup>-5</sup></td><td>5.34×10<sup>-5</sup></td><td>5.38×10<sup>-5</sup></td><td>4.21×10<sup>-5</sup></td><td>δ = 10<sup>-5</sup></td><td>5.22×10<sup>-4</sup></td><td>1.20×10<sup>-3</sup></td><td>7.42×10<sup>-4</sup></td></tr><tr><td>δ = 10<sup>-10</sup></td><td>5.34×10<sup>-10</sup></td><td>5.39×10<sup>-10</sup></td><td>4.21×10<sup>-10</sup></td><td>δ = 10<sup>-10</sup></td><td>1.76×10<sup>-4</sup></td><td>1.93×10<sup>-4</sup></td><td>3.48×10<sup>-4</sup></td></tr><tr><td>δ = 10<sup>-15</sup></td><td>2.99×10<sup>-16</sup></td><td>5.39×10<sup>-15</sup></td><td>4.21×10<sup>-15</sup></td><td>δ = 10<sup>-15</sup></td><td>4.62×10<sup>-5</sup></td><td>6.06×10<sup>-5</sup></td><td>1.87×10<sup>-4</sup></td></tr></tbody></table>

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 $\delta$, we must use $\epsilon_{tol}^p < \delta$. Then, we use here $\epsilon_{tol}^p = \delta/100$. Notice that we recover the expected well-balanced property of the scheme even in the limit case $\partial_X (\tilde{b} + b + h) = \mu_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^{-5}$ for any value of $\delta = 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 $\delta$ in this test because it is a steady at rest solution. In general transient problem we use a large value ($\delta = 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 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 $\theta$ 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 $\tilde{b}$ and 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 $p_{\alpha \pm 1/2}^{nh} Z_{\alpha \pm 1/2}$ in the horizontal momentum equations (9b). Concretely, $Z_{\alpha \pm 1/2}$ is defined in terms of the local bottom bottom $b$ and not $\tilde{b}$ in 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 \in [-3, 3]$ m and 640 nodes. The material is initially at rest, and we set a final time $t_f = 10$ s and CFL = 0.5. In this case we use free-outflow conditions at both boundaries. Figure 6 shows the initial conditions for the bottom and the mass thickness in both coordinates system. In the Cartesian case, they are

$$b_{Cart}(x) = 1 - \tanh(x) + \kappa e^{-10(x-1)^2} + 0.5e^{-10(x-3)^2},$$
and

$$h_{Cart}(x) = \max(\eta(x) - b_{Cart}(x), 0) + h_i, \quad \text{with} \quad \eta(x) = \begin{cases} y_0 - 1 + e^{-0.5(x-x_0)^2}, & \text{if } x \le 0, \\ 0, & \text{otherwise,} \end{cases}$$

with $\kappa = 0.3$, $h_i = 10^{-4}$ m, $x_0 = -1.53$ m and $y_0 = 1.91$ m. In that (Cartesian) case, we fix $\theta = 0^\circ$ and $\tilde{b}_{Cart}(x) = 0$. Note also that a small thickness $h_i$ has 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 $\tilde{b}_{loc} = -0.7 - \tan\theta(x - 3)$, with $\theta = 25^\circ$, is considered. Now the local bottom $b_{loc}(X)$ is the distance from $\tilde{b}_{loc}(x)$ to the bottom measured in the direction normal to the plane $\tilde{b}_{loc}(x)$. Analogously, $h_{loc}(X)$ is measured as the distance from $b_{loc}(X)$ to the free surface.

<figure><img src="image_7.png" /><figcaption>Figure 6: Sketch of the initial conditions in both Cartesian and local coordinates. $\kappa = 0.3$. The variables defined in the Cartesian system are in blue and those in the Local system are in red.</figcaption></figure>

### 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 7 shows 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 $b_{Cart}$ is also represented and varies from $-15^\circ$ to $50^\circ$ (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 $\kappa = 0.45$ here. Figure 9 shows 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
<figure><img src="image_8.png" /><figcaption>Figure 7: Deposit (at $t = 10$ s) for the hydrostatic (H) and non-hydrostatic (NH) models in both Cartesian and local coordinates with $\kappa = 0.3$. Upper figure (a): General view and slope of the bottom $b_{Cart}(x)$ (positive angles correspond to negative slopes). Lower figures: Zooms ((b): zone located around the position of the initial mass; (c): zone located downslope before and after the bump).</figcaption></figure>

town would be safe, while the H-Cartesian model would indicate that this town would be largely affected. Notice that this is an academic test showing the importance of having a model as accurate as possible in order to act against natural hazards. In any case, it is expected that NH- models are more accurate than H- models since they incorporate this very relevant physical contribution. It can be seen in Subsection 5.5 when comparing with laboratory experiments for a simpler configuration. In addition, local models are more accurate than Cartesian models in these situations (strong topography gradients) since they compute the velocity in the downslope/normal direction to the reference plane, and not in the $(x, z)$-Cartesian directions. A more accurate approach would be to consider curvilinear coordinates where the velocity vector in measured in the downslope/normal direction to the topography at every location (see [14, 35]).

In the rest of the paper, we only consider local models, with both hydrostatic and weakly non-hydrostatic pressure.

## 5.5 Comparison with granular collapse experiments

We use here the proposed model to simulate laboratory experiments of granular column collapse performed in [61] and compare numerical and experimental results.

<table><caption>Table 7: Summary of notation of the models considered in this section, all of them in local coordinates.</caption><thead><tr><td>Notation</td><td>Model description</td></tr></thead><tbody><tr><td>NH-k l.</td><td>Non-hydrostatic pressure, k vertical layers</td></tr><tr><td>H-k l.</td><td>Hydrostatic pressure, k vertical layers</td></tr><tr><td>NH-μ<sub>w</sub>-k l.</td><td>Non-hydrostatic pressure, k vertical layers, sidewalls friction (36)</td></tr></tbody></table>

All the models considered in this section are summarized in Table 7. We remind that all the models are considered
<figure><img src="image_9.png" /><figcaption>Figure 8: (a) Differences between the hydrostatic local (H-Local) and the non-hydrostatic Cartesian (NH-Cartesian) models, and (b) influence of the coordinate system for hydrostatic (H) and non-hydrostatic (NH) models, measured in $l_2$-norm. Results for the model with 1 layer (squares), 5 (crosses), 10 (circles) and 20 (no symbols) layers.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Deposit at rest ($t = 10$ s) for the hydrostatic (H) and non-hydrostatic (NH) models in both Cartesian and local coordinates, and for a higher bump ($\kappa = 0.45$).</figcaption></figure>

in local coordinates in this subsection. The experimental setup consists of a granular column, whose thickness is $h_0 = 14$ cm, length 20 cm and width $W = 10$ cm. The flow is confined between two walls and initially at rest, on a plane inclined at a slope $\theta$, which in practice takes the values $\theta = 0°, 9.78°, 16°, 19°, 22°$. Once the gate is opened, the material flows until it forms a deposit.

We take the computational domain and initial conditions as in Subsection 5.1, that is $X \in [-0.2, 3]$ m with 640 horizontal cells and the initial height and bottom defined by (34), where $h_i$ is the thickness of a possible erodible bed made of the same material that could cover the plane. For the discretization in the direction normal to the reference plane, we consider 10 layers. We will see later that we have certain layer-independence of the results for large and small slopes with this number of vertical layers (see Table 9). The stability condition is set to $CFL = 0.1$. We need to use this small value in order to recover the correct behaviour in the influence of the erodible bed over the runout for small slopes, as it is explained later. Here a wall boundary condition is used upstream and free-outflow condition downstream. Note
that the topography is an inclined plane and there is no deviation from this reference frame.

### 5.5.1 Comparison of hydro/non-hydrostatic and mono/multilayer models

Figures 10, 11 and 12 show the evolution of the granular mass and the final deposits for all the slopes considered and for the rigid case ($h_i = 0$ mm). We show the results of four different models: the hydrostatic one-layer, hydrostatic multilayer, non-hydrostatic one-layer and non-hydrostatic multilayer models, denoted H-1 l., H-10 l., NH-1 l. and NH-10 l., respectively (see Table 7). We see that the non-hydrostatic approach allows to be closer to the experiment during the first instants of the collapse for all slopes, as expected. This is because at short times, when the material starts flowing, the non-hydrostatic effects are important. We see that for small values of the slope, there is no significant difference between the multilayer and the one-layer models, whereas these models differ more clearly for $\theta = 19^\circ, 22^\circ$. This suggests that in-depth variations solved with the multilayer model have more influence for large slopes while non-hydrostatic effects are more important at smaller slopes. Concerning the deposits, we see an excellent agreement of the NH-10l. model with the experiments for $\theta = 0^\circ, 9.78^\circ, 16^\circ$. For $\theta = 19^\circ$ the NH-1 l. model is the most accurate, and for $\theta = 22^\circ$ it is the H-10 l. model. In general, the final runout is too short with the NH-10 l. model for $\theta = 19^\circ, 22^\circ$.

We have also computed the errors (measured in $l_2$-norm) made by all the models, for all the slopes and for rigid ($h_i = 0$ mm) and erodible ($h_i > 0$) beds, when compared with the laboratory experiments. Figure 13 shows the results for the final deposits and also at a time close to the initial mass release. At this time, we see that the non-hydrostatic models produce smaller errors in all cases, and the non-hydrostatic multilayer model always has a smaller associated error. Concerning the deposits, the error depends on the slope, as discussed above. In general, for small slopes the NH-10 l. model gives smaller errors, but it is not the case for large slopes, mainly because the final run-out is underestimated with this model. This can be observed, for example, for high $h_i$ high and high slopes, where the 1 layer models are better since the simulated runouts are larger. Nevertheless, it is also possible that other effects are involved in the evolution of granular flows such as the release of the initial gate or wall effects, see for example [10] as discussed below.

In the experiments, it was observed that the runout distance $r_f$ increases when the thickness of the erodible bed $h_i$ increases for $\theta \ge 16^\circ$. In [36] it was shown that the combination of the multilayer approach and the variable friction coefficient $\mu(I)$ allows to recover this behaviour, contrary to one-layer models. Figure 14 shows that this behaviour is qualitatively reproduced by both the NH-10 l. and H-10 l multilayer models. On the contrary, Figure 11a shows that the runout distance decreases with increasing $h_i$ for simulations with the one-layer models H-1 l. and NH-1 l. for $\theta = 16^\circ$. We must remark that this behaviour of the multilayer models is not achieved for $\theta = 16^\circ$ if a value CFL =0.5 is used.

The main advantage of multilayer models is to recover the in-depth variations of the flow and in particular the flow/no flow transition as illustrated in Figure 15. Indeed, during the collapse an upper layer of material is flowing over a lower

Slope angle $\theta = 0^\circ$

<figure><img src="image_11.png" /><figcaption>Figure 10: Time evolution of the granular mass thickness over a plane of slope $\theta = 0^\circ$ for laboratory experiments (solid-circle blue line), and the hydrostatic (H) and non-hydrostatic (NH) models when considering 1 or 10 layers.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>(b) Slope angle $\theta = 16^{\circ}$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: Time evolution of the granular mass thickness over a plane of slope (a) $\theta = 9.78^{\circ}$ and (b) $16^{\circ}$, for laboratory experiments (solid-circle blue line), and the hydrostatic (H) and non-hydrostatic (NH) models when considering 1 or 10 layers.</figcaption></figure>

layer of static materials. Obviously, this behaviour cannot be reproduced by one-layer models. Figure 16 illustrates the ability of multilayer models to recover the different shapes of the normal velocity profiles from Bagnold-type to S-shape profiles, as also shown in [36]. The flow/no-flow interface has a similar shape with the hydrostatic and non-hydrostatic models at $\theta = 22^{\circ}$. They are more different at $\theta = 0^{\circ}$, again suggesting that non-hydrostatic effects have more influence for flows on smaller slopes. At certain times, we have also represented the results of [63] obtained by solving the full Navier-Stokes equations. We see both hydrostatic and non-hydrostatic multilayer models do not reproduce very well the flow/no-flow interface. Concretely, the thickness of static layer is underestimated. This may result from the fact that viscous contributions in the non-hydrostatic pressure are not accounted for.

Non-hydrostatic models much better reproduce the acceleration phase of the front than hydrostatic models as shown in Figure 17 for all the slopes and when considering a rigid bed ($h_i = 0$). From a qualitative analysis we see how the
<figure><img src="image_14.png" /><figcaption>(b) Slope angle $\theta = 22^\circ$</figcaption></figure>

<figure><img src="image_15.png" /></figure>

Figure 12: Time evolution of the granular mass thickness over a plane of slope (a) $\theta = 19^\circ$ and (b) $\theta = 22^\circ$, for laboratory experiments (solid-circle blue line), and the hydrostatic (H) and non-hydrostatic (NH) models when considering 1 or 10 layers.

front velocity. Indeed, with hydrostatic models (one-layer and multilayer), the initial velocity is much too high and the front has a maximal velocity at the initial instant and then the front decelerate until it stops, contrary to lab-experiments where a first acceleration phase is observed before the deceleration phase. These acceleration/deceleration phases are qualitatively reproduced with non-hydrostatic models. This is an important advantage of the NH-10 l. model compared to the H-10 l. model.

### 5.5.2 Sidewall effects

It is well-known that sidewall friction has a strong influence on the flow dynamics in narrow channels. To roughly account for these effects, the coefficients of the $\mu(I)$-rheology have been increased by 0.1 (Table 1), following [52]. We now investigate if adding non-hydrostatic effects makes it possible to obtain reasonable results when accounting more
appropriately sidewall friction as in [37], without having to increase the friction coefficients by 0.1. In practice, it consists of adding the term

$$M_{\alpha} = -l_{\alpha}h \frac{2}{W} \mu_{w} p_{\alpha} \frac{u_{\alpha}}{|u_{\alpha}|} \qquad (36)$$

to the right-hand side term in the momentum equation (9b), where W in the channel width, $\mu_w$ the friction coefficient at the lateral walls, and $p_\alpha$ the total pressure in the middle of the layer $\Omega_\alpha$. This term is discretized in Step 2 together with the other viscous terms as in [37]. Concretely, we add the contribution

$$-\frac{2\Delta t}{W\rho} \mu_w h_i^{n+1} p_{\alpha,i}^{n} \frac{u_{\alpha,i}^{n+2/3}}{\sqrt{|u_{\alpha,i}^n|^2 + \delta^2}}$$

<figure><img src="image_16.png" /><figcaption>Figure 13: Relative $l_2$-errors on the height along the domain computed with the hydrostatic (H) and non-hydrostatic (NH) models. Square-symbols lines correspond to the errors at final deposits, while cross-symbols lines are the errors at small times $t_{ini} = 0.24 (0^{\circ}), 0.32 (22^{\circ})$ s.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 14: Influence of the thickness $h_i$ of the erodible bed on the final runout $r_f$ for slopes $\theta = 16^{\circ}, 19^{\circ}, 22^{\circ}$ observed in the experiments (solid-circle blue line) and obtained with different simulations using hydrostatic (H) and non-hydrostatic (NH) models.</figcaption></figure>
<table><caption>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.</caption><thead><tr><th>W (cm)</th><th>model</th><th>μ<sub>s</sub></th><th>μ<sub>2</sub></th><th>μ<sub>w</sub></th></tr></thead><tbody><tr><td>10</td><td>NH-μ<sub>w</sub></td><td>tan(20.9°) ≈ 0.38</td><td>tan(32.76°) ≈ 0.64</td><td>tan(10.5°)</td></tr><tr><td>20</td><td>NH</td><td>tan(23.4°) ≈ 0.43</td><td>tan(34.74°) ≈ 0.69</td><td>0</td></tr><tr><td>20</td><td>NH-μ<sub>w</sub></td><td>tan(20.9°) ≈ 0.38</td><td>tan(32.76°) ≈ 0.64</td><td>tan(10.5°)</td></tr></tbody></table>

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 − μ<sub>w</sub>, and the rheological parameters used are summarized in Table 8 for 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(μ<sub>s</sub>) = 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-μ<sub>w</sub> 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><img src="image_18.png" /><figcaption>Figure 15: Flow/no-flow interface at different times at θ = 9.78° (left-hand side) and θ = 22° (right-hand side) and h<sub>i</sub> = 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].</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 16: Normal profiles of the downslope normalized velocity ($10u_{\alpha,i}/(\max_{\beta,j}|u_{\beta,j}|)$, $\beta = 1, \dots, N$, $j = 1, \dots, N_x$) obtained with the NH-20 l. model for $\theta = 22^\circ$ and $h_i = 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.</figcaption></figure>

<table><caption>Table 9: $l_2$-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 $t_{ini} = 0.24$ s ($\theta = 0°$), $0.32$ s ($\theta = 22°$), and at final time $t_{fin} = 1.5$ s ($\theta = 0°$), $3$ s ($\theta = 22°$). $t_{comp}$ is the computational time needed to reach $t_{fin}$.</caption><thead><tr><th colspan="3">(a) θ = 0°, h<sub>i</sub> = 1 mm</th><th colspan="3">(b) θ = 22°, h<sub>i</sub> = 1.82 mm</th></tr><tr><th>Model</th><th>Err(%)<br>(t<sub>ini</sub>)</th><th>Err(%)<br>(t<sub>fin</sub>)</th><th>t<sub>comp</sub><br>(s)</th><th>Model</th><th>Err(%)<br>(t<sub>ini</sub>)</th><th>Err(%)<br>(t<sub>fin</sub>)</th><th>t<sub>comp</sub><br>(s)</th></tr></thead><tbody><tr><td>NH-20 l.</td><td>-</td><td>-</td><td>1096.7 (19.08 min)</td><td>NH-20 l.</td><td>-</td><td>-</td><td>1517.4 (25.3 min)</td></tr><tr><td>NH-10 l.</td><td>0.30</td><td>0.33</td><td>368.2 (6.1 min)</td><td>NH-10 l.</td><td>0.92</td><td>1.92</td><td>544.6 (9.1 min)</td></tr><tr><td>NH-8 l.</td><td>0.45</td><td>0.47</td><td>197.7 (3.3 min)</td><td>NH-8 l.</td><td>1.31</td><td>2.87</td><td>290.6 (4.8 min)</td></tr><tr><td>NH-5 l.</td><td>0.79</td><td>0.88</td><td>60.1 (1.0 min)</td><td>NH-5 l.</td><td>2.19</td><td>5.58</td><td>89.2 (1.5 min)</td></tr><tr><td>NH-2 l.</td><td>1.68</td><td>1.93</td><td>14.2</td><td>NH-2 l.</td><td>4.70</td><td>14.18</td><td>21.0</td></tr><tr><td>NH-1 l.</td><td>3.05</td><td>2.30</td><td>6.4</td><td>NH-1 l.</td><td>6.47</td><td>22.80</td><td>9.7</td></tr><tr><td>H-20 l.</td><td>10.01</td><td>8.32</td><td>56.1</td><td>H-20 l.</td><td>7.36</td><td>11.21</td><td>79.4</td></tr><tr><td>H-10 l.</td><td>10.41</td><td>8.58</td><td>28.8</td><td>H-10 l.</td><td>8.23</td><td>12.90</td><td>41.0</td></tr><tr><td>H-1 l.</td><td>11.65</td><td>11.25</td><td>4.7</td><td>H-1 l.</td><td>11.78</td><td>31.93</td><td>7.1</td></tr></tbody></table>

### 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 9 we give the difference between the mass thickness
<figure><img src="image_20.png" /><figcaption>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 hᵢ = 0 mm. Here h₀ = 0.14 m, v₀ = √g cos θh₀ m/s and τc = √h₀/(g cos θ) s.*</figcaption></figure>

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 (hᵢ = 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 4 and 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
<figure><img src="image_21.png" /><figcaption>Figure 18: *Granular mass over a plane of slope $\theta = 0°$ ($h_i = 0$ mm), 16° ($h_i = 1.4$ mm), 22° ($h_i = 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 1 and 8). Note that at $\theta = 22°$, the granular mass is not at rest yet.*</figcaption></figure>

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) \times N$, being $M$ the number of cells and $N$ the 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.
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.

<figure><img src="image_22.png" /><figcaption>Figure 19: Time evolution of the normalized velocity of the front computed with non-hydrostatic models with and without the correction in the $\mu(I)$ coefficients accounting for the sidewall friction, and experimental data (solid-circle blue lines), for $W = 10, 20$ cm. Here $h_0 = 0.14$ m, $v_0 = \sqrt{g \cos \theta h_0}$ m/s and $\tau_c = \sqrt{h_0/(g \cos \theta)}$ s.</figcaption></figure>
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

$$
\mathcal{R} = \begin{pmatrix} \mathcal{R}_1 \\ \vdots \\ \mathcal{R}_N \end{pmatrix}, \quad \mathcal{R}_\alpha = \begin{pmatrix} 0 \\ 0 \\ \mathcal{C}_\alpha \left( \mathcal{U}_{5/2}^{n+2/3}, (\partial_X \mathcal{U}^{n+2/3})_{5/2}, (\partial_X \mathcal{Z})_{5/2} \right) \\ \vdots \\ \mathcal{C}_\alpha \left( \mathcal{U}_{N_x-5/2}^{n+2/3}, (\partial_X \mathcal{U}^{n+2/3})_{N_x-5/2}, (\partial_X \mathcal{Z})_{N_x-5/2} \right) \\ 0 \\ 0 \end{pmatrix}, \quad \alpha = 1, \dots, N.
$$

The tridiagonal matrices $M_{\alpha,\alpha-1}$, $M_{\alpha,\alpha}$, and $M_{\alpha,\alpha+1}$ for $\alpha = 2, \dots, N$ in (32) are described, for $i,j = 3, \dots, N_x - 2$, by

$$
\begin{align*}
(M_{\alpha,\alpha-1})_{i,j} &= \begin{cases}
\frac{h_{i+1/2}^2}{\Delta x^2} - h_{i+1/2} \frac{\xi_{\alpha,i+1/2} + \xi_{\alpha+1/2,i+1/2}}{\Delta x} + N^2, & \text{for } j=i-1, \\
-2 \frac{h_{i+1/2}^2}{\Delta x^2} + 4 \xi_{\alpha,i+1/2}^2 + 2h_{i+1/2} \partial_x \xi_{\alpha,i+1/2} + 2N^2, & \text{for } j=i, \\
\frac{h_{i+1/2}^2}{\Delta x^2} + h_{i+1/2} \frac{\xi_{\alpha,i+1/2} + \xi_{\alpha+1/2,i+1/2}}{\Delta x} + N^2, & \text{for } j=i+1,
\end{cases} \\
(M_{\alpha,\alpha})_{i,j} &= \begin{cases}
2 \frac{h_{i+1/2}^2}{\Delta x^2} - h_{i+1/2} \frac{\xi_{\alpha,i+1/2} + \xi_{\alpha+3/2,i+1/2} - \xi_{\alpha+1,i+1/2} - \xi_{\alpha-1/2,i+1/2}}{\Delta x} - 2N^2, & \text{for } j=i-1, \\
-4 \frac{h_{i+1/2}^2}{\Delta x^2} - 4 \left( \xi_{\alpha,i+1/2}^2 + \xi_{\alpha+1,i+1/2}^2 \right) - 2h_{i+1/2} \left( \partial_x \xi_{\alpha,i+1/2} + \partial_x \xi_{\alpha+1,i+1/2} \right) - 4N^2, & \text{for } j=i, \\
2 \frac{h_{i+1/2}^2}{\Delta x^2} + h_{i+1/2} \frac{\xi_{\alpha,i+1/2} + \xi_{\alpha+3/2,i+1/2} - \xi_{\alpha+1,i+1/2} - \xi_{\alpha-1/2,i+1/2}}{\Delta x} - 2N^2, & \text{for } j=i+1
\end{cases} \\
(M_{\alpha,\alpha+1})_{i,j} &= \begin{cases}
\frac{h_{i+1/2}^2}{\Delta x^2} + h_{i+1/2} \frac{\xi_{\alpha+1,i+1/2} + \xi_{\alpha+1/2,i+1/2}}{\Delta x} + N^2, & \text{for } j=i-1, \\
-2 \frac{h_{i+1/2}^2}{\Delta x^2} + 4 \xi_{\alpha+1,i+1/2}^2 - 2h_{i+1/2} \partial_x \xi_{\alpha+1,i+1/2} + 2N^2, & \text{for } j=i, \\
\frac{h_{i+1/2}^2}{\Delta x^2} - h_{i+1/2} \frac{\xi_{\alpha+1,i+1/2} + \xi_{\alpha+1/2,i+1/2}}{\Delta x} + N^2, & \text{for } j=i+1,
\end{cases}
\end{align*}
$$

with $M_{N,N+1} = 0$, and

$$
\xi_{\alpha,i+1/2} = (\alpha - 1/2) \frac{h_{i+1} - h_i}{\Delta x} + N \frac{h_{i+1/2}^{+} - h_{i+1/2}^{-} - (h_{i+1} - h_i)}{\Delta x},
$$

that is proportional to the derivative of the midpoint of the layer $Z_\alpha$, and we approximate their second order derivatives
as follows

$$
\partial_x \xi_{\alpha, i+1/2} = \frac{1}{\Delta x} \text{minmod} \left( \xi_{\alpha, i+3/2} - \xi_{\alpha, i+1/2}, \xi_{\alpha, i+1/2} - \xi_{\alpha, i-1/2}, \frac{1}{2} (\xi_{\alpha, i+3/2} - \xi_{\alpha, i-1/2}) \right),
$$

minmod being the usual slope limiter function.
The special cases of $M_{1,1}$, and $M_{1,2}$, corresponds to the discretization at the lowest layer, where the constraint (27) is simply given by $2hNN_1 = 2hNE_1$, and reads

$$ (M_{1,1})_{i,j} = \begin{cases} \frac{h_{i+1/2}^2}{\Delta x^2} - h_{i+1/2} \frac{\xi_{\alpha+3/2,i+1/2} - \xi_{\alpha+1,i+1/2}}{\Delta x} - N^2, & \text{for } j = i-1, \\ -2 \frac{h_{i+1/2}^2}{\Delta x^2} - 4\xi_{\alpha+1,i+1/2}^2 + 2h_{i+1/2}\partial_x\xi_{\alpha+1,i+1/2} - 2N^2, & \text{for } j=i, \\ \frac{h_{i+1/2}^2}{\Delta x^2} + h_{i+1/2} \frac{\xi_{\alpha+3/2,i+1/2} - \xi_{\alpha+1,i+1/2}}{\Delta x} - N^2, & \text{for } j=i+1 \end{cases} $$ $$ (M_{1,2})_{i,j} = \begin{cases} \frac{h_{i+1/2}^2}{\Delta x^2} + h_{i+1/2} \frac{\xi_{\alpha+1,i+1/2} + \xi_{\alpha+1/2,i+1/2}}{\Delta x} + N^2, & \text{for } j = i-1, \\ -2 \frac{h_{i+1/2}^2}{\Delta x^2} + 2N^2, & \text{for } j=i, \\ \frac{h_{i+1/2}^2}{\Delta x^2} - h_{i+1/2} \frac{\xi_{\alpha+1,i+1/2} + \xi_{\alpha+1/2,i+1/2}}{\Delta x} + N^2, & \text{for } j=i+1. \end{cases} $$

Regarding the horizontal wall/free-outflow boundary conditions along the independent variable $x$, we consider, for $\alpha = 1, 2, \dots, N$,

$$ \begin{aligned} (M_{\alpha,\alpha})_{1,1} &= (M_{\alpha,\alpha})_{2,2} = (M_{\alpha,\alpha})_{N_x-1,N_x-1} = (M_{\alpha,\alpha})_{N_x,N_x} = 1, \\ (M_{\alpha,\alpha})_{1,2} &= (M_{\alpha,\alpha})_{2,3} = (M_{\alpha,\alpha})_{N_x-1,N_x-2} = (M_{\alpha,\alpha})_{N_x,N_x-1} = -1, \end{aligned} \quad (37) $$

and for the rest of matrices, $M_{\alpha,\alpha \pm 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 1 concerning the energy balance for model (9). Some of the steps of this proof are similar to the proof of the energy inequality verified by LDNH$_0$ model 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

$$ \begin{aligned} & h_{\alpha} \left( \partial_t u_{\alpha} + \partial_x \left( \frac{u_{\alpha}^2}{2} \right) \right) + u_{\alpha} (G_{\alpha+1/2} - G_{\alpha-1/2}) + \partial_x \left( \frac{h_{\alpha} p_{\alpha}^{nh}}{\rho} \right) + g \cos \theta h_{\alpha} \partial_x (\tilde{b} + b + h) \\ & = \frac{1}{\rho} (K_{\alpha-1/2} - K_{\alpha+1/2}) + \frac{p_{\alpha+1/2}^{nh}}{\rho} \partial_x z_{\alpha+1/2} - \frac{p_{\alpha-1/2}^{nh}}{\rho} \partial_x z_{\alpha-1/2} + u_{\alpha+1/2} G_{\alpha+1/2} - u_{\alpha-1/2} G_{\alpha-1/2}. \end{aligned} \quad (38) $$

Now, we sum up (38) multiplied by $u_{\alpha}$ with the equation obtained by multiplying the mass conservation equation by $u_{\alpha}^2/2 + g \cos \theta(\tilde{b} + b + h)$, getting

$$ \begin{aligned} & \partial_t \left( h_{\alpha} \frac{u_{\alpha}^2}{2} \right) + g \cos \theta (\tilde{b} + b + h) \partial_t h_{\alpha} + \partial_x \left( h_{\alpha} u_{\alpha} \left( \frac{u_{\alpha}^2}{2} + g \cos \theta (\tilde{b} + b + h) \right) \right) \\ & - g \cos \theta (\tilde{b} + b + h) (G_{\alpha+1/2} - G_{\alpha-1/2}) + \frac{1}{\rho} u_{\alpha} \partial_x (h_{\alpha} p_{\alpha}^{nh}) - \frac{1}{\rho} u_{\alpha} \left( p_{\alpha+1/2}^{nh} \partial_x z_{\alpha+1/2} - p_{\alpha-1/2}^{nh} \partial_x z_{\alpha-1/2} \right) \\ & - G_{\alpha+1/2} \left( \frac{u_{\alpha}^2}{2} - u_{\alpha} u_{\alpha+1/2} \right) + G_{\alpha-1/2} \left( \frac{u_{\alpha}^2}{2} - u_{\alpha} u_{\alpha-1/2} \right) + \frac{1}{\rho} u_{\alpha} (K_{\alpha-1/2} - K_{\alpha+1/2}) \end{aligned} \quad (39) $$

Similar steps are followed for the equations of the vertical velocity (9c), and using the mass equation, we get

$$ \begin{aligned} & \partial_t \left( h_{\alpha} \frac{w_{\alpha}^2}{2} \right) + \partial_x \left( h_{\alpha} u_{\alpha} \frac{w_{\alpha}^2}{2} \right) + \frac{1}{\rho} w_{\alpha} \left( p_{\alpha+1/2}^{nh} - p_{\alpha-1/2}^{nh} \right) = \\ & \qquad - G_{\alpha+1/2} \left( \frac{w_{\alpha}^2}{2} - w_{\alpha} \tilde{w}_{\alpha+1/2} \right) + G_{\alpha-1/2} \left( \frac{w_{\alpha}^2}{2} - w_{\alpha} \tilde{w}_{\alpha-1/2} \right). \end{aligned} \quad (40) $$
By defining

$$E_{\alpha} = h_{\alpha} \left( \frac{u_{\alpha}^{2} + w_{\alpha}^{2}}{2} + g \cos \theta \left( \tilde{b} + b + \frac{h}{2} \right) \right),$$

and summing up equations (39) and (40), after some straightforward computations we obtain

$$\partial_t E_{\alpha} + \partial_x \left( E_{\alpha} u_{\alpha} + g \cos \theta h_{\alpha} \frac{h}{2} u_{\alpha} \right) + \frac{1}{\rho} P_{NH,\alpha} = M T_{\alpha} + \frac{1}{\rho} u_{\alpha} (K_{\alpha-1/2} - K_{\alpha+1/2}) + \frac{g \cos \theta}{2} (h_{\alpha} \partial_t h - h \partial_t h_{\alpha}), \quad (41)$$

where we have used that

$$g \cos \theta (\tilde{b} + b + h) \partial_t h_{\alpha} = \partial_t \left( h_{\alpha} g \cos \theta \left( \tilde{b} + b + \frac{h}{2} \right) \right) + \frac{g \cos \theta}{2} (h \partial_t h_{\alpha} - h_{\alpha} \partial_t h).$$

In equation (41), $P_{NH,\alpha}$ and $MT_{\alpha}$ collect the terms related to non-hydrostatic pressure and mass transference, respectively, in the layer $\Omega_{\alpha}$. The non-hydrostatic contributions are

$$\begin{aligned} P_{NH,\alpha} &= u_{\alpha} \partial_x (h_{\alpha} p_{\alpha}^{nh}) + p_{\alpha+1/2}^{nh} (-u_{\alpha} \partial_x z_{\alpha+1/2} + w_{\alpha}) - p_{\alpha-1/2}^{nh} (-u_{\alpha} \partial_x z_{\alpha-1/2} + w_{\alpha}) \\ &= \partial_x (h_{\alpha} u_{\alpha} p_{\alpha}^{nh}) - p_{\alpha+1/2}^{nh} \sum_{\beta=1}^{\alpha} \partial_x (h_{\beta} u_{\beta}) + p_{\alpha-1/2}^{nh} \sum_{\beta=1}^{\alpha-1} \partial_x (h_{\beta} u_{\beta}) \end{aligned}$$

where we have used restriction (4) and (9d). Notice that in the particular cases $\alpha = 1$ and $\alpha = N$, using the non-penetration condition and $p_{|b+h}^{nh} = 0$, we have

$$P_{NH,1} = \partial_x (h_1 u_1 p_1^{nh}) - p_{3/2}^{nh} \partial_x (h_1 u_1), \quad P_{NH,N} = \partial_x (h_N u_N p_N^{nh}) + p_{N-1/2}^{nh} \sum_{\beta=1}^{N-1} \partial_x (h_\beta u_\beta).$$

It leads to a conservative contribution since the terms on $p_{\alpha \pm 1/2}^{nh}$ vanish when summing up from $\alpha = 1, \dots, N$:

$$\sum_{\beta=1}^{N} P_{NH,\beta} = \sum_{\beta=1}^{N} \partial_x (h_{\beta} p_{\beta}^{nh} u_{\beta}).$$

Concerning the mass transfer terms, they are

$$MT_{\alpha} = g \cos \theta (\tilde{b} + b + h) (G_{\alpha+1/2} - G_{\alpha-1/2}) + \frac{G_{\alpha+1/2}}{2} [u_{\alpha} u_{\alpha+1} + w_{\alpha} w_{\alpha+1}] - \frac{G_{\alpha-1/2}}{2} [u_{\alpha} u_{\alpha-1} + w_{\alpha} w_{\alpha-1}].$$

It is easy to see that, since we have set $G_{1/2} = G_{N+1/2} = 0$, it holds

$$\sum_{\beta=1}^{N} MT_{\beta} = 0, \quad \text{and} \quad \frac{g}{2} \sum_{\beta=1}^{N} (h_{\alpha} \partial_t H - H \partial_t h_{\alpha}) = 0,$$

and thus, summing up (41) in all the layers, the proof is finished.

Finally, by using the definitions of $K_{\alpha+1/2}$ for $\alpha = 0, \dots, N$, the viscous terms satisfy

$$\begin{aligned} \sum_{\beta=1}^{N} u_{\alpha} (K_{\beta-1/2} - K_{\beta+1/2}) &= u_1 K_{1/2} - u_N K_{N+1/2} + \sum_{\beta=1}^{N-1} K_{\beta+1/2} (u_{\beta+1} - u_{\beta}) \\ &= -\mu(I_{1/2})p_{1/2} \frac{u_1^2}{|u_1|} - \sum_{\beta=1}^{N-1} \frac{\eta_{\beta+1/2}}{2} \frac{(u_{\beta+1} - u_{\beta})^2}{h_{\beta+1/2}}, \end{aligned}$$

that is then a dissipative contribution, as expected.

□
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és Madroñal, 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 $\mu(\cdot)$-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 $\mu(\cdot)$-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 $\mu(\cdot)$-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äuser Basel, 2004. 14, 18

[10] F. Bouchut, J. M. Delgado-Sánchez, E. D. Fernández-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ández-Nieto, E. H. Koné, A. Mangeney, and G. Narbona-Reina. Dilatancy in dry granular flows with a compressible $\mu(\cdot)$ rheology. *Journal of Computational Physics*, 429:110013, mar 2021. 2

[12] F. Bouchut, E. D. Fernández-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éorie 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ément, and B. Andreotti. Nonlocal rheology of granular flows across yield conditions. *Physical Review Letters*, 111(23), dec 2013. 2

[17] F. Boyer, É. 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
[19] M. Brunet, L. Moretti, A. Le Friant, A. Mangeney, E. D. Fernández Nieto, and F. Bouchut. Numerical simulation of the 30–45 ka debris avalanche flow of Montagne Pelée 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ález-Vida, and C. Parés. 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ández-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édale. A three-dimensional numerical model for dense granular flows based on the $\mu(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öner, 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ández-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ánchez, E. D. Fernández-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ández-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ández-Nieto, J. Garres-Díaz, A. Mangeney, and G. Narbona-Reina. A multilayer shallow model for dry granular flows with the $\mu(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ández-Nieto, J. Garres-Díaz, A. Mangeney, and G. Narbona-Reina. 2D granular flows with the $\mu(I)$ rheology and side walls friction: a well-balanced multilayer discretization. *Journal of Computational Physics*, 356:192–219, 2018. 2, 14, 15, 31
[38] E. D. Fernández-Nieto, J. Garres-Díaz, and P. Vigneaux. Multilayer models for hydrostatic Herschel-Bulkley viscoplastic flows. *Computers & Mathematics with Applications*, 139:99–117, jun 2023. 24

[39] E. D. Fernández-Nieto, E. H. Koné, and T. Chacón Rebollo. A multilayer method for the hydrostatic Navier-Stokes equations: A particular weak solution. *Journal of Scientific Computing*, 60(2):408–437, nov 2013. 2, 6, 11

[40] E. D. Fernández-Nieto, M. Parisot, Y. Penel, and J. Sainte-Marie. A hierarchy of dispersive layer-averaged approximations of Euler equations for free surface flows. *Communications in Mathematical Sciences*, 16(5):1169–1202, 2018. 1, 3, 7, 10, 11, 38

[41] I.A. Frigaard and C. Nouar. On the usage of viscosity regularisation methods for visco-plastic fluid flow computation. *Journal of Non-Newtonian Fluid Mechanics*, 127(1):1–26, apr 2005. 19

[42] J. Garres-Díaz, E. D. Fernández-Nieto, A. Mangeney, and T. Morales de Luna. A weakly non-hydrostatic shallow model for dry granular flows. *Journal of Scientific Computing*, 86(2), jan 2021. 1, 3, 4, 5, 9, 11, 13, 14, 15, 19, 22, 24, 25, 35

[43] GDR MiDi. On dense granular flows. *The European Physical Journal E*, 14(4):341–365, 2004. 1, 2

[44] L. Gesenhues, J. J. Camata, A. M.A. Côrtes, F. A. Rochinha, and A.L.G.A. Coutinho. Finite element simulation of complex dense granular flows using a well-posed regularization of the $\mu(I)$-rheology. *Computers & Fluids*, 188:102–113, jun 2019. 2

[45] R. Glowinski, J.-L. Lions, and R. Trémolières. *Analyse numérique des inéquations variationnelles* (Tomes 1 et 2). Bordas (Dunod), 1976. 4

[46] J. M. N. T. Gray and A. N. Edwards. A depth-averaged $\mu(I)$-rheology for shallow granular free-surface flows. *Journal of Fluid Mechanics*, 755:503–534, aug 2014. 2, 5

[47] G. Grosso, M. Antuono, and M. Brocchini. Dispersive nonlinear shallow-water equations: some preliminary numerical results. *Journal of Engineering Mathematics*, 67(1-2):71–84, September 2009. 3

[48] E. Guazzelli and O. Pouliquen. Rheology of dense granular suspensions. *Journal of Fluid Mechanics*, 852:P1, 2018. 1

[49] J.-L. Guermond, B. Popov, E. Tovar, and C. Kees. Robust explicit relaxation technique for solving the Green-Naghdi equations. *Journal of Computational Physics*, 399:108917, December 2019. 3

[50] J.L. Guermond, P. Minev, and J. Shen. An overview of projection methods for incompressible flows. *Computer Methods in Applied Mechanics and Engineering*, 195(44-47):6011–6045, sep 2006. 3, 13

[51] J. Heyman, R. Delannay, H. Tabuteau, and A. Valance. Compressibility regularizes the $\mu(I)$-rheology for dense granular flows. *Journal of Fluid Mechanics*, 830:553–568, oct 2017. 2

[52] I. R. Ionescu, A. Mangeney, F. Bouchut, and R. Roche. Viscoplastic modeling of granular column collapse with pressure-dependent rheology. *Journal of Non-Newtonian Fluid Mechanics*, 219(0):1–18, 2015. 1, 2, 19, 30

[53] P. Jop, Y. Forterre, and O. Pouliquen. A constitutive law for dense granular flows. *Nature*, 441(7094):727–730, jun 2006. 1, 4

[54] P. Jop, Y. Forterre, and O. Pouliquen. Initiation of granular surface flows in a narrow channel. *Physics of Fluids*, 19(8):088102, aug 2007. 1

[55] P.-Y. Lagrée, L. Staron, and S. Popinet. The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes with a $\mu(I)$-rheology. *Journal of Fluid Mechanics*, 686:378–408, sep 2011. 1, 2, 4

[56] C. Lusso, F. Bouchut, A. Ern, and A. Mangeney. A free interface model for static/flowing dynamics in thin-layer flows of granular materials with yield: Simple shear simulations and comparison with experiments. *Applied Sciences*, 7(4), 2017. 2

[57] C. Lusso, F. Bouchut, A. Ern, and A. Mangeney. Explicit solutions to a free interface model for the static/flowing transition in thin granular flows. *ESAIM: Mathematical Modelling and Numerical Analysis*, 55:S369–S395, 2021. 2
[58] C. Lusso, A. Ern, F. Bouchut, A. Mangeney, M. Farin, and O. Roche. Two-dimensional simulation by regularization of free surface viscoplastic flows with Drucker-Prager yield stress and application to granular collapse. *Journal of Computational Physics*, 333:387-408, 2017. 19, 24

[59] J. Macías, C. Escalante, and M.J. Castro. Multilayer-HySEA model validation for landslide generated tsunamis. Part I Rigid slides. *Nat. Hazards Earth Syst. Sci. Discuss.*, 2021. 18

[60] A. Mangeney, F. Bouchut, N. Thomas, J. P. Vilotte, and M. O. Bristeau. Numerical modeling of self-channeling granular flows and of their levee-channel deposits. *Journal of Geophysical Research: Earth Surface*, 112(F2):n/a-n/a, 2007. F02017. 2, 24

[61] A. Mangeney, O. Roche, O. Hungr, N. Mangold, G. Faccaroni, and A. Lucas. Erosion and mobility in granular collapse over sloping beds. *Journal of Geophysical Research*, 115(F3), sep 2010. 2, 26

[62] A. Mangeney-Castelnau, J.-P. Vilotte, M. O. Bristeau, B. Perthame, F. Bouchut, C. Simeoni, and S. Yerneni. Numerical modeling of avalanches based on Saint Venant equations using a kinetic scheme. *Journal of Geophysical Research: Solid Earth*, 108(B11), 2003. eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2002JB002024. 2

[63] N. Martin, I. R. Ionescu, A. Mangeney, F. Bouchut, and M. Farin. Continuum viscoplastic simulation of a granular column collapse on large slopes: $\mu(I)$ rheology and lateral wall effects. *Physics of Fluids*, 29(1):013301, 2017. 1, 2, 29, 32

[64] Eike Müller, Xu Guo, Robert Scheichl, and Sinan Shi. Matrix-free GPU implementation of a preconditioned conjugate gradient solver for anisotropic elliptic PDEs. *Comput. Vis. Sci.*, 16(2):41-58, apr 2013. 18

[65] M. Pailha, M. Nicolas, and O. Pouliquen. Initiation of underwater granular avalanches: Influence of the initial volume fraction. *Physics of Fluids*, 20(11):111701, 2008. 1

[66] T.C. Papanastasiou. Flows of Materials with Yield. *Journal of Rheology*, 31(5):385-404, 1987. 4

[67] D.H. Peregrine. Long waves on a beach. *J. Fluid Mech.*, 27(04):815-827, 1967. 2

[68] M. Peruzzetto, J.-C. Komorowski, A. Le Friant, M. Rosas-Carbajal, A. Mangeney, and Y. Legendre. Modeling of partial dome collapse of La Soufrière of Guadeloupe volcano: implications for hazard assessment and monitoring. *Scientific Reports*, 9(1), sep 2019. 25

[69] M. Peruzzetto, C. Levy, Y. Thiery, G. Grandjean, A. Mangeney, A.-M. Lejeune, A. Nachbaur, Y. Legendre, Benoit Vittecoq, Jean-Marie Saurel, Valérie Clouard, Thomas Dewez, Fabrice R. Fontaine, M. Mergili, S. Lagarde, J.C. Komorowski, A. Le Friant, and A. Lemarchand. Simplified simulation of rock avalanches and subsequent debris flows with a single thin-layer model: Application to the Prêcheur river (Martinique, Lesser Antilles). *Engineering Geology*, 296:106457, jan 2022. 25

[70] M. Peruzzetto, A. Mangeney, F. Bouchut, G. Grandjean, C. Levy, Y. Thiery, and A. Lucas. Topography curvature effects in thin-layer models for gravity-driven flows without bed erosion. *Journal of Geophysical Research: Earth Surface*, 126(4), apr 2021. 24

[71] M. Peruzzetto, A. Mangeney, G. Grandjean, C. Levy, Y. Thiery, J. Rohmer, and A. Lucas. Operational Estimation of Landslide Runout: Comparison of Empirical and Numerical Methods. *Geosciences*, 10(11):424, oct 2020. 25

[72] M. Pirulli and A. Mangeney. Results of back-analysis of the propagation of rock avalanches as a function of the assumed rheology. *Rock Mechanics and Rock Engineering*, 41(1):59-84, 2008. 2

[73] O. Pouliquen. On the shape of granular fronts down rough inclined planes. *Physics of Fluids*, 11(7):1956-1958, 1999. 2

[74] O. Pouliquen. Scaling laws in granular flows down rough inclined planes. *Physics of Fluids*, 11(3):542-548, mar 1999. 2

[75] O. Pouliquen and Y. Forterre. Friction law for dense granular flows: application to the motion of a mass down a rough inclined plane. *Journal of Fluid Mechanics*, 453:133-151, feb 2002. 2
[76] O. Pouliquen and Y. Forterre. A non-local rheology for dense granular flows. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 367(1909):5091–5107, dec 2009. 2

[77] M. Rauter. The compressible granular collapse in a fluid as a continuum: validity of a Navier–Stokes model with $\mu(I)$-rheology. *Journal of Fluid Mechanics*, 915, mar 2021. 2

[78] M. Rauter, T. Barker, and W. Fellin. Granular viscosity from plastic yield surfaces: The role of the deformation type in granular flows. *Computers and Geotechnics*, 122:103492, jun 2020. 1

[79] S. B. Savage and K. Hutter. The motion of a finite mass of granular material down a rough incline. *Journal of Fluid Mechanics*, 199:177–215, 1989. 2, 11

[80] D. G. Schaeffer, T. Barker, D. Tsuji, P. Gremaud, M. Shearer, and J. M. N. T. Gray. Constitutive relations for compressible granular flow in the inertial regime. *Journal of Fluid Mechanics*, 874:926–951, jul 2019. 2

[81] M. Thomadakis and M. Leschziner. A pressure-correction method for the solution of incompressible viscous flows on unsstructured grids. *International Journal for Numerical Methods in Fluids*, 22(7):581–601, April 1996. 3

[82] Y. Yamazaki, Z. Kowalik, and K. F. Cheung. Depth-integrated, non-hydrostatic model for wave breaking and run-up. *Numerical Methods in Fluids*, 61:473–497, 2008. 2, 3