High-order well-balanced schemes for shallow models for dry avalanches

Appl. Numer. Math. 215 (2025) 138–156 Available online 17 April 2025 0168-9274/© 2025 The Author(s). Published by Elsevier B.V. on behalf of IMACS. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Applied Numerical Mathematics journal homepage: www.elsevier.com/locate/apnum Research Paper High-order well-balanced schemes for shallow models for dry avalanches M.J. Castro Díaza, C. Escalanteb,, J. Garres-Díazc, ,∗, T. Morales de Lunaa, aDpto. Análisis Matemático, Estad. e I.O. y Matemática Aplicada, Universidad de Málaga, 29071 Málaga, Spain bDpto. Matemática Aplicada, Universidad de Málaga, 29071 Málaga, Spain cDpto. Matemática Aplicada II, Universidad de Sevilla, 41092 Sevilla, Spain A R T I C L E I N F O A B S T R A C T Keywords: Finite volume method Well-balanced schemes High-order scheme Savage-Hutter model Granular flows In this work we consider a depth-averaged model for granular flows with a Coulomb-type friction force described by the 𝜇(𝐼)rheology. In this model, the so-called lake-at-rest steady states are of special interest, where velocity is zero and the slope is under a critical threshold defined by the angle of repose of the granular material. It leads to a family with an infinite number of lakeat-rest steady states. We describe a well-balanced reconstruction procedure that allows to define well-balanced finite volume methods for such problem. The technique is generalized to highorder space/time schemes. In particular, the second and third-order schemes are considered in the numerical tests section. An accuracy test is included showing that second and third-order are achieved. A well-balanced test is also considered. The proposed scheme is well-balanced for steady states with non-constant free surface, and it is exactly well-balanced for those steady states given by a simple characterization. 1. Introduction Geophysical flows have been deeply studied in last decades due to their importance designing early warning systems or studies against natural hazards (floods, landslides, tsunamis, snow avalanches, volcanic eruptions, etc.). The mathematical modeling and numerical simulation of such flows is a very active research topic, especially since popularization of advanced technological resources, such as high-performance computers. Nevertheless, many challenges remain: the physical understanding and the design of complex rheologies; the development of sophisticated models whose associated computational cost is reasonable for practical purposes; the design of efficient numerical schemes to solve complex models including viscous (viscoplastic, viscoelastic,...) terms. The complexity of these flows requires progress on all these research lines, where different scientists (geophysicists, mathematicians, computer scientists,...) may contribute with valuable advances. Here, we focus on granular flows, concretely aerial avalanches, and we tackle the task of defining efficient numerical schemes for some depth-averaged avalanche models previously introduced in the literature. Concerning the modeling of aerial avalanches, one of the most accepted rheological laws describing their dynamics is the so-called 𝜇(𝐼)-rheology (see [1]). This rheology considers a variable friction coefficient depending on the velocity and pressure. It is based on a Drucker-Prager plasticity criterion defining the deviatoric tensor as *Corresponding author. E-mail address: [email protected] (J. Garres-Díaz). https://doi.org/10.1016/j.apnum.2025.04.008 Received 17 September 2024; Received in revised form 17 March 2025; Accepted 15 April 2025 Applied Numerical Mathematics 215 (2025) 138–156 139 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. ⎧ ⎪ ⎨ ⎪ ⎩ 𝝉=𝜇(𝐼)𝑝 ‖𝑫‖ 𝑫if ‖𝑫‖≠0, ‖𝝉‖≤𝜇𝑠𝑝if ‖𝑫‖=0, (1) where 𝑫is the strain rate tensor, 𝑝is the pressure and 𝜇𝑠=tan𝜃0, with 𝜃0the angle of repose of the granular material. This rheology was considered in [2] where it was included in a continuous 2D Navier-Stokes solver and the granular material was considered as a fluid with a 𝜇(𝐼)-viscosity. In order to avoid the singularity when the strain rate vanishes, the authors considered a regularization of the 𝜇(𝐼)-viscosity. The regularization technique was compared with the Augmented Lagrangian (AL) method in [3] for these kind of flows. It was shown that, for a proper choice of the regularization parameter and with a suitable finite element discretization, that the error associated to this regularization is negligible with respect to the time/space discretization errors. Other authors have considered finite element discretizations combined with AL methods to simulate granular flows and analyze different aspects associated to this rheology, including comparisons with laboratory experiments (see e.g. [4,5]). However, the computational cost associated to these discretizations is very high. Therefore, simplified models, based on a depth-averaging process, have become very popular for these flows. This type of models will be the starting point for this work. Let us remark that some other more complex models have been proposed in the last decade for granular flows. For instance, in [6–9] one can find models with non-local effects or a compressible version of the 𝜇(𝐼)-rheology, where the solid volume fraction is no more constant. As previously said, depth-averaged models have been widely used to study granular flows. Since the pioneering work of SavageHutter [10], where the friction force with the bottom is modeled through a Coulomb friction term, several more sophisticated depthaveraged models have been introduced. For instance, Gray and Edwards [11] proposed a shallow model, which includes secondorder viscous terms, for dry granular flows with 𝜇(𝐼)-rheology (see also [12]). An extension of this model, neglecting second-order viscous terms and using the multilayer (or layer-averaged) framework, has been proposed in [13]. Moreover, dispersive effects can be incorporated as in [14], where a weakly non-hydrostatic shallow model was proposed, which takes into account the contribution of the vertical acceleration of the flow. This model was also extended to the multilayer case in [15]. Most of these models were obtained using a tilted reference system over a plane with constant slope and performing an average on the normal direction. This means that in these models the average horizontal component of the velocity is measured in the downslope direction, which is an advantage with respect to Cartesian models. Moreover, when validating such models, some of the laboratory experiments used can be only reproduced for these local models due to the definition of their initial conditions. Interestingly, in [16] authors showed that the motion criterion is wrong in all previous local models due to the change of variables (to tilted coordinates), and they proposed a correction to properly capture this motion criterion. From the numerical point of view, one of the major difficulties of discretizing these models is the numerical treatment of the multi-evaluated stress tensor (1)when the velocity vanishes (case 𝑫=0). On the one hand, in the motion phase of the flow (|𝑢|>0), (semi-)implicit discretizations of the viscous term have been considered for friction terms in order to avoid restrictive time-steps (see e.g. [17–19]). On the other hand, different treatments have been considered to deal with the multi-evaluated definition of the stress tensor when the flow is at rest (𝑢=0). In [13,15] a regularization technique was considered, which makes the velocity tend to zero in stop-motion situations, although null velocity cannot be obtained, since a residual velocity remains depending on the regularization parameter. In [14,20] a different treatment is used based on the interpretation of the friction force made in [21]. It assumes that the friction force opposes to the motion of the flow, but it cannot change the sign of the velocity. This assumption means that the friction force should be bounded in practice. Although it is difficult to apply this technique to multilayer systems, it is a valid approach for single-layer models. It has the advantage that a zero velocity is exactly recovered when the material stops and remains at rest. Moreover, this leads to an explicit discretization of the model with a posteriori (friction) correction of the velocity. In particular, a similar approach will be used in this work. An important property concerning numerical schemes for geophysical flows is the so-called well-balanced property. Indeed, many practical applications of fluid or granular flows correspond to small perturbation of steady states. In such situations, the use of a well-balanced numerical method is mandatory since the amplitudes may be of the order (or bigger than) the truncation error of the scheme. This is for instance the case of a dry avalanche propagating over a slope which is initially steady or when it converges to a steady state after a large time evolution. Remark that it is not always possible to refine the mesh so that the truncation error of the method is lower than this perturbation. Therefore, it is important to preserve, all or a subfamily of, the steady states of the model. See for instance the works [22–28] for balance laws, in particular, shallow water systems, and [29,20,30–34] for avalanche and granular models, among many others. Focusing on granular flows models, which is the topic of this work, the so called lake-at-rest are of special interest. These equilibria correspond to situations where the velocity is zero and the non-flat free surface of the material has a slope smaller than a critical value fixed by the repose angle. Remark that, in general, such steady states are characterized by a inequality instead of an equality. This fact makes it difficult to design well-balanced numerical schemes preserving such steady solutions. In previous works ([31,14,15]) first-order well-balanced schemes were proposed. The well-balanced property was achieved by means of a modified hydrostatic reconstruction which includes the Coulomb friction term (see [35]). In [20] a second-order MUSCL reconstruction combined with central-upwind scheme was proposed for a one-dimensional Savage-Hutter type model of submarine landslides and generated tsunami waves. The Coulomb friction term was treated using [29] and the proposed scheme is well-balanced for lake-at-rest steady states However, designing high-order well-balanced schemes is not a trivial task and, up to our knowledge, there are no previous works showing a high order well-balanced scheme for granular or complex flows including such a friction coefficient. Higher-order schemes capture gradients and variations in the solution more accurately than low-order methods. They reduce numerical dissipation and dispersion errors, which are crucial for accurately capturing shocks, discontinuities, and smooth solution features. Moreover, for a Applied Numerical Mathematics 215 (2025) 138–156 140 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. given level of accuracy, high-order schemes often require fewer grid points compared to low-order methods. Therefore, high-order schemes allow to reduce computational costs, especially in large-scale simulations, which results in a better efficiency. The goal of this work is to propose a high order well-balanced discretization for shallow models for granular flows including the multi-evaluated definition of the stress tensor when the flow is at rest. Based on the ideas in [36], we develop a novel procedure to obtain a well-balanced scheme for lake-at-rest solutions this kind of granular flows. As already mentioned, the main difficulty in this case lies in the expression of the stationary solutions, which are defined in terms of an inequality. The main novelty of this paper is to describe a procedure to obtain a well-balanced scheme for granular flows, which can be easily generalized to the arbitrary high-order case. The paper is organized as follows: Section 2is devoted to the description of the model and to introduce the basic notation used in the paper. Section 3will describe the numerical scheme proposed here. In particular, in Subsection 3.1 the general procedure to obtain a well-balanced scheme based on a reconstruction operator is presented, and Subsection 3.2 is devoted to describe the reconstruction procedure for the granular model. Subsection 3.3 deals with the time discretization. Some numerical tests will be presented in Section 4. In particular, we will show that the expected order of accuracy is achieved for the second and third-order case, and also that the scheme is indeed well-balanced. Finally, the conclusions are presented in Section 5. The paper is completed with Appendix A, where some particular cases of steady states are studied. 2. Shallow model for granular flows For the sake of simplicity and clarity in the exposition, we focus here on the hydrostatic model introduced in [14] for granular flows, which coincides in fact with the one in [11] when the second-order viscous term is neglected. We will consider Cartesian coordinates, although everything can be trivially extended to the tilted coordinates case. Adding non-hydrostatic effects is no major problem either, following [14]. Thus, we consider a bidimensional domain, (𝑥,𝑧)∈ℝ2, and denote by ℎ(𝑥),𝑢(𝑥)∈ℝthe total fluid depth and the depth-integrated horizontal velocity, respectively. For simplicity, we consider that the bottom topography, 𝑏(𝑥), is a continuous and differentiable function. The hydrostatic model for granular flows is then written as 𝜕𝑡ℎ+𝜕𝑥(ℎ𝑢)=0, 𝜕𝑡(ℎ𝑢)+𝜕𝑥(ℎ𝑢2+1 2𝑔ℎ2)+𝑔ℎ𝜕𝑥𝑏=−𝜏𝑥𝑧|𝑏∕𝜌, where 𝜌is the density of the material, assumed to be constant, 𝜏𝑥𝑧 is the deviatoric stress tensor, which includes the friction with the bottom, defined by 𝜏𝑥𝑧|𝑏=⎧ ⎪ ⎨ ⎪ ⎩ 𝜌𝑔ℎ𝜇|𝑏 𝑢 |𝑢|if |𝑢|≠0, |||𝜏𝑥𝑧|𝑏|||≤𝜌𝑔ℎ𝜇𝑠if |𝑢|=0, (2a) where 𝜇is the variable friction coefficient 𝜇|𝑏=𝜇𝑠+𝜇2−𝜇𝑠 𝐼0+𝐼|𝑏 𝐼|𝑏, with 𝐼|𝑏= 𝑑𝑠|||(𝜕𝑧𝑢)|𝑏||| √𝜑𝑠𝑔ℎ , and (𝜕𝑧𝑢)|𝑏=𝑢 ℎ,(2b) being 𝑑𝑠the mean grain size and 𝜇2,𝜇𝑠,𝐼0,𝜑𝑠constant rheological parameters depending on the considered granular material. When looking for the steady solutions of system (2), we focus on lake-at-rest solutions, that is, those with zero velocity. Then, looking at the momentum conservation equation with 𝑢=0, and denoting by 𝜂=𝑏+ℎthe free surface level, we get 𝑔ℎ||𝜕𝑥𝜂||=|||𝜏𝑥𝑧|𝑏,𝑢=0 ||| 𝜌 ≤𝑔ℎ𝜇𝑠. Therefore, the steady states corresponding to water at rest, are given by 𝑢=0, and ||𝜕𝑥𝜂||≤𝜇𝑠,(3) where we recall that 𝜇𝑠=tan𝜃0, with 𝜃0the angle of repose of the material. We will now proceed to design a high order numerical scheme for (2)that preserves lake-at-rest steady states in a sense to be precised. 3. Numerical scheme First of all, let us write system (2)in a compact form 𝜕𝑡𝑼+𝜕𝑥𝑭(𝑼)+𝑺(𝑼)𝜕𝑥𝑏=𝑻(𝑼),(4) where Applied Numerical Mathematics 215 (2025) 138–156 141 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. 𝑼=(ℎ ℎ𝑢), 𝑭(𝑼)=(ℎ𝑢 ℎ𝑢2+1 2𝑔ℎ2), 𝑺(𝑼)=(0 𝑔ℎ), 𝑻(𝑼)=(0 −𝜏𝑥𝑧|𝑏∕𝜌). Recall that we are assuming here a continuous and differentiable topography 𝑏(𝑥). The extension to the case of a discontinuous bottom could be addressed following the ideas of [36,37]. We consider now, as usual in finite volume schemes, a spatial discretization into control volumes or cells 𝑉𝑖=[𝑥𝑖−1∕2,𝑥𝑖+1∕2], for 𝑖∈={1,...,𝑁}. For the sake of simplicity, we consider a uniform discretization with cell sizes Δ𝑥=𝑥𝑖+1∕2 −𝑥𝑖−1∕2. Let us denote by 𝑥𝑖=(𝑥𝑖−1∕2 +𝑥𝑖+1∕2)∕2 the center of 𝑉𝑖. We keep, for now, the time variable continuous and denote by 𝑼𝑖(𝑡)the cell average of 𝑼(𝑥,𝑡), that is, for 𝑖∈: 𝑼𝑖(𝑡)= 1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 𝑼(𝑥,𝑡) 𝑑𝑥. (5) Following [38,28,37], given a sequence of states {𝑼𝑖(𝑡)}𝑖∈, we consider, for every cell 𝑉𝑖, a 𝑠-order reconstruction operator 𝑷𝑡 𝑖(𝑥)=𝑷𝑖(𝑥;{𝑼𝑗(𝑡)}𝑗∈𝑖), where 𝑗is the stencil of cells of dependence for the operator, satisfying 1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 𝑷𝑡 𝑖(𝑥) 𝑑𝑥 =𝑼𝑖(𝑡)=𝑼(𝑥𝑖,𝑡)+(Δ𝑥𝑠). We then consider the reconstructed values at the inter-cells 𝑼𝑡,+ 𝑖−1∕2 = lim 𝑥→𝑥+ 𝑖−1∕2 𝑷𝑡 𝑖(𝑥), 𝑼− 𝑖+1∕2 = lim 𝑥→𝑥− 𝑖+1∕2 𝑷𝑡 𝑖(𝑥), which are assumed to be approximations of order 𝑠of the solution at the interfaces: 𝑼𝑡,± 𝑖+1∕2 =𝑼(𝑥𝑖+1∕2,𝑡)+(Δ𝑥𝑠). Then, a semi-discrete high-order path conservative method is written as 𝑑 𝑑𝑡𝑼𝑖=− 1 Δ𝑥(𝑖+1∕2(𝑡)−𝑖−1∕2(𝑡) + 𝑺𝑖−𝑻𝑖),(6) with 𝑖+1∕2(𝑡)=𝔽(𝑼𝑡,− 𝑖+1∕2,𝑼𝑡,+ 𝑖+1∕2), where 𝔽(⋅,⋅)is a first order consistent numerical flux (𝔽(𝑼,𝑼)=𝑭(𝑼)) and 𝑺𝑖−𝑻𝑖≈ 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 𝑺(𝑷𝑡 𝑖(𝑥))𝜕𝑥𝑏 𝑑𝑥 − 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 𝑻(𝑷𝑡 𝑖(𝑥)) 𝑑𝑥. In practice, we consider in this work a numerical flux, which may be written as (see [39]): 𝔽(𝑼𝑡,− 𝑖+1∕2,𝑼𝑡,+ 𝑖+1∕2)=1 2[𝑭(𝑼𝑡,+ 𝑖+1∕2)+𝑭(𝑼𝑡,− 𝑖+1∕2)−(𝛽0(𝑼𝑡,+ 𝑖+1∕2 −𝑼𝑡,− 𝑖+1∕2)+𝛽1(𝑭(𝑼𝑡,+ 𝑖+1∕2)−𝑭(𝑼𝑡,− 𝑖+1∕2)))].(7) Notice that the definition of the numerical scheme (6)assumes a continuous bottom (or a continuous reconstruction of the bottom). Otherwise, a contribution of the jumps related to the discontinuities of the topography at the interfaces should be taken into account. 3.1. Well-balanced scheme We now adapt scheme (6)to preserve the lake-at-rest equilibria. To do so, let us recall first the definitions introduced in [40]: we shall distinguish between exactly well-balanced or well-balanced schemes. The former preserves the cell-averages of the exact steady state while the later preserves the cell-averages of a discrete approximation of the stationary solutions. More precisely, Definition 1. Let 𝑼𝑒(𝑥)be any stationary solution of (4)and the corresponding cell averages {𝑼𝑒 𝑖}𝑖∈obtained using (5)(or any given high order quadrature formula). Then, a numerical method written in the form (6)is said to be: •exactly well-balanced for 𝑼𝑒(𝑥)if the cell-averages {𝑼𝑒 𝑖}𝑖∈are an equilibrium of the ODE system (6). •well-balanced with order 𝑟≥𝑠if, for every Δ𝑥, there exists a discrete approximation 𝑼𝑒 Δ𝑥,𝑖 =𝑼𝑒 𝑖+(Δ𝑥𝑟) such that {𝑼𝑒 Δ𝑥,𝑖}𝑖∈is an equilibrium of the ODE system (6). Applied Numerical Mathematics 215 (2025) 138–156 142 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. In order to get well-balanced schemes, either exactly or not, we shall describe a general algorithm based on a reconstruction procedure following the ideas proposed in [36]. Let us start by assuming that for any time 𝑡, given the sequence {𝑼𝑖(𝑡)}, we select for every control volume 𝑉𝑖a continuous stationary solution 𝑼𝑡,𝑒 𝑖(𝑥)of system (4)satisfying the conservation property 1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 𝑼𝑡,𝑒 𝑖(𝑥) 𝑑𝑥 =𝑼𝑖(𝑡). In particular, it holds 𝜕𝑥(𝑭(𝑼𝑡,𝑒 𝑖))+𝑺(𝑼𝑡,𝑒 𝑖)𝜕𝑥𝑏=𝑻(𝑼𝑡,𝑒 𝑖).(8) Remark that although stationary solutions do not depend on time, the selected stationary solution at each time may not be necessarily the same, hence the notation 𝑼𝑡,𝑒(𝑥). Remark 1. Notice that one of the main difficulties is to give an explicit expression of the term 𝑻(𝑼𝑡,𝑒 𝑖). Since here we will focus on lake-at-rest steady states (3), the operator 𝑻is multi-evaluated and a priori we only know that |||𝜏𝑥𝑧|𝑏∕𝜌|||≤𝑔ℎ𝜇𝑠. Nevertheless, since 𝑼𝑡,𝑒 𝑖should satisfy (8), this will actually fix the definition of 𝑻(𝑼𝑡,𝑒 𝑖). More explicitly, the friction should balance (8)and therefore we have: 𝑻(𝑼𝑡,𝑒 𝑖(𝑥))=(0,𝑔ℎ𝑡,𝑒 𝑖(𝑥)𝜕𝑥𝜂𝑡,𝑒 𝑖(𝑥))𝑇. Then, any semi-discrete scheme (6)can be rewritten as 𝑑 𝑑𝑡𝑼𝑖=− 1 Δ𝑥(𝑖+1∕2(𝑡)−𝑭(𝑼𝑡,𝑒 𝑖(𝑥𝑖+1∕2)) − 𝑖−1∕2(𝑡)+𝑭(𝑼𝑡,𝑒 𝑖(𝑥𝑖−1∕2))) −1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 (𝑺(𝑷𝑡 𝑖(𝑥)) − 𝑺(𝑼𝑡,𝑒 𝑖(𝑥)))𝜕𝑥𝑏 𝑑𝑥 + 1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 (𝑻(𝑷𝑡 𝑖(𝑥)) − 𝑻(𝑼𝑡,𝑒 𝑖(𝑥)))𝑑𝑥, where we have set 𝑺𝑖−𝑻𝑖=−𝑭(𝑼𝑡,𝑒 𝑖(𝑥𝑖+1∕2))+𝑭(𝑼𝑡,𝑒 𝑖(𝑥𝑖−1∕2))+ 1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 (𝑺(𝑷𝑡 𝑖(𝑥)) − 𝑺(𝑼𝑡,𝑒 𝑖(𝑥)))𝜕𝑥𝑏 𝑑𝑥 − 1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 (𝑻(𝑷𝑡 𝑖(𝑥)) − 𝑻(𝑼𝑡,𝑒 𝑖(𝑥)))𝑑𝑥. In practice, the integrals of the source terms will be computed by means of a quadrature formula with order 𝑟≥𝑠, with 𝑠is the order of the reconstruction operator: 𝑑 𝑑𝑡𝑼𝑖=− 1 Δ𝑥(𝑖+1∕2(𝑡)−𝑭(𝑼𝑡,𝑒 𝑖(𝑥𝑖+1∕2)) − 𝑖−1∕2(𝑡)+𝑭(𝑼𝑡,𝑒 𝑖(𝑥𝑖−1∕2))) − 𝑘 ∑ 𝑗=1 𝜔𝑖 𝑗(𝑺(𝑷𝑡 𝑖(𝑥𝑖 𝑗)) − 𝑺(𝑼𝑡,𝑒 𝑖(𝑥𝑖 𝑗)))(𝜕𝑥𝑏)𝑥=𝑥𝑖 𝑗+ 𝑘 ∑ 𝑗=1 𝜔𝑖 𝑗(𝑻(𝑷𝑡 𝑖(𝑥𝑖 𝑗)) − 𝑻(𝑼𝑡,𝑒 𝑖(𝑥𝑖 𝑗))), (9) where 𝑥𝑖 𝑗and 𝜔𝑖 𝑗, 𝑗=1,…,𝑘are the quadrature nodes in 𝑉𝑖and the corresponding weights, respectively, of the considered quadrature formula. Now, it is easy to prove the following result: Theorem 1. Let 𝑼𝑒(𝑥)be a continuous stationary solution of (4), and the corresponding set of cell averages {𝑼𝑒 𝑖}𝑖∈defined as in (5). Assume that the corresponding selected in-cell steady states verify either one of the properties (i) 𝑼𝑡,𝑒 𝑖(𝑥)=𝑼𝑒(𝑥), for all 𝑥∈𝑉𝑖and for all 𝑖∈, (ii) 1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 𝑼𝑡,𝑒 𝑖(𝑥)𝑑𝑥 =𝑼𝑒 𝑖+(Δ𝑥𝑟), for all 𝑥∈𝑉𝑖and for all 𝑖∈and the global solution defined piece-wisely by 𝑼𝑡,𝑒 𝑖(𝑥)in 𝑉𝑖is continuous. Suppose that the reconstruction operator 𝑷𝑡 𝑖(𝑥)is well-balanced for {𝑼𝑡,𝑒 𝑖}𝑖∈in the sense 𝑷𝑡 𝑖(𝑥;{𝑼𝑒 𝑗}𝑗∈𝑖)=𝑼𝑡,𝑒 𝑖(𝑥), ∀𝑥∈𝑉𝑖, ∀𝑖∈. Then, the numerical scheme (9)is either exactly well-balanced (if (i) is satisfied) or well-balanced with order 𝑟(if (ii) is satisfied) for 𝑼𝑒(𝑥). Applied Numerical Mathematics 215 (2025) 138–156 143 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. Proof. Thanks to the continuity of the solution 𝑼𝑡,𝑒(𝑥), we automatically get 𝑖+1∕2(𝑡)=𝑭(𝑼𝑡,𝑒 𝑖(𝑥𝑖+1∕2)) and all the terms on the right hand side of (9)cancel, so the result follows. □ 3.2. Continuous stationary solution and reconstruction operator In this section we will describe an algorithm that allows to define the in-cell steady states 𝑼𝑡,𝑒 𝑖(𝑥)for lake-at-rest equilibrium (2) (either exactly or approximated). We will define as well well-balanced second-order reconstruction operator 𝑷𝑡 𝑖so that the resulting scheme (9)is well-balanced with order 2. For the sake of simplicity, we shall describe first in Subsection 3.2.1 the technique for the second-order scheme, although it can be extended to a high-order version as it is shown in Subsection 3.2.2. Given a collection of cell-averages {𝑼𝑖(𝑡)}, we shall describe first the construction of the in-cell steady states {𝑼𝑡,𝑒 𝑖(𝑥)} and then the reconstruction operator 𝑷𝑡 𝑖(𝑥). In what follows, we will not write explicitly the dependence on time in order to make the notation less cumbersome. 3.2.1. Second order method We start by describing the stationary solution and the reconstruction procedure in the case of a second-order well-balanced scheme. In-cell steady states. Assume a given function 𝑼(𝑥)and the corresponding average values {𝑼𝑖}computed as in (5)either exactly or using a high order quadrature formula. We then define the in-cell steady state functions 𝑼𝑒 𝑖(𝑥)=(𝜂𝑒 𝑖(𝑥),0) as follows: 1. Consider first  𝑷𝑖(𝑥)=(𝑝ℎ 𝑖,𝑝𝑞 𝑖)𝑇with 𝑝ℎ 𝑖,𝑝𝑞 𝑖∈ℙ1[𝑥]the usual MUSCL reconstruction operator with the minmod slope limiter. As it is usual in shallow water framework, we use the reconstruction of the free surface values 𝜂𝑖=ℎ𝑖+𝑏𝑖to define the height reconstruction, that is, we define: 𝑝𝑤 𝑖(𝑥)=𝑤𝑖+(𝛿𝑤)𝑖 Δ𝑥 (𝑥−𝑥𝑖), with (𝛿𝑤)𝑖=minmod(𝑤𝑖−𝑤𝑖−1,𝑤𝑖+1 −𝑤𝑖) for 𝑤∈{𝜂,𝑞}, and 𝑝ℎ 𝑖(𝑥)=𝑝𝜂 𝑖(𝑥)−𝑏(𝑥). Remark that strictly speaking 𝑝ℎ 𝑖(𝑥)is not a polynomial in general, but this does not affect to what is said afterwards, since everything is defined in terms of (𝛿𝜂)𝑖. 2. Define the operator 𝑟𝜂 𝑖∈ℙ1[𝑥]by 𝑟𝜂 𝑖(𝑥)=𝜂𝑖+𝛼𝑖(𝑝𝜂 𝑖(𝑥)−𝜂𝑖)with 𝛼𝑖=⎧ ⎪ ⎨ ⎪ ⎩ 1if |(𝛿𝜂)𝑖|≤Δ𝑥𝜇𝑠, 𝜇𝑠Δ𝑥 (𝛿𝜂)𝑖 otherwise. Notice that 𝑟𝜂 𝑖(𝑥)is a linear piecewise function with slope less or equal to the tangent of the angle of repose of the material 𝜇𝑠. Moreover, 𝑟𝜂 𝑖satisfies the conservation property 1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 𝑟𝜂 𝑖(𝑥)𝑑𝑥 =𝜂𝑖. Remark that in the case that  𝑷𝑖(𝑥)defines a lake-at-rest stationary state in the cell 𝑉𝑖, then 𝑟𝜂 𝑖(𝑥)=𝑝𝜂 𝑖(𝑥). However, one may find in general discontinuities at the interfaces 𝑥𝑖+1∕2, i.e., the piecewise function 𝑟𝜂(𝑥)given by 𝑟𝜂 𝑖(𝑥)in 𝑉𝑖is not a global stationary solution. 3. Consider now at each interface 𝑥𝑖+1∕2 the left (𝜂− 𝑖+1∕2) and right (𝜂+ 𝑖+1∕2) limits, as well as the mean inter-cell value (𝜂𝑖+1∕2) given by 𝜂− 𝑖+1∕2 =𝑟𝜂 𝑖(𝑥𝑖+1∕2), 𝜂+ 𝑖+1∕2 =𝑟𝜂 𝑖+1(𝑥𝑖+1∕2), 𝜂𝑖+1∕2 =1 2(𝜂− 𝑖+1∕2 +𝜂+ 𝑖+1∕2),(10) and consider 𝑟𝜂 𝑖(𝑥)∈ℙ2[𝑥]the quadratic polynomial verifying 𝑟𝜂 𝑖(𝑥𝑖)=𝜂𝑖, and 𝑟𝜂 𝑖(𝑥𝑖±1∕2)=𝜂𝑖±1∕2. Remark that this quadratic polynomial satisfies the conservation property up to second order: 1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 𝑟𝜂 𝑖(𝑥)𝑑𝑥 =𝜂𝑖+(Δ𝑥2). Then, consider the maximum slope in the volume 𝑀𝑖=max 𝑥∈𝑉𝑖|||(𝑟𝜂 𝑖)′(𝑥)|||and set 𝜂𝑒 𝑖(𝑥)={𝑟𝜂 𝑖(𝑥)if 𝑀𝑖≤𝜇𝑠, 𝑟𝜂 𝑖(𝑥)otherwise. Applied Numerical Mathematics 215 (2025) 138–156 144 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. Fig. 1. In-cell stationary steady states for the function −𝜇𝑠𝑥2∕2 + 1. Left figure shows the linear limited slope function 𝑟𝑒 𝑖and right figure shows the final steady state 𝜂𝑒 𝑖, including a zoom. Fig. 2. In-cell stationary steady states for the function 𝑥4. Left figure shows the linear limited slope function 𝑟𝑒 𝑖and right figure shows the final steady state 𝜂𝑒 𝑖. In blue cells where 𝜂𝑒 𝑖=𝑟𝑒 𝑖and in red where 𝜂𝑒 𝑖=𝑟𝑒 𝑖. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) Suppose that the internal slopes within the cells, 𝑀𝑖, are under the limiting tangent of the repose angle 𝜇𝑠. Then, the piece-wise function given by 𝜂𝑒 𝑖(𝑥)in 𝑉𝑖is a stationary solution restricted to the interior of each cell 𝑉𝑖and continuous at the interfaces, therefore it is a global continuous steady state. To illustrate the construction of the in-cell steady state, we show in Fig. 1the case of the initial function 𝜂(𝑥, 0) = −𝜇𝑠𝑥2∕2 + 1 in the interval [−1,1]. In this case, since 𝑀𝑖≤1, 𝜂𝑒 𝑖corresponds to the parabola 𝑟𝜂 𝑖. Remark that in this case, the global function defined piece-wise as 𝜂𝑒 𝑖inside each cell coincides (up to second order) with the initial function at discrete level, that is, they cell-averaged values coincide. In a similar way, Fig. 2shows the case for initial function 𝜂(𝑥, 0) = 𝑥4in the interval [0,1.2]. In this case, the slope on the right part of the function is greater than 𝜇𝑠, which means that we select the limited slope linear function 𝑟𝑖, while on the left part 𝜂𝑒 𝑖coincides with the parabola 𝑟𝜂 𝑖. Reconstruction operator. Using the same notation defined for the in-cell steady states, the reconstruction operator 𝑷𝑖(𝑥)inside the volume 𝑉𝑖will be set as: 𝑷𝑖(𝑥)={(𝜂𝑒 𝑖(𝑥)−𝑏(𝑥), 𝑝𝑞 𝑖(𝑥))𝑇,if 𝛼𝑖=1, (𝑝𝜂 𝑖(𝑥)−𝑏(𝑥), 𝑝𝑞 𝑖(𝑥))𝑇,otherwise. Notice that the step 3 in previous algorithm characterize the solutions we are going to preserve: those stationary solution whose discrete representation by a quadratic polynomial 𝑟𝜂(𝑥)are stationary (||(𝑟𝜂)′(𝑥)||<𝜇 𝑠). In particular, it holds for linear piecewise solutions whose slopes are lower than 𝜇𝑠in each volume, provided that we are able to have a exact approximation of (𝛿𝜂)𝑖. Applied Numerical Mathematics 215 (2025) 138–156 145 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. Theorem 2. Let 𝑼(𝑥)=(ℎ(𝑥),0)𝑇such that 𝜂(𝑥)=ℎ(𝑥)+𝑏(𝑥)is sufficiently smooth and satisfies |𝜂′(𝑥)|≤𝑀<𝜇 𝑠. Then the scheme (9) is well-balanced for 𝑼(𝑥). Proof. We will make use of Theorem 1. First, remark that for Δ𝑥sufficiently small, we may assume |(𝑝𝜂 𝑖)′(𝑥)|=|(𝛿𝜂)𝑖| Δ𝑥 <𝜇 𝑠. Therefore we have 𝑟𝜂 𝑖(𝑥)=𝑝𝜂 𝑖(𝑥). Moreover, assuming the function sufficiently smooth, we will have 𝜂𝑖+1∕2 =𝜂(𝑥𝑖+1∕2)+(Δ𝑥2). The quadratic polynomial 𝑟𝜂 𝑖(𝑥)can be written as 𝑟𝜂 𝑖(𝑥)=𝜂𝑖−1∕2 +2 Δ𝑥(𝜂𝑖−𝜂𝑖−1∕2)(𝑥−𝑥𝑖−1∕2)+2 Δ𝑥2(𝜂𝑖+1∕2 −2𝜂𝑖+𝜂𝑖−1∕2)(𝑥−𝑥𝑖−1∕2)(𝑥−𝑥𝑖). Since (𝑟𝜂 𝑖)′(𝑥)is a linear polynomial, the maximum 𝑀𝑖is reached at one of the interfaces of the volume 𝑉𝑖: (𝑟𝜂 𝑖)′(𝑥𝑖+1∕2)= 1 Δ𝑥(𝜂𝑖−1∕2 −𝜂𝑖+3(𝜂𝑖+1∕2 −𝜂𝑖))=𝜂′(𝑥𝑖)+(Δ𝑥) (𝑟𝜂 𝑖)′(𝑥𝑖−1∕2)= 1 Δ𝑥(3(𝜂𝑖−𝜂𝑖−1∕2)+(𝜂𝑖−𝜂𝑖+1∕2))=𝜂′(𝑥𝑖)+(Δ𝑥). Then for Δ𝑥sufficiently small we get 𝑀𝑖<𝜇 𝑠and we have that the in-cell steady state corresponds to the parabola, 𝜂𝑒 𝑖(𝑥)=𝑟𝜂 𝑖(𝑥), which is an approximation of 𝜂(𝑥). Moreover, 𝑷𝑖(𝑥)=(𝜂𝑒 𝑖(𝑥)−𝑏(𝑥),0)𝑇and the result follows. □ Remark 2. Although we have proved that the scheme (9)is in general well-balanced, we may have the exactly well-balanced property in some situations. For instance, if 𝑼(𝑥)=(ℎ(𝑥),0)𝑇is such that ℎ(𝑥)+𝑏(𝑥)corresponds to a linear polynomial with slope smaller than 𝜇𝑠, then the scheme is exactly well-balanced. Indeed, in such cases the slope approximations (𝛿𝜂)𝑖correspond exactly to 𝜂′(𝑥) and it is easy to check that 𝜂𝑖(𝑥)=𝜂(𝑥)=𝑝𝜂 𝑖(𝑥). Of special interest is the case of two connected linear functions (or a continuous polygonal in general) with slopes smaller than 𝜇𝑠. This is studied in Appendix A. Another particular case is when the initial condition produces a sequence of cell averages {𝜂𝑖}𝑖∈such that the reconstruction procedure results in 𝑟𝜂 𝑖inside the cells and they connect continuously at the interface. In that case, at discrete level, both the initial condition and the reconstructed function coincide. This initial condition is therefore preserved. 3.2.2. Generalization to higher order The algorithm introduced previously, can be adapted to obtain high order well-balanced schemes as follows: In-cell steady states. 1. We consider now  𝑷𝑖(𝑥)=(𝑝𝜂 𝑖−𝑏(𝑥),𝑝𝑞 𝑖)𝑇with 𝑝ℎ 𝑖,𝑝𝑞 𝑖∈ℙ𝑠[𝑥]a high order polynomial reconstruction operator (WENO, CWENO, etc.). Again, we rather use the reconstructions on the free surface 𝑝𝜂 𝑖. 2. Define the operator 𝑟𝜂 𝑖∈ℙ𝑠[𝑥]by 𝑟𝜂 𝑖(𝑥)=𝜂𝑖+𝛼𝑖(𝑝𝜂 𝑖(𝑥)−𝜂𝑖)with 𝛼𝑖={1if |(𝑝𝜂 𝑖)′(𝑥)|≤𝜇𝑠, 𝜇𝑠∕𝑀otherwise, with 𝑀=max 𝑥∈𝑉𝑖|(𝑝𝜂 𝑖)′(𝑥)|. Notice that, since 𝑝𝜂 𝑖(𝑥)is assumed to satisfy the conservation property, so does 𝑟𝜂 𝑖(𝑥): 1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 𝑟𝜂 𝑖(𝑥)𝑑𝑥 =𝜂𝑖. 3. Consider now at each interface 𝜂± 𝑖+1∕2,𝜂𝑖+1∕2 as in (10), and a Gauss quadrature formula of degree 𝑟≥𝑠in 𝑉𝑖with nodes 𝑥𝑖 𝑗∈ (𝑥𝑖−1∕2,𝑥𝑖+1∕2)and weights 𝜔𝑖 𝑗, 𝑗=1,…,𝑘. Consider 𝑟𝜂 𝑖(𝑥)∈ℙ𝑘+1[𝑥]the polynomial verifying 𝑟𝜂 𝑖(𝑥𝑖 𝑗)=𝑝𝜂 𝑖(𝑥𝑖 𝑗), 𝑗=1,…,𝑘, and 𝑟𝜂 𝑖(𝑥𝑖±1∕2)=𝜂𝑖±1∕2. Remark that this polynomial satisfies the conservation property up to order 𝑟≥𝑠: 1 Δ𝑥 𝑥𝑖+1∕2 ∫ 𝑥𝑖−1∕2 𝑟𝜂 𝑖(𝑥)𝑑𝑥 = 𝑘 ∑ 𝑗=1 𝜔𝑖 𝑗𝑟𝜂 𝑖(𝑥𝑖 𝑗)+(Δ𝑥𝑟)= 𝑘 ∑ 𝑗=1 𝜔𝑖 𝑗𝑝𝜂 𝑖(𝑥𝑖 𝑗)+(Δ𝑥𝑟)=𝜂𝑖+(Δ𝑥𝑟). Applied Numerical Mathematics 215 (2025) 138–156 146 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. Then, consider the maximum slope in the volume 𝑀𝑖=max 𝑥∈𝑉𝑖|||(𝑟𝜂 𝑖)′(𝑥)|||and set 𝜂𝑒 𝑖(𝑥)={𝑟𝜂 𝑖(𝑥)if 𝑀𝑖≤𝜇𝑠, 𝑟𝜂 𝑖(𝑥)otherwise. Reconstruction operator. Similarly as in the second-order case, the reconstruction operator 𝑷𝑖(𝑥)inside the volume 𝑉𝑖will be set as: 𝑷𝑖(𝑥)={(𝜂𝑒 𝑖(𝑥)−𝑏(𝑥), 𝑝𝑞 𝑖(𝑥))𝑇,if 𝛼𝑖=1, (𝑝𝜂 𝑖(𝑥)−𝑏(𝑥), 𝑝𝑞 𝑖(𝑥))𝑇,otherwise. Then the well-balance property of the reconstruction operator follows analogously. Theorem 3. Let 𝑼(𝑥)=(ℎ(𝑥),0)𝑇such that 𝜂(𝑥)=ℎ(𝑥)+𝑏(𝑥)is sufficiently smooth and satisfies |𝜂′(𝑥)|≤𝑀<𝜇 𝑠. Assume a reconstruction operator 𝑝𝜂 𝑖(𝑥)such that 𝜂(𝑥)=𝑝𝜂 𝑖(𝑥)+(Δ𝑥𝑝), 𝜂′(𝑥)=(𝑝𝜂 𝑖)′(𝑥)+(Δ𝑥𝑚), ∀𝑥∈𝑉𝑖 for some integers 1≤𝑚<𝑝. Then the scheme (9)is well-balanced for 𝑼(𝑥). Proof. We will make use of Theorem 1. First, remark that for Δ𝑥sufficiently small, we may assume |(𝑝𝜂 𝑖)′(𝑥)|≤|𝜂′(𝑥)|+(Δ𝑥𝑚)<𝜇 𝑠. Therefore, we have 𝑟𝜂 𝑖(𝑥)=𝑝𝜂 𝑖(𝑥). Moreover, assuming the function sufficiently smooth, we will have 𝜂𝑖+1∕2 =𝜂(𝑥𝑖+1∕2)+(Δ𝑥𝑝). Now, consider 𝑟𝜂 𝑖(𝑥)∈ℙ𝑘+1[𝑥]the polynomial verifying 𝑟𝜂 𝑖(𝑥𝑖 𝑗)=𝑝𝜂 𝑖(𝑥𝑖 𝑗), 𝑗=1,…,𝑘, and 𝑟𝜂 𝑖(𝑥𝑖±1∕2)=𝜂𝑖±1∕2, and we shall assume 𝑥𝑖−1∕2 <𝑥𝑖 1<…<𝑥𝑖 𝑘<𝑥 𝑖+1∕2. Let us use the notation 𝑥𝑖 0=𝑥𝑖−1∕2, 𝑥𝑖 𝑘+1 =𝑥𝑖+1∕2. Now, let 𝑥∈𝑉𝑖and consider 0≤𝑗≤𝑘such that 𝑥∈[𝑥𝑖 𝑗,𝑥𝑖 𝑗+1]. Then, we have (𝑟𝜂 𝑖)′(𝑥)= 𝑟𝜂 𝑖(𝑥𝑖 𝑗+1)−𝑟𝜂 𝑖(𝑥𝑖 𝑗) 𝑥𝑖 𝑗+1 −𝑥𝑖 𝑗 +(Δ𝑥)= 𝜂(𝑥𝑖 𝑗+1)−𝜂(𝑥𝑖 𝑗)+(Δ𝑥𝑝) 𝑥𝑖 𝑗+1 −𝑥𝑖 𝑗 +(Δ𝑥)=𝜂′(𝑥)+(Δ𝑥). Therefore, for Δ𝑥sufficiently small we get 𝑀𝑖<𝜇 𝑠and we have that the in-cell steady state corresponds to the polynomial, 𝜂𝑒 𝑖(𝑥)= 𝑟𝜂 𝑖(𝑥), which is an approximation of 𝜂(𝑥). Moreover, 𝑷𝑖(𝑥)=(𝜂𝑒 𝑖(𝑥)−𝑏(𝑥),0)𝑇and the result follows. □ 3.3. Time discretization We detail in this subsection the time discretization of the semi-discrete system (9). Let us recall that the main difficulty here is the definition of the friction term. It is usual, when dealing with friction terms (Coulomb, Manning, Darcy,...), to use a semi-implicit approach. In particular, for the case of single-layer models for granular flows, the technique introduced in [21] can be used. The basis idea is to assume that the friction force is opposed to the motion of the flow, but it cannot change the sign of the velocity vector (direction of movement). Moreover, if the magnitude of the friction force is greater than the rest of forces, then the flow must stop. For the sake of clarity, let us summarize this technique for a very simple case. Consider the ODE 𝑑𝑢 𝑑𝑡 =−𝑔𝜇ℎ sgn(𝑢),(11) where sgn(⋅)is the sign function and 𝜇>0a friction coefficient. Then, applying the semi-implicit method with a linearization of the pressure term, we get 𝑢𝑛+1 =𝑢𝑛−Δ𝑡𝑔𝜇ℎ𝑛sgn(𝑢𝑛+1), where the superscripts 𝑛,𝑛 +1denote the solutions at times 𝑡𝑛and 𝑡𝑛+1 =𝑡𝑛+Δ𝑡. Since the friction term cannot change the sign of the velocity, then sgn(𝑢𝑛+1)=sgn(𝑢𝑛)is assumed, and the following correction step on the update of 𝑢is performed: 𝑢𝑛+1 ={𝑢𝑛−Δ𝑡𝑔𝜇ℎ𝑛sgn(𝑢𝑛)if |𝑢𝑛|>Δ𝑡𝑔ℎ𝑛𝜇 0otherwise. (12) Thus, looking at previous equation, we observe that in practice (11)is explicitly discretized, except for the correction step, which sets 𝑢𝑛+1 =0when needed. This explicit treatment of the velocity is an advantage when compared to other techniques such as a regularization of the sign function. In particular, notice that null velocity is obtained exactly when the friction surpasses a critical level, reproducing properly the case of material at rest and stopping criterion. Applied Numerical Mathematics 215 (2025) 138–156 153 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. Fig. A.9. Sketch of two connected linear steady states with slopes 𝛽𝐿, 𝛽𝑅. 5. Conclusions We have introduced a well-balanced reconstruction procedure for shallow water granular flows model with Coulomb-type friction term, which allows to preserve lake-at-rest steady states. This procedure, together with an appropriate time discretization, is immediately generalized to an arbitrary order of accuracy. The resulting scheme is well-balanced in general for such steady states, although it is exactly well-balanced for some particular cases. Concretely, those steady states whose discrete representation, defined in the form of cell averaged values, corresponds to a global continuous stationary solution described inside the cells by the stationary reconstruction procedure described in Subsection 3.2 are exactly preserved. Remark that the work presented here is easily extended to the case of tilted coordinates and weakly non-hydrostatic framework (see [14]). All this has been shown in the numerical tests section. In particular, it is worth noticing that when the scheme is not exactly well-balanced for some initial steady state, the material slightly flows for a very short time and then an approximation of the initial condition is obtained, which is then preserved. That is, the scheme is well-balanced. Moreover, it has been shown that the expected order of accuracy is obtained for the second and third-order schemes. Let us remark that, up to our knowledge, this is the first time that a high-order well-balanced finite volume method has been proposed for such granular flow model. Future works could study the extension of this technique to systems including the vertical discretization of the granular system, where the vertical effects are essential to obtain realistic simulation in complex viscoplastic flows. CRediT authorship contribution statement M.J. Castro Díaz: Writing – review & editing, Software, Methodology, Funding acquisition, Conceptualization. C. Escalante: Writing – review & editing, Writing – original draft, Software, Methodology, Conceptualization. J. Garres-Díaz: Writing – review & editing, Writing – original draft, Software, Methodology, Conceptualization. T. Morales de Luna: Writing – review & editing, Software, Methodology, Funding acquisition, Conceptualization. Acknowledgements This work is supported by projects PID2022-137637NB-C21 and PID2022-137637NB-C22 funded by MICIU/AEI/10.13039/ 501100011033/ and FEDER, UE, and project PDC2022-133663-C21 funded by MICIU/AEI/10.13039/501100011033/ and European Union NextGenerationEU/PRTR. The work has also been partially supported by project PROYEXCEL_00525 funded by Junta de Andalucía. Appendix A. The case of connected stationary linear states In this Appendix we consider first the case of two connected linear polynomials and then the more general case of a polygonal, as they may be of particular interest. Let 𝛽𝐿,𝛽𝑅,𝛾 ∈ℝ, with |𝛽𝐿|≤𝜇𝑠, |𝛽𝑅|≤𝜇𝑠and consider 𝜂(𝑥)={𝛽𝐿(𝑥−𝑐)+𝛾, for 𝑥≤𝑐, 𝛽𝑅(𝑥−𝑐)+𝛾, for 𝑥>𝑐. (A.1) We will consider a space discretization 𝑉𝑖=[𝑥𝑖−1∕2,𝑥𝑖+1∕2]such that 𝑥𝑗+1∕2 =𝑐for some index 𝑗(see Fig. A.9). Moreover, we shall identify the cell averages with the values at the center of the cells, since they are second-order approximations. Remark that the constructions described in Subsection 3.2 are invariant by vertical translation in the sense that if we add a constant 𝛾to the cell averages, then the reconstructions coincide with those for the original values translated by the constant 𝛾. Therefore, we Applied Numerical Mathematics 215 (2025) 138–156 154 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. assume without loss of generality 𝛾=0. Moreover, we suppose |𝛽𝐿|>|𝛽𝑅|, otherwise we may consider a symmetry with respect to the vertical axis 𝑥=𝑥𝑗+1∕2 and everything follows analogously. Remark that the reconstructed slopes (𝛿𝜂)𝑖∕Δ𝑥correspond exactly to 𝛽𝐿for 𝑖≤𝑗−1and to 𝛽𝑅for 𝑖≥𝑗+2. The only particular cases are then the cells 𝑗and 𝑗+1, where we get (𝛿𝜂)𝑗 Δ𝑥 =minmod(𝛽𝐿+𝛽𝑅 2 ,𝛽𝐿), (𝛿𝜂)𝑗+1 Δ𝑥 =minmod(𝛽𝐿+𝛽𝑅 2 ,𝛽𝑅). A systematic study of the different possibilities when |𝛽𝐿|>|𝛽𝑅|let us show that minmod(𝛽𝐿+𝛽𝑅 2 ,𝛽𝐿)=𝛽𝐿+𝛽𝑅 2 , minmod(𝛽𝐿+𝛽𝑅 2 ,𝛽𝑅)={𝛽𝑅,if 𝛽𝐿𝛽𝑅≥0, 0,otherwise. It is clear then that all slopes are smaller than the critical value 𝜇𝑠and 𝛼𝑖=1. Now, we easily check that at interface 𝑥𝑗+1∕2 we have 𝜂− 𝑗+1∕2 =Δ𝑥 4 (𝛽𝑅−𝛽𝐿), 𝜂+ 𝑗+1∕2 =⎧ ⎪ ⎨ ⎪ ⎩ 0,if 𝛽𝐿𝛽𝑅≥0, Δ𝑥 2 𝛽𝑅,otherwise, 𝜂𝑗+1∕2 =⎧ ⎪ ⎨ ⎪ ⎩ Δ𝑥 8 (𝛽𝑅−𝛽𝐿),if 𝛽𝐿𝛽𝑅≥0, Δ𝑥 8 (3𝛽𝑅−𝛽𝐿),otherwise, at the left interface, 𝑥𝑗−1∕2, 𝜂− 𝑗−1∕2 =−Δ𝑥𝛽𝐿, 𝜂+ 𝑗−1∕2 =−Δ𝑥(3 4𝛽𝐿+1 4𝛽𝑅), 𝜂𝑗−1∕2 =−Δ𝑥(7 8𝛽𝐿+1 8𝛽𝑅), and at the right, 𝑥𝑗+3∕2, 𝜂− 𝑗+3∕2 =⎧ ⎪ ⎨ ⎪ ⎩ Δ𝑥𝛽𝑅,if 𝛽𝐿𝛽𝑅≥0, Δ𝑥 2 𝛽𝑅,otherwise, 𝜂+ 𝑗+3∕2 =Δ𝑥𝛽𝑅, 𝜂𝑗+3∕2 =⎧ ⎪ ⎨ ⎪ ⎩ Δ𝑥𝛽𝑅,if 𝛽𝐿𝛽𝑅≥0, 3Δ𝑥 4 𝛽𝑅,otherwise. For any cell 𝑉𝑖with 𝑖≠𝑗−1,𝑗,𝑗+1,𝑗+2, the second-order polynomial 𝑟𝑖(𝑥)coincides with the original linear function inside the cell and we have 𝜂𝑖(𝑥)=𝑝𝜂 𝑖(𝑥)=𝜂(𝑥)for 𝑥∈𝑉𝑖. Now, in order to consider the particular case 𝑘∈{𝑗−1,𝑗,𝑗+1,𝑗+2}, let us write the general polynomial 𝑟𝑘(𝑥)∈ℙ2[2] in the form: 𝑟𝑘(𝑥)=𝜂𝑘−1∕2 +2𝜂𝑘−𝜂𝑘−1∕2 Δ𝑥 (𝑥−𝑥𝑘−1∕2)+2𝜂𝑘+1∕2 −2𝜂𝑘+𝜂𝑘−1∕2 Δ𝑥2(𝑥−𝑥𝑘−1∕2)(𝑥−𝑥𝑘), which gives 𝑟′ 𝑘(𝑥)= 𝜂𝑘+1∕2 −𝜂𝑘−1∕2 Δ𝑥 +4𝜂𝑘+1∕2 −2𝜂𝑘+𝜂𝑘−1∕2 Δ𝑥2(𝑥−𝑥𝑘). Since 𝑟′ 𝑘(𝑥)is linear, we get max 𝑥∈𝑉𝑘|𝑟′ 𝑘(𝑥)|=max{|𝑟′ 𝑘(𝑥𝑘−1∕2)|,|𝑟′ 𝑘(𝑥𝑘+1∕2)|}= 1 Δ𝑥max {|𝜂𝑘−𝜂𝑘+1∕2 +3(𝜂𝑘−𝜂𝑘−1∕2)|,|3(𝜂𝑘+1∕2 −𝜂𝑘)+𝜂𝑘−1∕2 −𝜂𝑘|}. After some easy calculations we have: •In the cell 𝑉𝑗−1 two possibilities arise: –If |11𝛽𝐿−3𝛽𝑅|≤8𝜇𝑠, then 𝑀𝑗−1 ≤𝜇𝑠and we set 𝜂𝑒 𝑗−1(𝑥)=𝑟𝑗−1(𝑥). –If |11𝛽𝐿−3𝛽𝑅|>8𝜇𝑠, then 𝑀𝑗−1 >𝜇 𝑠and we set 𝜂𝑒 𝑗−1(𝑥)=𝑟𝑗−1(𝑥)=𝑝𝜂 𝑗−1(𝑥). •In the cell 𝑉𝑗we have 𝑀𝑗≤𝜇𝑠and 𝜂𝑒 𝑗(𝑥)=𝑟𝑗(𝑥). Moreover, if 𝛽𝐿𝛽𝑅≥0, then 𝑟𝑗(𝑥)is linear. •In the cell 𝑉𝑗+1 we have 𝑀𝑗+1 ≤𝜇𝑠and 𝜂𝑒 𝑗+1(𝑥)=𝑟𝑗+1(𝑥). •In the cell 𝑉𝑗+2 two possibilities arise: –If 𝛽𝐿𝛽𝑅≥0, then 𝑟𝑗+2(𝑥)is linear and coincides with the original function inside the cell: 𝜂𝑒 𝑗+2(𝑥)=𝑟𝑗+2(𝑥)=𝜂(𝑥). –If 𝛽𝐿𝛽𝑅<0, then we only get 𝑀𝑗+2 <𝜇 𝑠in the case |𝛽𝑅|≤4 7𝜇𝑠. 𝜂𝑒 𝑗+2(𝑥)=⎧ ⎪ ⎨ ⎪ ⎩ 𝑟𝑗+2(𝑥),if 𝛽𝐿𝛽𝑅<0and |𝛽𝑅|≤4 7𝜇𝑠, 𝑟𝑗+2(𝑥)=𝑝𝜂 𝑗+2(𝑥),if 𝛽𝐿𝛽𝑅<0and |𝛽𝑅|>4 7𝜇𝑠. Applied Numerical Mathematics 215 (2025) 138–156 155 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. Therefore, we conclude that when 𝛽𝐿and 𝛽𝑅are such that {|11𝛽𝐿−3𝛽𝑅|≤8𝜇𝑠}and {𝛽𝐿𝛽𝑅≥0or |𝛽𝑅|≤4 7𝜇𝑠},(A.2) the reconstruction procedure of 𝜂(𝑥)is such that it is (exactly) preserved by the scheme. If the previous conditions are not satisfied, then we potentially get two cells 𝑉𝑗−1 and 𝑉𝑗+2 where the function is not necessarily preserved. We recall that we have assumed the case ||𝛽𝐿||>||𝛽𝑅||. Otherwise, by symmetry reasons one obtains {|11𝛽𝑅−3𝛽𝐿|≤8𝜇𝑠}and {𝛽𝐿𝛽𝑅≥0or |𝛽𝐿|≤4 7𝜇𝑠}. Remark that this can be easily generalized. Consider now the general case, where the initial condition corresponds to a continuous piece-wise linear function with slopes under the critical value 𝜇𝑠. Concretely, let us consider the piece-wise linear function connecting nodes with abscissas 𝑋𝑗for 𝑗=1,..., 𝑀. Denoting by 𝛿=min 𝑗=2,...,𝑀(𝑋𝑗+1 −𝑋𝑗), if Δ𝑥<𝛿∕3, which means that there are at least three consecutive cells with the same slope, then we locally find two connected linear polynomials as depicted in Fig. A.9. Therefore, all that has been said previously remains valid. In particular, in the numerical simulations we observe that this initially function is usually preserved with the exception of the joining cells. In those cells, the initial condition is modified during the first iterations and then a discrete steady state is obtained which is an approximation of the original functions. That is, the scheme is well-balanced, but not exactly well-balanced. This behavior will be analyzed in the numerical tests section (see Subsection 4.2). References [1] P. Jop, Y. Forterre, O. Pouliquen, A constitutive law for dense granular flows, Nature 441 (2006) 727–730. [2] P.-Y. Lagrée, L. Staron, S. Popinet, The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes with a 𝜇(I)-rheology, J. Fluid Mech. 686 (2011) 378–408. [3] C. Lusso, A. Ern, F. Bouchut, A. Mangeney, M. Farin, O. Roche, Two-dimensional simulation by regularization of free surface viscoplastic flows with DruckerPrager yield stress and application to granular collapse, J. Comput. Phys. 333 (2017) 387–408. [4] I.R. Ionescu, A. Mangeney, F. Bouchut, R. Roche, Viscoplastic modeling of granular column collapse with pressure-dependent rheology, J. Non-Newton. Fluid Mech. 219 (2015) 1–18. [5] N. Martin, I.R. Ionescu, A. Mangeney, F. Bouchut, M. Farin, Continuum viscoplastic simulation of a granular column collapse on large slopes: 𝜇(I) rheology and lateral wall effects, Phys. Fluids 29 (2017) 013301. [6] O. Pouliquen, Y. Forterre, A non-local rheology for dense granular flows, Philos. Trans. R. Soc. A, Math. Phys. Eng. Sci. 367 (2009) 5091–5107. [7] F. Boyer, É. Guazzelli, O. Pouliquen, Unifying suspension and granular rheology, Phys. Rev. Lett. 107 (2011). [8] M. Bouzid, M. Trulsson, P. Claudin, E. Clément, B. Andreotti, Nonlocal rheology of granular flows across yield conditions, Phys. Rev. Lett. 111 (2013). [9] M. Rauter, The compressible granular collapse in a fluid as a continuum: validity of a Navier–Stokes model with 𝜇(I)-rheology, J. Fluid Mech. 915 (2021). [10] S.B. Savage, K. Hutter, The motion of a finite mass of granular material down a rough incline, J. Fluid Mech. 199 (1989) 177–215. [11] J.M.N.T. Gray, A.N. Edwards, A depth-averaged 𝜇(I)-rheology for shallow granular free-surface flows, J. Fluid Mech. 755 (2014) 503–534. [12] A.N. Edwards, J.M.N.T. Gray, Erosion-deposition waves in shallow granular free-surface flows, J. Fluid Mech. 762 (2015) 35–67. [13] E.D. Fernández-Nieto, J. Garres-Díaz, A. Mangeney, G. Narbona-Reina, A multilayer shallow model for dry granular flows with the 𝜇(𝐼)-rheology: application to granular collapse on erodible beds, J. Fluid Mech. 798 (2016) 643–681. [14] J. Garres-Díaz, E.D. Fernández-Nieto, A. Mangeney, T.M. de Luna, A weakly non-hydrostatic shallow model for dry granular flows, J. Sci. Comput. 86 (2021). [15] C. Escalante, E.D. Fernández-Nieto, J. Garres-Díaz, A. Mangeney, Multilayer shallow model for dry granular flows with a weakly non-hydrostatic pressure, J. Sci. Comput. 96 (2023). [16] F. Bouchut, J.M. Delgado-Sánchez, E.D. Fernández-Nieto, A. Mangeney, G. Narbona-Reina, A bed pressure correction of the friction term for depth-averaged granular flow models, Appl. Math. Model. 106 (2022) 627–658. [17] E. Audusse, M. Bristeau, B. Perthame, J. Sainte-Marie, A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation, ESAIM: Math. Model. Numer. Anal. 45 (2010) 169–200. [18] E.D. Fernández-Nieto, E.H. Koné, T.C. Rebollo, A multilayer method for the hydrostatic Navier-Stokes equations: a particular weak solution, J. Sci. Comput. 60 (2013) 408–437. [19] S. Boscarino, G. Russo, M. Semplice, High order finite volume schemes for balance laws with stiff relaxation, Comput. Fluids 169 (2018) 155–168. [20] A. Kurganov, J. Miller, Central-upwind scheme for Savage–Hutter type model of submarine landslides and generated tsunami waves, Comput. Methods Appl. Math. 14 (2014) 177–201. [21] A. Mangeney-Castelnau, J.-P. Vilotte, M.O. Bristeau, B. Perthame, F. Bouchut, C. Simeoni, S. Yerneni, Numerical modeling of avalanches based on Saint Venant equations using a kinetic scheme, J. Geophys. Res., Solid Earth 108 (2003), https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2002JB002024. [22] A. Bermudez, M.E. Vazquez, Upwind methods for hyperbolic conservation laws with source terms, Comput. Fluids 23 (1994) 1049–1071. [23] J.M. Greenberg, A.Y. Leroux, A well-balanced scheme for the numerical processing of source terms in hyperbolic equations, SIAM J. Numer. Anal. 33 (1996) 1–16. [24] L. Gosse, A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms, Comput. Math. Appl. 39 (2000) 135–159. [25] E. Audusse, F. Bouchut, M. Bristeau, R. Klein, B. Perthame, A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows, SIAM J. Sci. Comput. 25 (2004) 2050–2065. [26] M.J. Castro Díaz, J.A. López-García, C. Parés, High order exactly well-balanced numerical methods for shallow water systems, J. Comput. Phys. 246 (2013) 242–264. [27] C. Berthon, C. Chalons, A fully well-balanced, positive and entropy-satisfying Godunov-type method for the shallow-water equations, Math. Comput. 85 (2015) 1281–1307. Applied Numerical Mathematics 215 (2025) 138–156 156 M.J. Castro Díaz, C. Escalante, J. Garres-Díaz et al. [28] M.J. Castro, T. Morales de Luna, C. Parés, Well-balanced schemes and path-conservative numerical methods, in: R. Abgrall, C.-W. Shu (Eds.), Handbook of Numerical Analysis, in: Handbook of Numerical Methods for Hyperbolic Problems Applied and Modern Issues, vol. 18, Elsevier, 2017, pp. 131–175. [29] E.D. Fernández-Nieto, F. Bouchut, D. Bresch, M.J. Castro Díaz, A. Mangeney, A new Savage-Hutter type model for submarine avalanches and generated tsunami, J. Comput. Phys. 227 (2008) 7720–7754. [30] J. Zhai, W. Liu, L. Yuan, Solving two-phase shallow granular flow equations with a well-balanced noc scheme on multiple gpus, Comput. Fluids 134–135 (2016) 90–110. [31] E.D. Fernández-Nieto, J. Garres-Díaz, A. Mangeney, G. Narbona-Reina, 2D granular flows with the 𝜇(𝐼)rheology and side walls friction: a well-balanced multilayer discretization, J. Comput. Phys. 356 (2018) 192–219. [32] F. Campos, M. Sepúlveda, R. Abarca del Río, D. Issler, Estudio de modelos de avalancha usando esquemas de volúmenes finitos bien balanceados, Obras y Proyectos, 2023. [33] M. Sanz-Ramos, E. Bladé, P. Oller, G. Furdada, Numerical modelling of dense snow avalanches with a well-balanced scheme based on the 2D shallow water equations, J. Glaciol. (2023) 1–17. [34] E. Fernández-Nieto, J. Garres-Díaz, P. Vigneaux, Multilayer models for hydrostatic Herschel-Bulkley viscoplastic flows, Comput. Math. Appl. 139 (2023) 99–117. [35] F. Bouchut, Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws, Birkhäuser, Basel, 2004. [36] M.J. Castro, C. Parés, Well-balanced high-order finite volume methods for systems of balance laws, J. Sci. Comput. 82 (2020). [37] M.J. Castro Diaz, I. Gómez Bueno, C. Parés, Hyperbolic Problems: Theory, Numerics, Applications. Volume I, Springer, 2022. [38] C. Parés, Numerical methods for nonconservative hyperbolic systems: a theoretical framework, SIAM J. Numer. Anal. 44 (2006) 300–321. [39] M.J. Castro Díaz, E.D. Fernández-Nieto, A class of computationally fast first order finite volume solvers: PVM methods, SIAM J. Sci. Comput. 34 (2012) A2173–A2196. [40] I. Gómez-Bueno, M.J.C. Díaz, C. Parés, G. Russo, Collocation methods for high-order well-balanced methods for systems of balance laws, Mathematics 9 (2021) 1799. [41] C.A. Kennedy, M.H. Carpenter, Additive Runge–Kutta schemes for convection–diffusion–reaction equations, Appl. Numer. Math. 44 (2003) 139–181. [42] L. Pareschi, G. Russo, Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation, J. Sci. Comput. 25 (2005) 129–155. [43] S. Boscarino, F. Filbet, G. Russo, High order semi-implicit schemes for time dependent partial differential equations, J. Sci. Comput. 68 (2016) 975–1001. [44] S. Gottlieb, C.-W. Shu, Total variation diminishing Runge-Kutta schemes, Math. Comput. 67 (1996). [45] X. Zhang, C.-W. Shu, Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments, Proc. Royal Soc. A, Math. Phys. Eng. Sci. 467 (2011) 2752–2776. [46] I. Cravero, M. Semplice, G. Visconti, Optimal definition of the nonlinear weights in multidimensional central WENOZ reconstructions, SIAM J. Numer. Anal. 57 (2019) 2328–2358.