# Multilayer models for hydrostatic Herschel-Bulkley viscoplastic flows

Enrique Domingo Fernández-Nieto, J Garres-Díaz, Paul Vigneaux

► To cite this version:

Enrique Domingo Fernández-Nieto, J Garres-Díaz, Paul Vigneaux. Multilayer models for hydrostatic Herschel-Bulkley viscoplastic flows. Computers & Mathematics with Applications, 2023, 139, pp.99-117. 10.1016/j.camwa.2023.03.018. hal-03679369

HAL Id: hal-03679369

https://hal.science/hal-03679369

Submitted on 3 Apr 2023

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 models for hydrostatic
Herschel-Bulkley viscoplastic flows

E.D. Fernández-Nieto,\* J. Garres-Díaz,† P. Vigneaux‡§

Preprint / Authors' version as of: March 7, 2023. Now in:
Computers and Mathematics with Applications 139 (2023)
Pages 99-117. DOI: 10.1016/j.camwa.2023.03.018

**Abstract**

Starting from Navier-Stokes' equation we derive two shallow water multilayer models for yield stress fluids, depending on the asymptotic analysis. One of them takes into account the normal stress contributions, making possible to recover a pseudoplug layer instead of a purely plug zone. A specific numerical scheme is designed to solve this model thanks to a finite volume discretization. It involves well-balancing techniques to be able to compute accurately the transitions between yielded and unyielded (or pseudoplug) zones, an important feature of the original partial differential equations' model. We perform numerical simulations on various test cases relevant to this physics: analytical solution of a uniform flow, steady solutions for arrested state, and a viscoplastic dam break. Simulations agree well when we perform comparisons with physical experiments of the group of Christophe Ancey (EPFL) and we make a comparative study including shallow water models and lubrication models that they present in [Ancey et al. Viscoplastic dam break waves: Review of simple computational approaches and comparison with experiments. Advances in Water Resources 48 (2012) 79-91]. Thanks to the multilayer structure of our model, we can go further on the description of the vertical structure associated to the (bottom) sheared layer and the top (pseudo-)plug layer.

**Key words:** multilayer shallow water model; finite volume; well-balanced; lubrication theory; dam break flow; comparison with physical experiments

<footnote>*Dpto. Matemática Aplicada I, Universidad de Sevilla, Spain (edofer@us.es)</footnote>

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

<footnote>‡École Normale Supérieure de Lyon, UMPA CNRS UMR 5669, 46 allée d'Italie, 69364 Lyon Cedex 07, France</footnote>

<footnote>§LAMFA, CNRS UMR 7352 Université de Picardie Jules Verne, 33 rue Saint Leu, F-80039 Amiens Cedex
1, France. (paul.vigneaux@math.cnrs.fr)</footnote>
# 1 Introduction

Flows where the material can be either in a fluid state or in a solid state are ubiquitous in the nature or in the industry. One class of models to describe such materials is the yield stress fluids formalism, which goes back to the turn of the 20th century with Schwedoff [47], Bingham [13] and Herschel & Bulkley [37]. For these models, also called viscoplastic models, if the stress of the material is above the yield stress threshold $\tau_y$, it behaves like a fluid. While the material is rigid when its stress is below $\tau_y$. Such flow behavior can be encountered in many practical situations such as food pastes, cosmetics creams, heavy oils, mud and clays, lava flows and avalanches [10, 23, 33].

In this article, we are interested in the derivation of multilayer integrated Herschel-Bulkley models for shallow flows, as well as in the design of numerical algorithms to solve the resulting equations. The goal is to simulate the evolution of thin sheets of viscoplastic materials on inclined planes and, in particular, to be able to recover general velocity profiles and to compute precisely the transitions between fluid and rigid states. We have particularly in mind applications in geophysics.

As far as viscoplastic flows are concerned, the constitutive law for Herschel-Bulkley fluids is considered. It links the deviatoric stress tensor $\tau$ and the rate of deformation tensor

$$D(\mathbf{u}) = \frac{1}{2}(\nabla \mathbf{u} + (\nabla \mathbf{u})'), \quad (1)$$

by defining

$$\left\{ \begin{array}{ll} \tau = 2^n K \|D(\mathbf{u})\|^{n-1} D(\mathbf{u}) + \frac{\tau_y}{\|D(\mathbf{u})\|} D(\mathbf{u}) & \text{if } \|D(\mathbf{u})\| \neq 0, \\ \\ \|\tau\| \le \tau_y & \text{if } \|D(\mathbf{u})\| = 0, \end{array} \right. \quad (2)$$

where $\tau_y$ ($Pa$), $K$ ($Pa \cdot s^n$) and $n$ are the yield stress, consistency, and power index of the material, respectively. Notice that if $n = 1$ and $K$ is identified with $\eta$, the dynamic viscosity, we find the simplest and most emblematic constitutive law for viscoplastic fluids, the Bingham law. This latter law can be viewed as a generalization (adding a shift with the threshold $\tau_y$) of the Newtonian constitutive law leading to Navier-Stokes' equation, namely $\tau = 2\eta D(\mathbf{u})$. One of the difficulties of Bingham-type laws is that the deviatoric stress is not uniquely defined when the material is rigid ($\|D(\mathbf{u})\| = 0$). Mathematical tools belonging to non-smooth optimization need to be used in order to properly solve these models (e.g. duality methods like augmented Lagrangian). In particular this is crucial in order to capture a flow which evolves from a deformable state to a rigid state. Of note, nearly one century after Bingham, in the context of dense granular flows, the so called $\mu(I)$ rheology was introduced (see [21, 40]), which makes use of a variable coefficient for $\tau_y$ through the pressure and the state of the granular packing (encoded in the friction coefficient $\mu(I)$). In other words, this $\mu(I)$ constitutive law can be viewed has a yield criterion of the Drucker-Prager type [24]. We mention here Bingham and $\mu(I)$ together because they both share the threshold behavior and are both involved in some geophysical applications. However their physical background is different: while as said $\mu(I)$ is linked to granular material, Bingham is better suited to model more cohesive materials. To fix the ideas, one of the most used laboratory prototype of a viscoplastic material is the Carbopol mixture or the Kaolin mixture [14]. The differences in the equations' structure equally translate in the obtained vertical velocity profiles for both kind of rheologies: for Bingham, there is a strong shear zone at the bottom and a moving "plug"-like layer at the top of the material. While for the $\mu(I)$ rheology, there can be a significant static unsheared bottom layer with a sheared moving layer at the top. So in terms of (vertical) shear, these
two rheologies are totally opposed. This difference can be observed in real flows when the material is of the “more cohesive” type, as opposed to the “granular” type.

In geophysical flows, it is known that the full resolution of 3D Navier-Stokes free-surface flows is implementable but computationally very expensive. This is even worse for materials with more complex rheologies (going from Newtonian to non-Newtonian constitutive laws) because the models involve more non-linearities (power laws, threshold, time-space dependent laws, etc.) leading to supplementary algorithmic costs. As a consequence, computation times for 3D flows are not feasible in practice. Hence, a classical approach to lower the computational cost is to reduce the dimension of the problem from 3D to 2D through asymptotic analysis. Indeed, geophysical flows often verify that the characteristic horizontal length is much larger than the characteristic vertical height of the flowing material. This leads to the well-known shallow-water or Saint-Venant models, originally derived in the Newtonian context and subject to a vast amount of literature. The derivation of shallow-water models based on non-Newtonian constitutive laws is also very active. Some Saint-Venant type models for Herschel-Bulkley fluids are found in [41, 30, 1]. Lubrication models were introduced before the shallow-water models and were heavily used in practical simulations (see [44, 39, 8, 9]). We also refer the reader to the reviews [2, 48, 46]. The comparisons between these two approaches will be discussed further in the following of the article. Some other approaches to derive shallow-water models for other non-Newtonian fluids, namely granular flows, can be seen in [5, 45].

One of the key point in deriving these 2D reduced model is the description of the vertical velocity profiles. Indeed, it is often necessary to postulate a certain form (constant, power law, etc.) of this profile in order to perform the derivation of the asymptotic. One way of handling more general velocity profiles while keeping the numerical cost reasonable is to use the numerical multilayer approaches as initially derived in [6] for Navier-Stokes equations under the assumption of a hydrostatic pressure. In [29] a multilayer system is derived, whose solution is a particular weak solution of the Navier-Stokes system with a piecewise constant vertical profile of the horizontal velocity. Taking into account the normal jump conditions and the incompressibility condition the vertical velocity is deduced, being a piecewise linear profile. A generalization of the multilayer model with the $\mu(I)$-rheology was introduced in [27, 28]. This model was able to approximate changes in the typology of vertical profiles of the velocity, without prescribing it, as well as to approximate the static/flowing interface characterizing these flows. Similar results are expected when applying the multilayer approach to other viscoplastic fluids. In particular, it should be possible to approximate the yield surface characterizing Herschel-Bulkley fluids.

Another difficulty is the design of well-balanced numerical methods for depth-averaged and multilayer models. A well-balanced finite volume method combined with a duality technique is proposed in [1] for a Saint-Venant Herschel-Bulkley type model. The particularity of this model is that it is written as a variational inequality. The duality technique, namely Augmented Lagrangian or Bermúdez-Moreno methods, allows to rewrite the system in terms of an optimization problem with constraints. The Lagrange multiplier associated to the optimization problem is used in the finite volume step to obtain a well-balanced method, preserving non-trivial stationary solutions. The main drawback of this technique is that it is necessary to solve a fixed point problem at each time step. Then, an iterative algorithm is considered. A different alternative, although it is less accurate, is the use of a regularization of the stress tensor (see [36, 11, 43]). In [27] a well-balanced approximation of multilayer model is proposed, with a regularization of the stress tensor for the $\mu(I)$-rheology. Note that considering a multilayer model as a variational inequality would lead to solve a more complex optimization problem than for shallow one-layer models, then increasing the computational cost.

In this paper, we first focus on the deduction of a multilayer approach for Herschel-Bulkley fluids (see equation 3). Note that Bingham model can be seen as a particular case of Herschel-Bulkley model. A regularization of the stress tensor is also considered. Two
multilayer models are introduced, one of them including the normal stress contributions thanks to a particular asymptotic hypothesis in the pseudoplug layer. A well-balanced discretization of the multilayer model is then proposed. In the numerical tests a systematic comparison with Shallow Water type and lubrication models is presented. The ability of the model to capture different velocity profiles and the interface between sheared/(pseudo-)plug are also analyzed. Moreover, relevant results concerning the pseudoplug layer are recovered. For instance, the change of the convexity of the yield surface near the front when considering the normal stress contributions and parabolic vertical profile of the velocity there.

The paper is organized as follows. In the next section, the starting 3D governing equations for viscoplastic material are introduced. These equations have a structure of the Navier-Stokes' type but are extended to take into account the ability of the material to be either fluid or rigid, as sketched previously. In section 3, the derivation of the multilayer Herschel-Bulkley models (with and without normal stress contributions) and some relevant associated steady states are presented. Section 4 is dedicated to the design of a well-balanced finite volume discretization of previous model. Care is devoted to preserve aforementioned stationary states. Numerical tests are then presented in section 5. We begin by a study of the ability of the scheme to compute an analytical solution for a uniform flow and then to preserve a family of steady solutions without or with wet/dry fronts. Then, we revisit and extend, thanks to this scheme, the 1D viscoplastic dam break problem as presented by Ancey and colleagues in [4, 3]. This also leads to a comparative analysis with the results of lubrication models. Interesting outcomes are presented on the ability of the simulations to reproduce the sheared/(pseudo-)plug layers within the vertical.

## 2 Governing equations

For the sake of simplicity the two-dimensional case is only considered in the paper. Note that this could entail a limitation in the case of narrow channels with a strong side walls friction, or in real applications on 3D domains. The incompressible Navier-Stokes system describing the dynamics of a fluid with velocity $u \in \mathbb{R}^2$ and constant density $\rho \in \mathbb{R}$, together with an appropriate definition of the stress tensor accounting for the rheological behavior of the non-Newtonian fluid is considered. This system reads

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

where $g$ is the gravity force. The total stress tensor is

$$ \sigma = -p\mathcal{I} + \tau, $$

with $p \in \mathbb{R}$ the pressure, $\mathcal{I}$ the 2D identity tensor and $\tau$ the deviatoric stress tensor, which is defined by the rheology. For Herschel-Bulkley viscoplastic fluids the strain-rate and deviatoric tensors are defined by (1) and (2), respectively.

Focusing on the first case in (2), when $||D(u)|| \neq 0$, this definition leads to

$$ \tau = \eta D(u), \quad \text{with} \quad \eta = \frac{\tau_y}{||D(u)||} + 2^n K ||D(u)||^{n-1}, $$

being $\eta \in \mathbb{R}$ the generalized viscosity coefficient, which depends on the velocity. Notice that it is not defined if the strain rate vanishes. In this case, $\tau$ is a multivalued tensor. In
order to avoid this singularity, we consider a regularization technique (see [36, 43]), and we finally obtain

$$ \eta = \frac{\tau_y + 2^n K \|D(\mathbf{u})\|^n}{\sqrt{\|D(\mathbf{u})\|^2 + \delta_0^2}}, \quad (4) $$

where $\delta_0$ is the regularization parameter. Note that $\delta_0$ must be considered small enough, in such a way that the regularization error does not dominate over the space discretization error (see [42]). As mentioned in the introduction, a different alternative would be using a duality method such as Augmented Lagrangian [31] of Bermúdez-Moreno [12] methods. These methods were considered in [25, 26] for shallow Bingham fluids in the 1D and 2D cases. On the one hand, the main advantage of the regularization method, with respect to these duality methods, is its lower computational effort. On the other hand, its main disadvantage is that it is not possible to truly recover the unyielded zones or to detect precisely the yield surfaces. In practice, we recover $\|D(\mathbf{u})\| \sim \delta_0$ and some tolerance of the order of $\delta_0$ must be used to compute the yielded/unyielded interface, as it will be commented in the numerical tests (see section 5).

In the next section, the final model is derived as a combination of two ingredients: a dimensional analysis and a depth-integrated procedure in the framework of the multilayer (also called layer-averaged) approach.

## 3 Derivation of the model

First, we set the coordinate system. A tilted coordinate system $(x, z)$ is considered, with $x \in [x_0, x_{end}]$ and $z \in \mathbb{R}$, over a reference plane $\hat{b}(x)$ with a constant slope $\theta > 0$. We adopt the geophysical convention that considers negative slopes for $\theta > 0$ (see Figure 1), then we define $\hat{b}(x) = (x_{end} - x) \tan \theta$. In addition, $b(x)$ denotes a local bottom, measured in the normal direction to the reference plane $\hat{b}(x)$. The velocity vector is $\mathbf{u} = (u, w)$, where $u, w$ are its downslope and normal components, and $\nabla = (\partial_x, \partial_z)$ is the usual differential operator.

Then, system (3) reads

$$ \left\{ \begin{array}{l} \partial_x u + \partial_z w = 0, \\ \\ \rho(\partial_t u + u \partial_x u + w \partial_z u) + \partial_x p = \rho g \sin \theta + \partial_x \tau_{xx} + \partial_z \tau_{xz}, \\ \\ \rho(\partial_t w + u \partial_x w + w \partial_z w) + \partial_z p = -\rho g \cos \theta + \partial_x \tau_{zx} + \partial_z \tau_{zz}, \end{array} \right. \quad (5) $$

where $\tau_{xx}, \tau_{xz}, \tau_{zx}, \tau_{zz}$ denote the components of the deviatoric tensor $\tau$.

Concerning the boundary conditions, the usual kinematic and non-penetration conditions at the free surface and the bottom, respectively, are used

$$ \partial_t (b+h) + u_{|z=b+h} \partial_x (b+h) = w_{|z=b+h}, \quad \text{and} \quad u_{|z=b} \partial_x b = w_{|z=b}. $$

We also consider that there is no surface tension at the free surface

$$ \sigma \mathbf{n}^S = 0, \quad (6) $$

being $\mathbf{n}^S$ the normal vector at the free surface. At the bottom, a friction condition is assumed

$$ \sigma \mathbf{n}^b - ((\sigma \mathbf{n}^b) \cdot \mathbf{n}^b) \mathbf{n}^b = \left( \tau_b \frac{u_{|b}}{|u_{|b}|}, 0 \right)', \quad (7) $$
being $n^b$ the downward normal vector at the bottom, and $\tau_b$ (Pa) the friction stress, that is not necessarily constant. In this friction condition it appears the sign of the velocity at the bottom for the sake of generality, which is a difference with respect to other previous works (see e.g. [22, 3]), where it is always assumed a positive velocity, i.e., $sign(u) = u/|u| = 1$ constant. Let us also remark that a no-slip condition ($\mathbf{u} = 0$ at the bottom) could also be imposed in a weak sense, following the approach introduced in [28].

## 3.1 Dimensional analysis

The usual shallowness parameter $\epsilon = H/L$, the ratio of the characteristic height and length of the domain, is considered and the dimensional analysis in [30, 19] is assumed (tilde symbols (\tilde{\cdot})) denote dimensionless variables):

$$ 
\begin{aligned}
(x, z, t) &= (L\tilde{x}, H\tilde{z}, (L/U)\tilde{t}), \quad h = H\tilde{h}, \quad \rho = \rho_0\tilde{\rho}, \\
(u, w) &= (U\tilde{u}, \varepsilon U\tilde{w}), \quad p = \rho_0 g \cos \theta H\tilde{p},
\end{aligned}
$$

and

$$ D(\mathbf{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}. $$

It leads to

$$ \eta = K \left(\frac{U}{H}\right)^{n-1} \tilde{\eta} \quad \text{with} \quad \tilde{\eta} = \frac{\text{HB}}{\|\widetilde{D}(\mathbf{u})\|} + 2^n \|\widetilde{D}(\mathbf{u})\|^{n-1}, $$

where HB = $\tau_y H^n/(KU^n)$ is the Herschel-Bulkley number. We recall that $\tau_y$ and K have dimensions Pa and Pa$^n$ respectively. Note that previous hypotheses are equivalent to consider

$$ (\tau_{xx}, \tau_{xz}, \tau_{zz}) = K \left(\frac{U}{H}\right)^n (\varepsilon\tilde{\tau}_{xx}, \tilde{\tau}_{xz}, \varepsilon\tilde{\tau}_{zz}), $$

where

$$ \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}. $$

As mentioned earlier, the two regimes in the definition of the deviatoric stress tensor (2) result in a lower sheared layer and an upper non-sheared (or plug) layer in the flow. It should be mentioned that these assumptions hold in the lower sheared layer, where the strain-rate is $O(1)$. However, in transient non-uniform flows, a pseudoplug layer appears instead of a true plug layer (see [8]). Moreover, the strain-rate becomes $O(\varepsilon)$ in this pseudoplug layer. A more appropriate scaling can then be considered using the following hypothesis:

$$ \tau_{xz} = K \left(\frac{U}{H}\right)^n \varepsilon^\chi \tilde{\tau}_{xz}, \quad \text{and} \quad \tilde{\tau}_{xz} = \frac{\tilde{\eta}}{2} (\partial_{\tilde{z}} \tilde{u} + \varepsilon^2 \partial_{\tilde{x}} \tilde{w}), $$

with $\chi = 1$ in the pseudoplug layer and $\chi = 0$ otherwise.

Defining now the non-dimensional Reynolds and Froude numbers

$$ Re = \frac{\rho_0 U^{2-n} H^n}{K}, \quad Fr = \frac{U}{\sqrt{gH \cos \theta}}, $$
the non-dimensional Navier-Stokes system is written as (tildes have been dropped for simplicity)

$$
\left\{
\begin{aligned}
\partial_x u + \partial_z w &= 0, \\
\partial_t u + u \partial_x u + w \partial_z u + \frac{1}{Fr^2} \partial_x p &= \frac{1}{\varepsilon} \frac{1}{Fr^2} \tan \theta + \frac{1}{Re} \left( \varepsilon \partial_x \tau_{xx} + \frac{1}{\varepsilon^{(1-\chi)}} \partial_z \tau_{xz} \right), \\
\varepsilon^2 (\partial_t w + u \partial_x w + w \partial_z w) + \frac{1}{Fr^2} \partial_z p &= -\frac{1}{Fr^2} + \frac{1}{Re} \left( \varepsilon^{(1+\chi)} \partial_x \tau_{zx} + \varepsilon \partial_z \tau_{zz} \right).
\end{aligned}
\right.
$$

In practice, we consider two different models. The first is the leading order model up to $O(\varepsilon)$, which cannot reproduce the pseudoplug layer but a plug layer. In that case, we start from the first order system

$$
\left\{
\begin{array}{l}
\partial_x u + \partial_z w = 0, \\
\\
\partial_t u + u \partial_x u + w \partial_z u + \frac{1}{Fr^2} \partial_x p = \frac{1}{\varepsilon} \frac{1}{Fr^2} \tan \theta + \frac{1}{\varepsilon Re} \partial_z \left( \frac{\eta}{2} \partial_z u \right), \quad (8a) \\
\\
\partial_z p = -1,
\end{array}
\right.
$$

where we consider the leading order approximation $||D(\mathbf{u})|| \approx |\partial_z u|/2$ in the viscosity definition (4), yielding to

$$
\eta = \frac{\tau_y + K |\partial_z u|^n}{\sqrt{|\partial_z u|^2 / 4 + \delta_0^2}}. \qquad (8b)
$$

In order to obtain a model reproducing the pseudoplug layer, it is necessary to consider the model up to $O(\varepsilon^2)$, where normal stress contributions $(\tau_{xx}, \tau_{zz})$ are included [8, 30, 19]. It can be achieved by considering the model up to $O(\varepsilon^2)$ in both layers. However, notice that in the sheared layer there is no reason to neglect the term $\varepsilon^{(1+\chi)}\partial_x\tau_{zx}$, where $\chi=0$, which leads to a very complicated model with third order derivatives. Then, we shall consider an intermediate model including just the normal stress contributions within each layer, which satisfies a dissipative energy balance. This model can be seen as a second order perturbation of model (8) in the sheared layer, allowing us to recover the pseudoplug layer. In this case, the starting point is the first order system

$$
\left\{
\begin{array}{l}
\partial_x u + \partial_z w = 0, \\
\\
\partial_t u + u \partial_x u + w \partial_z u + \frac{1}{Fr^2} \partial_x p = \frac{1}{\varepsilon} \frac{1}{Fr^2} \tan \theta + \frac{\varepsilon}{Re} \partial_x (\eta \partial_x u) + \frac{1}{\varepsilon^{(1-\chi)} Re} \partial_z \left(\frac{\eta}{2} \partial_z u\right), \\
\\
\partial_z \left(p + \frac{\varepsilon F r^2}{Re} \eta \partial_x u\right) = -1,
\end{array}
\right.
\qquad (9a)
$$

where we have used the divergence free condition in the vertical momentum equation.
Concerning the strain-rate $||D||$, it reads, in non-dimensional form,

$$
\|D\| = \sqrt{0.5D : D} = \frac{1}{2} \sqrt{\varepsilon^{2\chi} \left( (\partial_z u)^2 + 2\varepsilon^2 \partial_x w \partial_z u + \varepsilon^4 (\partial_x w)^2 \right) + 4\varepsilon^2 (\partial_x u)^2}.
$$

Then, we consider here the enhanced first order approximation

$$
\|D_{\varepsilon}\| = \frac{1}{2} \sqrt{\varepsilon^{2\chi} (\partial_z u)^2 + 4\varepsilon^2 (\partial_x u)^2}. \qquad (9b)
$$
Note that this approximation follows the same idea of considering only those first order terms from the normal stress contributions. So, the viscosity coefficient is defined by

$$ \eta = \frac{\tau_y + 2^n K \|D_\epsilon\|^n}{\sqrt{\|D_\epsilon\|^2 + \delta^2}}. \qquad (9c) $$

In the following, the multilayer approach is applied to systems (8) and (9) to obtain the hydrostatic Multilayer Herschel-Bulkley models up to first and second order in the shallowness parameter $\epsilon$, respectively.

## 3.2 Multilayer Herschel-Bulkley models

Let us briefly recall the multilayer notations (see Figure 1). Following [29] the domain $\Omega$ is split into $N$ vertical layers $\Omega_\alpha$, for $\alpha = 1, \dots, N$. We denote by $z_{\alpha+1/2}$ the interface between layers $\Omega_\alpha$ and $\Omega_{\alpha+1}$, and then

$$ \Omega_\alpha = \{(x, z) \in [x_0, x_{end}] \times \mathbb{R}^+ / z_{\alpha-1/2} < z < z_{\alpha+1/2}\}. $$

Note that $z_{1/2}$ and $z_{N+1/2}$ are the bottom and free surface levels, respectively. We shall remark that they are virtual layers without a physical meaning, contrary to the stratified flows approach. Considering $h_\alpha = z_{\alpha+1/2} - z_{\alpha-1/2}$ the height of the layer $\Omega_\alpha$, and therefore $h = \sum_\alpha h_\alpha$ is the total height, the vertical mesh is defined by the positive coefficients $l_\alpha > 0$ satisfying

$$ h_\alpha = l_\alpha h, \quad \text{and} \quad \sum_{\alpha=1}^{N} l_\alpha = 1. $$

<figure><img src="image_2.png" /><figcaption>Figure 1: Sketch of the multilayer configuration and notation.</figcaption></figure>

In addition, given an arbitrary function $f$ we define

$$ f_{\alpha+1/2}^{-} = \lim_{\substack{z \to z_{\alpha+1/2} \\ z < z_{\alpha+1/2}}} f_{|\Omega_\alpha}, \quad f_{\alpha+1/2}^{+} = \lim_{\substack{z \to z_{\alpha+1/2} \\ z > z_{\alpha+1/2}}} f_{|\Omega_{\alpha+1}}, $$

the lower and upper limits of $f$ at the interface given by $z_{\alpha+1/2}$. In the case of a continuous function $f$, we simply denote by $f_{\alpha+1/2}$ its approximation at the interface. We also denote by $u_\alpha$ the averaged velocity in the layer $\Omega_\alpha$

$$ u_\alpha = \frac{1}{h_\alpha} \int_{z_{\alpha-1/2}}^{z_{\alpha+1/2}} u(t, x, z) dz. $$
Thus, the multilayer approach leads to a piecewise constant profile of velocity, which may be discontinuous at the interfaces $z_{\alpha+1/2}$. In order to approximate the viscosity and the components of $\tau$ at the interfaces, appropriate definitions of the deviatoric stress tensor components at the interfaces are introduced, based on the normal jump condition for the mass and momentum equations (see [27] for details).

It is also necessary to define $\mathcal{U}_{\bar{Z},\alpha+1/2}^H$, an approximation of $\partial_z u$ at the interface $z_{\alpha+1/2}$. In particular, we set

$$ \mathcal{U}_{\bar{Z},\alpha+1/2}^H = \frac{u_{\alpha+1} - u_{\alpha}}{h_{\alpha+1/2}}, \quad \text{for } \alpha = 1, \dots, N-1, $$

with $h_{\alpha+1/2} = (h_{\alpha} + h_{\alpha+1})/2$.

### 3.2.1 First order multilayer model

Following a layer-averaging procedure of system (8) as in [27] (see their Appendix A), the final multilayer model up to first order reads (in original variables)

$$ \left\{ \begin{array}{l} l_{\alpha} \left( \partial_t h + \partial_x (h u_{\alpha}) \right) = G_{\alpha+1/2} - G_{\alpha-1/2}, \\ l_{\alpha} \left( \partial_t (h u_{\alpha}) + \partial_x (h u_{\alpha}^2) + g \cos \theta h \partial_x (z_b + h) \right) = \frac{1}{\rho} (K_{\alpha-1/2} - K_{\alpha+1/2}) \quad (10) \\ + \frac{1}{2} G_{\alpha+1/2} (u_{\alpha+1} + u_{\alpha}) - \frac{1}{2} G_{\alpha-1/2} (u_{\alpha} + u_{\alpha-1}), \end{array} \right. $$

where $z_b = b + \hat{b}$ is the bottom topography. The term $G_{\alpha+1/2}$ is the mass transference between the layers $\Omega_\alpha$ and $\Omega_{\alpha+1}$, which by combining mass equations in previous system can be written as

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

Finally, the viscous term $K_{\alpha+\frac{1}{2}}$ is defined by

$$ K_{\alpha+1/2} = -\frac{1}{2} \eta_{\alpha+1/2} \mathcal{U}_{\bar{Z},\alpha+1/2}^H, \quad \alpha = 1, \dots, N-1, \quad (12) $$

where the viscosity coefficient at the interface $z_{\alpha+1/2}$ is given by

$$ \frac{1}{2} \eta_{\alpha+\frac{1}{2}} = \frac{\tau_y + K \left| \mathcal{U}_{\bar{Z},\alpha+1/2}^H \right|^n}{\sqrt{\left| \mathcal{U}_{\bar{Z},\alpha+1/2}^H \right|^2 + \delta^2}}, \quad (13) $$

with $\delta = 2\delta_0$.

The terms $K_{1/2}$ and $K_{N+1/2}$ are defined by the boundary conditions at the bottom and the free surface, respectively. Using the boundary condition (6) we set $K_{N+1/2} = 0$. In order to impose the friction condition (7), we use that $\text{sign}(u_1) = \text{sign}(\mathcal{U}_{\bar{Z},1/2}^H) = \mathcal{U}_{\bar{Z},1/2}^H / \sqrt{\mathcal{U}_{\bar{Z},1/2}^H + \delta^2}$, where the denominator in the sign function has been regularized, and the following definition is considered:

$$ K_{1/2} = -\frac{1}{2} \eta_{1/2} \mathcal{U}_{\bar{Z},1/2}^H, \quad \text{with} \quad \frac{1}{2} \eta_{1/2} = \frac{\tau_b}{\sqrt{\left| \mathcal{U}_{\bar{Z},1/2}^H \right|^2 + \delta^2}}, \quad (14) $$
where we recall that $\tau_b$ is the friction coefficient. For a general friction law, which is compatible with the rheology, we can consider $\tau_b$ defined by

$$ \tau_b = \tau_y + K_b \left| \mathcal{U}_{\mathcal{Z},1/2}^H \right|^n, \qquad (15) $$

being $K_b$ a calibrated coefficient, which may be variable. The goal of considering this general coefficient is the fact that the coefficients defining the friction law can be different from those defining the viscosity/consistency. This coefficient should be calibrated, very often by means of laboratory experiments. For instance, the Coussot's closure [22, 3] is defined by

$$ \tau_b = \tau_y \left(1 + 1.93G^{3/10}\right), \quad \text{with} \quad G = \left(\frac{K}{\tau_y}\right)^3 \left|\mathcal{U}_{\mathcal{Z},1/2}^H\right|, \qquad (16) $$

which is an experimental formula calibrated for shallow flows and for a kaolin fluid, where $n = 0.3$. This definition of $\tau_b$ coincides with (15) in the case $n = 0.3$ when setting $K_b = 1.93K (\tau_y/K)^{1/10}$.

In the numerical experiments, we will use a value $n = 0.33 \approx 1/3$. Note that, if we replace 0.3 by $n = 1/3$ in (16), we obtain $\tau_b$ given by (15) with $K_b = 1.93K$, that is,

$$ \tau_b = \tau_y + 1.93K \left| \mathcal{U}_{\mathcal{Z},1/2}^H \right|^n. \qquad (17) $$

In previous expressions for friction conditions, we define $\mathcal{U}_{\mathcal{Z},1/2}^H = u_1/h_1$. Alternatively, we can also consider a no-slip condition, where in this case

$$ \tau_b = \tau_y + K \left| \mathcal{U}_{\mathcal{Z},1/2}^H \right|^n, \quad \text{with} \quad \mathcal{U}_{\mathcal{Z},1/2}^H = \frac{2u_1}{h_1}. \qquad (18) $$

The factor 2 in previous equation is a consequence of the strategy followed in [28] to weakly impose friction or no-slip boundary conditions. This strategy is based on a ghost cell technique (see its subsection 2.1.2 and its Figure 7). Concretely, it consists in approximating $\mathcal{U}_{\mathcal{Z},1/2}^H = (u_1 - u_0)/h_1$ where $u_0$ is a virtual velocity that takes the value $u_0 = -u_1$ (no slip) or $u_0 = 0$ (friction). Then, the following definition is considered

$$ \mathcal{U}_{\mathcal{Z},1/2}^H = \frac{\lambda u_1}{h_1} \quad \text{where} \quad \begin{cases} \lambda = 2 & \text{if no-slip,} \\ \lambda = 1 & \text{if friction.} \end{cases} \qquad (19) $$

These friction or no-slip conditions will be used in the numerical tests of subsection 5.3.

Concerning the energy balance of the resulting model, it is easy to see that the following result holds:

**Theorem 1.** Defining the energy of the layer $\Omega_\alpha$ as

$$ E_\alpha = h_\alpha \left( \frac{|u_\alpha|^2}{2} + g \cos \theta \left( z_b + \frac{h}{2} \right) \right), \qquad (20a) $$

for $\alpha = 1, \dots, N$, system (10) satisfies the dissipative energy balance

$$ \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} \right) \right] \\ \leq -\frac{\eta_{1/2}}{2} \frac{\lambda u_1^2}{h_1} - \sum_{\alpha=1}^{N-1} \frac{\eta_{\alpha+\frac{1}{2}} (u_{\alpha+1} - u_\alpha)^2}{2 h_{\alpha+\frac{1}{2}}}, \qquad (20b) $$

where $\lambda = \{1, 2\}$ depending on the friction condition (19).
### 3.2.2 Second order multilayer model with pseudoplug

We remark that considering all the $O(\varepsilon)$ terms over the whole domain leads to a very complex model, involving third order spatial derivatives (see [16]). The treatment of such a resulting system, as well as the inclusion of non-hydrostatic effects, would be interesting, although it is beyond the scope of this paper.

In order to show here the ability of the multilayer approach to recover this pseudoplug layer in non-uniform flows, we discretize system (9), which includes the normal stress contributions, in the multilayer framework. Thus, the following system is obtained:

$$ egin{cases} l_{\alpha} \left( \partial_t h + \partial_x (hu_{\alpha}) \right) = G_{\alpha+1/2} - G_{\alpha-1/2}, \\ l_{\alpha} \left( \partial_t (hu_{\alpha}) + \partial_x (hu_{\alpha}^2) + g \cos \theta h \partial_x (z_b + h) \right) = \frac{1}{\rho} l_{\alpha} \partial_x (2h\eta_{\alpha} \partial_x u_{\alpha}) \\ + \frac{1}{\rho} (K_{\alpha-1/2} - K_{\alpha+1/2}) + \frac{1}{2} G_{\alpha+1/2} (u_{\alpha+1} + u_{\alpha}) - \frac{1}{2} G_{\alpha-1/2} (u_{\alpha} + u_{\alpha-1}), \end{cases} \quad (21) $$

where $G_{\alpha+1/2}$ and $K_{\alpha+1/2}$ are given by (11) and (12), respectively. Note that in this case definitions of the viscosity at the midpoint of each layer and at the interfaces are needed, which basically depend on the definition of $||D_{\varepsilon}||$ (9b). So, the viscosity coefficients are defined by

$$ \eta_{\alpha} = \frac{\tau_y + 2^n K ||D_{\varepsilon,\alpha}||^n}{\sqrt{||D_{\varepsilon,\alpha}||^2 + \delta^2}}, \quad \text{and} \quad \eta_{\alpha+\frac{1}{2}} = \frac{\tau_y + 2^n K ||D_{\varepsilon,\alpha+1/2}||^n}{\sqrt{||D_{\varepsilon,\alpha+1/2}||^2 + \delta^2}}. $$

Finally, the boundary condition $K_{1/2}$ is defined as in (14) where

$$ \eta_{1/2} = \frac{\tau_b}{\sqrt{||D_{\varepsilon,1/2}||^2 + \delta^2}}, $$

and $\tau_b$ defined as in previous subsection, where $||D_{\varepsilon,1/2}||$ instead of $|u_{\alpha}^H|_{1/2}$ must be considered.

Concerning the energy balance satisfied by this model, it reads

$$ \rho \partial_t \left( \sum_{\alpha=1}^{N} E_{\alpha} \right) + \partial_x \left[ \sum_{\alpha=1}^{N} u_{\alpha} \left( \rho E_{\alpha} + \rho g \cos \theta h_{\alpha} \frac{h}{2} - 2h_{\alpha} \eta_{\alpha} \partial_x u_{\alpha} \right) \right] \\ \leq -\frac{\eta_{1/2}}{2} \frac{\lambda u_1^2}{h_1} - \sum_{\alpha=1}^{N} 2h_{\alpha} \eta_{\alpha} (\partial_x u_{\alpha})^2 - \sum_{\alpha=1}^{N-1} \frac{\eta_{\alpha+\frac{1}{2}} (u_{\alpha+1} - u_{\alpha})^2}{h_{\alpha+\frac{1}{2}}}, $$

where $E_{\alpha}$ is given by (20a). In order to obtain such an energy balance for this model, it is essential to consider $K_{\alpha+1/2}$ defined by (12), although other terms of order $O(\varepsilon)$ are neglected here.

Model (21) will be considered in the numerical tests' section to illustrate the ability of the multilayer approach to reproduce the pseudoplug layer in non-uniform flows.

## 3.3 Steady states at rest

Let us analyze now the steady solutions of system (10). Concretely, we focus here on those steady states at rest, that is, with $u_{\alpha} = 0$ for $\alpha = 1, \dots, N$. Actually, they are also steady-at-rest states for system (21). In this case, we obtain

$$ l_{\alpha} \rho g h \cos \theta |\partial_x (z_b + h)| = K_{\alpha-1/2} - K_{\alpha+1/2}. $$
Summing up previous expression for $\alpha = 1, \dots, N$, and using boundary conditions ($K_{N+1/2} = 0$) we have

$$ \rho g h \cos \theta |\partial_x (z_b + h)| = K_{1/2} \le \tau_y, $$

where last inequality comes from the definition of the stress tensor (2) and the fact that $K_{1/2}$ approximates $\tau_{xz|_b}$ in the multilayer framework. Therefore, the steady solutions we are interested in, are those verifying

$$ u_{\alpha} = 0 \text{ for } \alpha = 1, \dots, N, \quad \text{and} \quad \rho g h \cos \theta |\partial_x (z_b + h)| \le \tau_y. \quad (22) $$

Note that they are also a family of steady states for the one-layer shallow model and not only in the multilayer case.

Let us consider now some particular cases that we will deal with in the numerical tests. Let us simplify by assuming a flat local bottom $\partial_x b = 0$.

First, we consider a uniform flow $\partial_x (z_b + h) = -\tan \theta$ (and therefore $\partial_x h = 0$), then condition (22) reads

$$ h \le \frac{\tau_y}{\rho g \sin \theta}. $$

It is essential to notice the dependence on $h$ in previous inequality, which is a difference from steady states for other materials, such as granular flows (see e.g. [35]). It induces that for fixed slope and rheological properties, the material flows or not depending on the thickness of the uniform layer.

Second, in the case where $\partial_x h \ne 0$ but the local bottom $b$ remains constant ($\partial_x b = 0$), we obtain the family of steady states at rest given by

$$ h \left| -\tan \theta + \partial_x h \right| \le \frac{\tau_y}{\rho g \cos \theta}. \quad (23) $$

Actually, in the limit case, the solution of the initial value problem

$$ \left\{ \begin{array}{l} h \left| -\tan \theta + \partial_x h \right| = \frac{\tau_y}{\rho g \cos \theta}, \\ h(x_0) = h_0, \end{array} \right. $$

where we recall that $x_0$ is the left boundary of the domain, is a steady state. Assuming $(-\tan \theta + \partial_x h) < 0$, this solution takes the form

$$ h(x) = \begin{cases} \frac{\gamma}{\kappa} \left(1 + \frac{\gamma}{\kappa} W_0(z)\right) & \text{if } -e^{-1} < z < 0, \\ 0 & \text{otherwise,} \end{cases} \quad \text{where} \quad z = \frac{-1}{\gamma} e^{\frac{\kappa^2(\zeta-x)}{\gamma}} - 1, \quad (24a) $$

with $W_0(z)$ the main branch of the Lambert's W-function and

$$ \gamma = \frac{-\tau_y}{\rho g \cos \theta}, \quad \kappa = -\tan \theta, \quad \text{and} \quad \zeta = \frac{\gamma}{\kappa^2} \left( 1 + \log \left( (\gamma - \kappa h_0) e^{\kappa h_0 / \gamma} - 1 \right) \right). \quad (24b) $$

# 4 A well-balanced numerical discretization

We focus here on the discretization of system (10). If the system (21) is considered, only the last step of the scheme should be modified accordingly (see Remark 1).

First, note that system (10) cannot be written as a hyperbolic system of conservation laws, in particular because of the viscous terms $K_{\alpha+1/2}$ in the momentum equations. Moreover, the mass transference terms lead to non-conservative products. Therefore, multilayer
models are usually discretized by a splitting approach where a finite volume method is considered in a first step, whereas viscous contributions are discretized semi-implicitly in a second stage for stability reasons (see [6, 29] among others).

The numerical scheme employed here is very similar to the one in [28] for a multilayer system for dry granular flows. Before giving details, let us summarize the key-points to obtain a well-balanced discretization, which is the main goal. Actually, the scheme is designed to preserve the steady solutions (22) in which we are interested in. It is achieved as a result of the combination of both steps in the splitting procedure. First, a hydrostatic reconstruction based on the friction term is considered to make the numerical diffusion vanishing when the steady state condition is satisfied, making possible to obtain $\partial_t h = 0$ in such cases. The second step is linked to the fact that one wants the scheme to produce a vanishing velocity when converging to the steady states at rest. However, by the structure of the problem, the linear system associated to the velocities (see (28) below) at steady states is such that it naturally leads to obtain very small values of the velocities. Indeed, remark that the coefficients of the system's matrix are of order $\delta^{-1}$ at cells where the strain rate is of order $\delta$. We remark that we cannot obtain zero velocity because of the regularization technique, but a residual velocity that depends on the regularization parameter $\delta_0$, as will be analyzed in subsection 5.2.

Let us define $\mathbf{w} = (h, q_1, q_2, ..., q_N)' \in \mathbb{R}^{N+1}$ the vector of the conservative variables, where $q_\alpha = hu_\alpha$ is the discharge of the layer $\Omega_\alpha$. Using this vector, system (10) can be rewritten in a compact form as

$$ \partial_t \mathbf{w} + \partial_x \mathbf{F}_c(\mathbf{w}) + \mathbf{S}(\mathbf{w})\partial_x (z_b + h) + \mathbf{B}(\mathbf{w})\partial_x \mathbf{w} = \mathbf{E}(\mathbf{w}), \quad (25) $$

where $\mathbf{F}_c, \mathbf{S} \in \mathbb{R}^{N+1}$ represent convective and pressure terms, respectively, and whose components, for $j = 0, 1, ..., N$, are

$$ F_{c,j} = \begin{cases} \sum_{\beta=1}^{N} l_{\beta} q_{\beta}, & \text{if } j = 0, \\ \frac{q_j^2}{h}, & \text{otherwise,} \end{cases} \qquad S_j = \begin{cases} 0, & \text{if } j = 0, \\ gh \cos \theta, & \text{otherwise.} \end{cases} $$

The matrix $\mathbf{B} = (b_{jk}) \in \mathcal{M}_{N+1}(\mathbb{R})$ contains all the non-conservative contributions, i.e., the terms coming from the mass transference terms. Concretely, it is defined by

$$ b_{jk} = \begin{cases} 0, & \text{if } jk = 0, \\ \frac{1}{2hl_j} (q_j + q_{j-1}) \xi_{j-1,k} - \frac{1}{2hl_j} (q_{j+1} + q_j) \xi_{j,k}, & \text{otherwise,} \end{cases} $$

with

$$ \xi_{j,k} = \begin{cases} (1 - (l_1 + \cdots + l_j))l_k, & \text{if } k \le j, \\ -(l_1 + \cdots + l_j)l_k, & \text{otherwise,} \end{cases} $$

for $j, k \in \{1, ..., N\}$. Finally, $\mathbf{E} \in \mathbb{R}^{N+1}$ collects the viscous terms, with

$$ E_j = \begin{cases} 0, & \text{if } j = 0, \\ \frac{1}{\rho l_j} (K_{j-1/2} - K_{j+1/2}), & \text{otherwise.} \end{cases} $$

Once we have the model as a system of conservation laws with non-conservative products and source terms, we move to the finite volume framework. We consider the usual
uniform subdivision of the domain in control volumes $V_i = [x_{i-1/2}, x_{i+1/2}]$, for $i \in \mathcal{I}$ (the set of ad hoc indices), where the constant mesh step is $\Delta x$. The center of each volume is denoted by $x_i = (x_{i-1/2} + x_{i+1/2})/2$. Then, for a fixed time $t^m$, we denote the averaged conservative variable vector

$$ \boldsymbol{w}_i^m = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} \boldsymbol{w}(t^m, x)dx. $$

In the first step, the viscous contributions are removed, i.e. $\boldsymbol{E} = \mathbf{0}$, and the discrete form of (25) in this case reads

$$ \boldsymbol{w}_i^{m+1/2} = \boldsymbol{w}_i^m - \frac{\Delta t}{\Delta x} \left( \mathcal{F}_{c,i+1/2}^m - \mathcal{F}_{c,i-1/2}^m + \frac{1}{2} \left( \mathcal{B}_{i+1/2}^m + \mathcal{B}_{i-1/2}^m + \mathcal{S}_{i+1/2}^m + \mathcal{S}_{i-1/2}^m \right) \right), $$

being

$$ \mathcal{B}_{i+1/2} = \frac{1}{2}(\mathbf{B}(\boldsymbol{w}_{i+1}) + \mathbf{B}(\boldsymbol{w}_i))(\boldsymbol{w}_{i+1} - \boldsymbol{w}_i), \quad \text{and} \quad \mathcal{S}_{i+1/2} = \frac{1}{2}(\mathbf{S}(\boldsymbol{w}_{i+1}) + \mathbf{S}(\boldsymbol{w}_i))\Delta\varphi_{i+1/2}, $$

where $\Delta\varphi_{i+1/2}$ is the free surface variation

$$ \Delta\varphi_{i+1/2} = (h_{i+1} + z_{b,i+1} - Z_M)_+ - (h_i + z_{b,i} - Z_M)_+, $$

with $Z_M = \max(z_{b,i}, z_{b,i+1})$ and $(\cdot)_+$ the positive part. Finally, the numerical flux $\mathcal{F}_{c,i+1/2}$ used for the convective terms, in general, is written as

$$ \mathcal{F}_{c,i+1/2} = \frac{1}{2} (\mathbf{F}_c(\boldsymbol{w}_i) + \mathbf{F}_c(\boldsymbol{w}_{i+1})) - \frac{1}{2} \mathcal{D}_{i+1/2}, $$

being $\mathcal{D}_{i+1/2}$ the numerical diffusion, which actually determines the scheme. We used here a HLL-type approximated Riemann solver in the framework of the PVM schemes introduced in [18], where the numerical diffusion $\mathcal{D}_{i+1/2}$ is a polynomial evaluation of the Roe matrix.

A key point of the scheme is the definition of the numerical diffusion to preserve the steady water-at-rest solutions described by (22), in particular those with $\partial_x (z_b + h) \neq 0$. The result we look for is to cancel the numerical diffusion when we predict a state-at-rest at $x_{i+1/2}$ at the next time step. To do so, we evaluate if the sum of the pressure and convective terms is lower than the friction force at the current time step. To this aim, a key point is the definition of $f_{i+1/2}$ (27). It leads to ensure $\partial_t h = 0$ and therefore the well-balanced property of the scheme. This will be shown numerically in subsection 5.2.

To this end, we consider the hydrostatic reconstruction with source term (see [15]) taking into account the friction term. So, denoting by $S_L$ and $S_R$ the approximations of the minimum and maximum wave speeds with a baroclinic approximation (see [28]), we define

$$ \mathcal{D}_{i+1/2} = a_0 (\hat{\boldsymbol{w}}_{i+1} - \hat{\boldsymbol{w}}_i) + a_1 (\mathbf{F}_c(\boldsymbol{w}_{i+1}) - \mathbf{F}_c(\boldsymbol{w}_i) + \mathcal{B}_{i+1/2}), \quad (26) $$

with

$$ a_0 = \frac{S_R |S_L| - S_L |S_R|}{S_R - S_L}, \quad a_1 = \frac{|S_R| - |S_L|}{S_R - S_L}, $$

and $\hat{\boldsymbol{w}}_i, \hat{\boldsymbol{w}}_{i+1}$ the reconstructed states. Notice that we have removed the pressure term in the first order contribution of $\mathcal{D}_{i+1/2}$ and therefore it slightly differs from a usual path-conservative scheme. The reason is that it is a more stable choice.

The reconstructed states $\hat{\boldsymbol{w}}_i, \hat{\boldsymbol{w}}_{i+1}$ are defined as $\boldsymbol{w}_i, \boldsymbol{w}_{i+1}$ except for the first components $h_i, h_{i+1}$, which are replaced by

$$ h_{i+1/2-} = (h_i - (\Delta Z_{i+1/2})_+)_+, \quad h_{i+1/2+} = (h_{i+1} - (-\Delta Z_{i+1/2})_+)_+, $$
where

$$ \Delta Z_{i+1/2} = z_{b,i+1} - z_{b,i} + \Delta C_{i+1/2}, $$

with $\Delta C_{i+1/2} = f_{i+1/2}\Delta x_{i+1/2}/(g \cos \theta)$ accounting for the well-balanced correction associated to the friction term (15) or (16). Concretely, we set

$$ f_{i+1/2} = \begin{cases} U^* & \text{if } h_{i+1/2} |U^*| \le \tau_b/\rho, \\ 0 & \text{otherwise,} \end{cases} \qquad (27) $$

with

$$ U^* = \frac{u_{1,i+1/2}}{\Delta t} - \frac{g \cos \theta \Delta \varphi_{i+1/2}}{\Delta x}, $$

and $u_{1,i+1/2} = (u_i + u_{i+1})/2$, $h_{i+1/2} = (h_i + h_{i+1})/2$. Note that $U^*$ is an estimation of the velocity after the first step of the scheme $u_{i+1/2}^{m+1/2}$. Then, in practice (27) makes the balance between convective/pressure and friction terms and defines $h_{i+1/2\pm}$ to have $h_{i+1/2+} - h_{i+1/2-} = 0$ in (26) if the steady state condition (22) holds at the discrete level at the interface $x_{i+1/2}$.

Once we have solved the hyperbolic system, obtaining the intermediate state $w_i^{m+1/2}$, $i \in \mathcal{I}$, the viscous terms are discretized semi-implicitly. From the definition of $\mathcal{E}$, we trivially have $h_i^{m+1} = h_i^{m+1/2}$, and for the discharge we solve at each volume the $N \times N$ tridiagonal linear system

$$ q_{\alpha,i}^{m+1} = q_{\alpha,i}^{m+1/2} + \frac{\Delta t}{\rho l_{\alpha} h_{i}^{m+1}} \left( \frac{\eta_{\alpha+\frac{1}{2}}^{m} q_{\alpha+1,i}^{m+1} - q_{\alpha,i}^{m+1}}{2 l_{\alpha+\frac{1}{2}} h_{i}^{m+1}} - \frac{\eta_{\alpha-\frac{1}{2}}^{m} q_{\alpha,i}^{m+1} - q_{\alpha-1,i}^{m+1}}{2 l_{\alpha-\frac{1}{2}} h_{i}^{m+1}} \right), \quad (28) $$

where $\eta_{\alpha+1/2}$ is given by (13). Notice that for the particular cases $\alpha = 1$ and $\alpha = N$, previous equation should be modified according to boundary conditions $K_{N+1/2} = 0$ and $K_{1/2}$ defined by (14).

Finally, let us remark that dealing with dry areas is always a hard issue from the numerical point of view. Thus, the second step of the scheme is modified in case of wet/dry fronts. In practice, if the height $h_i^{m+1}$ is lower than a tolerance, which we set as $10^{-6}$, then we discretize the viscous terms there as explained but considering a one-layer model. The motivation is that a multilayer description makes no sense in case of too small material thicknesses, then we simply treat it as the one-layer case.

**Remark 1.** In the case of considering the multilayer model (21) with the normal stress contributions, the second step must be accordingly modified in order to discretize also the term $\partial_x (2h_\alpha \eta_\alpha \partial_x u_\alpha)$ in a semi-implicit way. In that case, the $N \times N$ tridiagonal linear system (28), which is solved in each control volume $V_i$ for $i = 1, \dots, N_x$, is replaced by a pentadiagonal linear system with $N \times N_x$ equations and unknowns that is solved for the whole domain. It leads to a huge increase of the computational cost compared to the algorithm for model (10). The linear system is solved with a biconjugate gradient stabilized method, although other linear solvers can be used.

# 5 Numerical tests

In this section we present some numerical tests to evaluate the proposed multilayer model, comparing it with analytical solutions and experimental data. Unless otherwise stated, model (10) is considered. In subsection 5.1 a test with the analytical solution of a uniform flow is presented. The ability of multilayer approach to reproduce the sheared/pseudoplug
layers is shown. In subsection 5.2 the influence of the well-balanced treatment described in previous section is analyzed. We study the accuracy of the numerical scheme preserving the steady states at rest with non-planar free surface. Later, in subsection 5.3, we compare the multilayer results with the experimental data of [4, 3] viscoplastic dam breaks and with the results of lubrication models.

Unless otherwise specified, the material properties are

$$ \rho = 1000 \text{ Kg/m}^3, \quad \tau_y = 33 \text{ Pa}, \quad K = 26 \text{ Pa s}^n, \quad n = 0.33, \qquad (29) $$

corresponding to Carbopol ultrez 10 at a mass concentration of 0.15% (see [3, 4] for details of the whole experimental setup). We fix CFL = 0.5 and $\delta = 10^{-5}$.

Concerning our choice of this regularization parameter, it agrees with other works for viscoplastic fluids, as [34, 20, 42]. In these works, authors analyzed deeply the influence of $\delta$ on the convergence, and also its connection with the spatial discretization. In subsection 5.1 we use the first order HLL scheme described in section 4 to deal with a uniform flow. While in subsections 5.2 and 5.3, its second order version is considered in order to reduce the numerical diffusion in more general flow configurations. The second order discretizations are based on a MUSCL reconstruction for space (see e.g. [17]) and the Heun method for the time.

When showing errors, they are computed as follows. For any generic variable $w$, being $w_{\alpha,i}$ the approximated value in the volume $I_i$ and layer $\Omega_\alpha$ (or $w_{\alpha,i} \equiv w_i$ if the variable does not depend on $z$), we consider the $L_1$, $L_2$, and $L_\infty$ norms

$$ \|w\|_1 = \Delta x \sum_{i=1}^{N_x} \frac{h_i}{N} \sum_{\alpha=1}^{N} |w_{\alpha,i}|, \quad \|w\|_2 = \sqrt{\Delta x \sum_{i=1}^{N_x} \frac{h_i}{N} \sum_{\alpha=1}^{N} w_{\alpha,i}^2}, \quad \text{and} \quad \|w\|_\infty = \max_{\substack{\alpha=1,...,N \\ i=1,...,N_x}} |w_{\alpha,i}|, $$

being $N_x$ the number of cells in the horizontal discretization and $N$ the number of vertical layers. If we measure the error in the vertical direction for a fixed control volume $I_i$, we denote

$$ \|w\|_{I_i,1} = \frac{h_i}{N} \sum_{\alpha=1}^{N} |w_{\alpha,i}|, \quad \|w\|_{I_i,2} = \sqrt{\frac{h_i}{N} \sum_{\alpha=1}^{N} w_{\alpha,i}^2}, \quad \text{and} \quad \|w\|_{I_i,\infty} = \max_{\alpha=1,...,N} |w_{\alpha,i}|. $$

Then, absolute errors are computed as

$$ \text{Err}_p[w] = \|w - \hat{w}^{ref}\|_p, \quad \text{Err}_p^z[w] = \|w - \hat{w}^{ref}\|_{I_{i,p}}, \quad \text{for } p = 1, 2, \infty, $$

where $\hat{w}_{ref}$ is the projection of the reference solution into the discrete solution space. Unless otherwise specified, absolute errors are considered.

## 5.1 Analytical solution for a uniform flow

Let us first compare with a uniform flow (no variations in the $x$-direction) over a constant slope, starting from the rest. In this configuration, the analytical solution is known (see [22, 39, 4, 32, 19] among others). Let us briefly summarize how this solution is obtained.

We consider a uniform flow with height $h$ over an inclined plane with constant slope $\theta > 0$. Under the uniform assumption and assuming $b = 0$, system (5) yields

$$ p(z) = \rho g \cos \theta (h - z), \quad \tau_{xz}(z) = \rho g \sin \theta (h - z). $$

Taking into account the definition of $\tau$ and assuming $\text{sign}(\partial_z u) > 0$, we get

$$ \tau_{xz}(z) = \tau_y + K (\partial_z u)^n, \quad \text{if } |\partial_z u| > 0, \qquad (30a) $$
and in case $|\partial_z u| = 0$, we only know that $|\tau_{xz}| < \tau_y$ holds, or equivalently

$$z > h_c, \quad \text{with} \quad h_c = h - h_p,$$

where $h_p$ is the thickness of the unsheared (or pseudoplug) layer, also called the pseudo-yield surface [8]. Then, $h_p$ is given by the yield condition $\tau_{xz} = \tau_y$, getting therefore

$$h_p = \frac{\tau_y}{\lambda K}, \quad \text{with} \quad \lambda = \frac{\rho g \sin \theta}{K}. \qquad (30b)$$

Note that the $h_c$ level defines the interface between the sheared and the pseudoplug layers. Thus, from (30a) we have

$$\partial_z u = \lambda^{1/n} (h_c - z)^{1/n}, \quad \text{if } z < h_c, \qquad (30c)$$

which is integrated to obtain the analytical profile of the velocity

$$u(z) = \begin{cases} u_{\text{plug}} \left( 1 - \left( 1 - \frac{z}{h_c} \right)^{(n+1)/n} \right), & \text{if } z < h_c, \\ u_{\text{plug}}, & \text{otherwise,} \end{cases} \qquad (30d)$$

where

$$u_{\text{plug}} = u(h_c) = \frac{n}{(n+1)} \lambda^{1/n} h_c^{(n+1)/n} \qquad (30e)$$

is the velocity in the pseudoplug layer. Here, a no-slip condition at the bottom has been used, although a basal velocity different from zero could be considered (see [38]).

From a dimensional analysis, an analytical profile that accounts for the term $\partial_x h$ can be also deduced for a non-uniform flow. This is actually the approach of the lubrication theory, where convective and time derivatives are neglected in the momentum conservation equation, whereas pressure and viscous terms are kept (see the references mentioned above).

It takes the form (30) by replacing $\lambda$ in (30b) by

$$\lambda = \frac{\rho g \cos \theta}{K} (\tan \theta - \partial_x h). \qquad (31)$$

Notice that it has no influence in this test, since $\partial_x h = 0$ holds in a uniform flow. Nevertheless, this analytical profile will be considered in subsection 5.3 for a general flow.

<table><caption>Table 1: Numerical convergence results for the vertical profile of velocity in $L_1$, $L_2$, and $L_\infty$ norms for a uniform flow, and computation time ($t_{comp}$) needed to simulate up to $t_f = 25$ s.</caption><thead><tr><th>N. layers</th><th>t<sub>comp</sub> (s)</th><th>Err<sup>z</sup><sub>1</sub>[u]</th><th>Order<sub>1</sub>[u]</th><th>Err<sup>z</sup><sub>2</sub>[u]</th><th>Order<sub>2</sub>[u]</th><th>Err<sup>z</sup><sub>∞</sub>[u]</th><th>Order<sub>∞</sub>[u]</th></tr></thead><tbody><tr><td>4</td><td>1.8</td><td>9.02×10<sup>-3</sup></td><td>-</td><td>4.13×10<sup>-2</sup></td><td>-</td><td>2.45×10<sup>-1</sup></td><td>-</td></tr><tr><td>8</td><td>2.9</td><td>2.27×10<sup>-3</sup></td><td>1.99</td><td>1.04×10<sup>-2</sup></td><td>1.99</td><td>6.64×10<sup>-2</sup></td><td>1.88</td></tr><tr><td>16</td><td>5.4</td><td>5.68×10<sup>-4</sup></td><td>2.00</td><td>2.61×10<sup>-3</sup></td><td>2.00</td><td>1.73×10<sup>-2</sup></td><td>1.94</td></tr><tr><td>32</td><td>10.9</td><td>1.42×10<sup>-4</sup></td><td>2.00</td><td>6.53×10<sup>-4</sup></td><td>2.00</td><td>4.41×10<sup>-3</sup></td><td>1.97</td></tr><tr><td>64</td><td>23.9</td><td>3.55×10<sup>-5</sup></td><td>2.00</td><td>1.63×10<sup>-4</sup></td><td>2.00</td><td>1.11×10<sup>-3</sup></td><td>1.99</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Figure 2: Comparison of the computed $u$, $|\partial_z u|$ and $\tau_{xz}$ profiles using 64 layers with analytical solution (30).</figcaption></figure>

We consider here a slope $\theta = 20^\circ$ and a flow initially at rest, whose height is $h_0(x) = 0.05 \, m$, and open boundary conditions, for a domain $x \in [0, 1]$, with a no-slip condition (18) at the bottom. We take $100(= N_x)$ nodes for the horizontal discretization and $CFL = 0.8$.

Figure 2 shows a good agreement between the computed approximation ($N = 64$ layers) and the analytical solution (30) for the velocity, strain rate and deviatoric stress tensor. Note that for $z > h_c$ we obtain that $\partial_z u \approx 0$ (its order of magnitude is around $10^{-5}$, $10^{-6}$ due to the value of $\delta = 10^{-5}$ in the regularized viscosity coefficient), where the velocity profile is constant and the deviatoric tensor is undefined. In Table 1 we show the $L_1$, $L_2$, and $L_\infty$ errors with respect to the analytical solution and the numerical convergence rates. In this test, we see that the method is second order accurate for the velocity. Let us mention that similar results are obtained in Table 1 and Figure 2 when using a smaller regularization parameter $\delta = 10^{-8}$, where $|\partial_z u| \sim 10^{-9}$ in the plug zone. Table 1 also shows the computation time required to reach $t_f = 25$ s, where we see, approximately, a linear increase of the computational effort with the number of vertical layers, similarly to what was observed in [27].

## 5.2 Steady-at-rest solutions and well-balanced discretization

In this subsection, we check the ability of the proposed scheme to preserve steady solutions given by (22). Let us first remark that this inequality provides not only a steady state but a family of such steady states. Here we deal with the limit case given by the equality. Note that if such state is preserved then any of the states in the family (22) is also preserved.

We set the domain $[0, x_R]$ with wall boundary conditions. At the bottom, a friction condition (15) with $K_b = K$ is considered. The flow is at rest at the initial time. In this case, we take as slope $\theta = 10^\circ$ and the same material properties (29), except for the yield stress taken as $\tau_y = 66 \, Pa$.

In order to perform a concrete numerical test, we consider here the case $b = 0$ and $(-\tan \theta + \partial_x h) < 0$, so the initial condition is given by the steady state (24). In practice,
we need to set the value of $h(0) = h_0$ to find $h(x)$. We set

$$h_0 = \frac{\tau_y}{\rho g \cos \theta (\tan \theta + \epsilon)},$$

with $\epsilon = 10^{-3}$. We consider two tests: Test 1 if $x_R = 0.8$ and Test 2 if $x_R = 1$. The obtained analytical free surfaces are shown in Figure 3, where we see that it exhibits a regular shape for Test 1, whereas we have a dry front in Test 2 with an infinite slope at the front position $x_f$. This point can be exactly computed from (24), yielding to

$$x_f = \zeta - \frac{\gamma}{\kappa^2} \log \gamma \approx 0.919.$$

Let us remark that $h(x)$ is the same for Test 1 and Test 2, and the only difference is the fact that, in Test 2, we see the front position $x_f$ and dry areas since $x_R > x_f$.

<figure><img src="image_4.png" /><figcaption>Figure 3: Free surface solutions for (a) Test 1 and (b) Test 2 at the initial time (black dashed line), and at time $t = 10$ s with the proposed scheme (green dotted lines). Gray dashed lines in (b) are the solutions without the well-balanced treatment in (26). Remind that $\delta = 10^{-5}$, here.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Maximum velocity ($\log_{10}(u_N)$) for the considered steady state tests for different values of the regularization parameter ($\delta$) at $t = 10$ s, and also the well-balanced discretization in (26). Inset figure in (b) Test 2 is a zoom of the front position.</figcaption></figure>

For both tests we consider 10 vertical layers. We clarify that there is no difference between the one-layer model and the multilayer model when looking at the steady-at-rest solutions, since these steady states are the same for both models as stated in subsection
<table><caption>Test 1: No dry front - WB</caption><thead><tr><th></th><th>δ = 10<sup>-5</sup></th><th>δ = 10<sup>-7</sup></th><th>δ = 10<sup>-10</sup></th><th>δ = 10<sup>-12</sup></th><th>δ = 10<sup>-15</sup></th></tr></thead><tbody><tr><td>Err<sub>1</sub>[h]</td><td>1.07×10<sup>-7</sup></td><td>1.17×10<sup>-9</sup></td><td>1.56×10<sup>-12</sup></td><td>1.98×10<sup>-14</sup></td><td>2.40×10<sup>-18</sup></td></tr><tr><td>Err<sub>2</sub>[h]</td><td>1.21×10<sup>-6</sup></td><td>1.30×10<sup>-8</sup></td><td>1.35×10<sup>-11</sup></td><td>1.44×10<sup>-13</sup></td><td>6.97×10<sup>-18</sup></td></tr><tr><td>Err<sub>∞</sub>[h]</td><td>2.63×10<sup>-5</sup></td><td>2.85×10<sup>-7</sup></td><td>2.86×10<sup>-10</sup></td><td>2.85×10<sup>-12</sup></td><td>4.86×10<sup>-17</sup></td></tr><tr><td>Err<sub>1</sub>[u]</td><td>1.07×10<sup>-8</sup></td><td>1.63×10<sup>-10</sup></td><td>3.52×10<sup>-13</sup></td><td>5.98×10<sup>-15</sup></td><td>8.28×10<sup>-18</sup></td></tr><tr><td>Err<sub>2</sub>[u]</td><td>6.42×10<sup>-8</sup></td><td>9.69×10<sup>-10</sup></td><td>2.08×10<sup>-12</sup></td><td>3.54×10<sup>-14</sup></td><td>4.89×10<sup>-17</sup></td></tr><tr><td>Err<sub>∞</sub>[u]</td><td>4.71×10<sup>-7</sup></td><td>6.78×10<sup>-9</sup></td><td>1.39×10<sup>-11</sup></td><td>2.39×10<sup>-13</sup></td><td>3.38×10<sup>-16</sup></td></tr></tbody></table>

<table><caption>Test 2: dry front - WB</caption><thead><tr><th></th><th>δ = 10<sup>-5</sup></th><th>δ = 10<sup>-10</sup></th><th>δ = 10<sup>-15</sup></th><th>NoWB<br>δ = 10<sup>-10</sup></th></tr></thead><tbody><tr><td>Err<sub>1</sub>[h]</td><td>1.37×10<sup>-7</sup></td><td>4.90×10<sup>-11</sup></td><td>3.70×10<sup>-11</sup></td><td>6.63×10<sup>-4</sup></td></tr><tr><td>Err<sub>2</sub>[h]</td><td>1.68×10<sup>-6</sup></td><td>1.13×10<sup>-9</sup></td><td>9.67×10<sup>-10</sup></td><td>2.20×10<sup>-3</sup></td></tr><tr><td>Err<sub>∞</sub>[h]</td><td>5.87×10<sup>-5</sup></td><td>3.45×10<sup>-8</sup></td><td>2.96×10<sup>-8</sup></td><td>1.79×10<sup>-2</sup></td></tr><tr><td>Err<sub>1</sub>[u]</td><td>1.12×10<sup>-8</sup></td><td>3.74×10<sup>-13</sup></td><td>1.16×10<sup>-15</sup></td><td>4.45×10<sup>-7</sup></td></tr><tr><td>Err<sub>2</sub>[u]</td><td>6.53×10<sup>-8</sup></td><td>2.27×10<sup>-12</sup></td><td>8.03×10<sup>-13</sup></td><td>1.69×10<sup>-5</sup></td></tr><tr><td>Err<sub>∞</sub>[u]</td><td>4.71×10<sup>-7</sup></td><td>5.65×10<sup>-10</sup></td><td>5.63×10<sup>-10</sup></td><td>1.49×10<sup>-3</sup></td></tr></tbody></table>

Table 2: Errors with respect to the initial state for the Test 1 and Test 2 for different values of the regularization parameter $\delta$. WB and NoWB denote a well-balanced or not well-balanced discretization.

3.3. For Test 1, we consider $N_x = 400$ cells, whereas $N_x = 1600$ is used in Test 2 in order to better capture the complex initial profile at the front. In Figure 3, we also see the computed multilayer solutions with the well-balanced numerical scheme proposed here at final time $t = 10$ s. In Figure 3b we also see the solution obtained for the scheme without a well-balanced correction (NoWB hereinafter) of the numerical diffusion (26). Despite the fact that the velocity is quite small, we see the movement of the free surface, that continues moving along the time. This is due to the numerical diffusion and not to a physical effect. On the contrary, we do not see differences between the initial and final states when the well-balanced discretization is used.

Let us remark again that in our scheme, it is not possible to obtain zero velocity or null variation of height because of the regularization method. It always leads to a residual velocity. Interestingly, there is a strong relation between the order of magnitude of both the residual velocity and the variation of the height with respect to the initial condition, and the value of the regularization parameter ($\delta$) used in (13). We show in Figure 4 the maximum velocity, which corresponds to the closest velocity to the free surface, for different values of $\delta$ and for both tests. We see that the velocity is about $10^{-6}$ for $\delta = 10^{-5}$ whereas its order of magnitude is approximately $10^{-16}$ for $\delta = 10^{-15}$. For Test 2 (Figure 4b, inset figure) we see a small spike of order $10^{-10}$ just at the dry front position. It is seen for $\delta = 10^{-10}, 10^{-15}$ but not for $\delta = 10^{-5}$ since the residual velocity is greater than this spurious spike.

In Table 2, we give the errors on the height and the full velocity field for both tests and all the tested values of $\delta$, and also for the NoWB discretization. Here errors are measured against the initial condition, which is the analytical steady state. We see that for all the cases the WB discretization allows us to preserve the steady state up to the precision of the
regularization parameter $\delta$. We also see that the errors are greater for Test 2 for $\delta = 10^{-10}$ and $10^{-15}$ due to the spike at the front position, whereas they are similar for $\delta = 10^{-5}$. We also show the errors for the NoWB discretization. Remark that the variations on the height in that case are approximately of the order of $h(x)$. We must also mention that it is possible to use very small values of $\delta$ in this test because it is a steady-at-rest case. In subsections 5.1 and 5.3 corresponding to transient problems, we cannot use too small values (e.g. $\delta = 10^{-12}$) for stability and convergence reasons (see [42]).

## 5.3 Viscoplastic dam break

The goal of this subsection is to show the advantages of using the multilayer approach for dam break problems. Mainly, we are able to get information about the vertical structure of the fluid at different stages of the flow, in contrast to reconstructing the velocity profile from a (quasi)uniform analytical configuration (see e.g. [7, 32, 19]). We also compare here with the laboratory experiments of [3, 4], and with the results of the lubrication model in these references. Lubrication models are deduced by neglecting convective and time derivative terms in the horizontal momentum equation, where the analytical profile (30d) is obtained. Next, this velocity profile is integrated to obtain the depth-averaged velocity $\bar{u}$, which is introduced in the mass equation

$$ \partial_t h + \partial_x (h\bar{u}) = 0 $$

to get a parabolic PDE describing the free surface evolution. We also compare with experimental vertical profiles of velocity of [4]. Finally, some numerical aspects related to the well-balancing and the space discretization order of the scheme will be also presented.

We consider the dam break test presented in [3], which consists of experiments performed with Carbopol ultrez 10 at a mass concentration of 0.15%, as commented on previously. The material properties are those of (29). They put 6 liters of material in a deposit above a platform, which is inclined to a slope of $\theta$ degrees. Then, the gate is quickly opened and the material flows out. In order to simulate these experiments, we consider $x \in [0, 3.5]$ with $N_x = 1400$ cells and 16 vertical layers. In order to be more accurate in the vertical direction, we used 32 layers in some particular figures showing vertical effects. If the previous default case is not used, this will be clarified in the figure legend. The initial height is

$$ h(x) = \begin{cases} 0.12 + (x - 0.25) \tan \theta & \text{if } x \le 0.5, \\ 0 & \text{otherwise,} \end{cases} $$

the material is at rest at $t = 0$ and the local bottom is $b = 0$. A wall boundary condition is used at $x = 0$ and open boundary at $x = 3.5$ m. Finally, in order to compare with these laboratory experiments and the results of the one-layer model, where the Coussot's formula is used as in [3], we consider the friction condition given by (17), whose results are similar to the Coussot's friction case (16). Results for different friction conditions will be commented below (see Figure 8).

Figure 5 shows the time evolution of the free surface with $\theta = 15^\circ$ and $25^\circ$ for the multilayer model (10) and model (21) with normal stress contributions, and the one-layer model, where we see an initial “inertial” phase of the dam break ($t \le 2$ approximately), with a fast movement of the material, and after that a phase, which we call “spreading” phase, with a slow movement with a small velocity. The distinction between these two phases is clearly visible also in Figure 7 where the evolution of the front position and its velocity are shown. The results of models (10) and (21) are similar for both slopes, except near the front position, where the diffusive term in (21) produces a rounded shape, specially for the case $\theta = 25^\circ$. We see in Figure 5 how the material continues to flow for a long time. In the first part of the domain we see a plateau effect at large times. It is due to
<figure><img src="image_6.png" /><figcaption>Figure 5: Time evolution of the local free surface, for $\theta = 15^\circ$ and $25^\circ$ at different times, with the multilayer model (10) (red dash-dotted lines) and model (21) (orange solid lines), and the one-layer model (cyan dashed lines).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Upper figure: local free surface $h(x)$ at different times; Lower figures: time evolution of the local free surface at point $x = 0$ m, for $\theta = 15^\circ$ with the multilayer model (10) (red dash-dotted lines) and the lubrication model (green (+) lines).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Time evolution of the front position $x_f$ (left-hand side) and the front velocity $v_f$ (right-hand side) in semilog scale, for $\theta = 15^\circ$ and $25^\circ$, with the one-layer (cyan dashed lines) and multilayer model (10) (red dash-dotted lines) and multilayer model (21) (orange solid lines). Gray dotted lines are multilayer (10) solutions by a refinement of the mesh on $x$-direction $\Delta x/2$, and circles are data experiments. Green (+) lines are the solution of the lubrication model. Inset figures in the right-hand side are zooms.</figcaption></figure>

the decreasing of the height, since looking at the sufficient condition to be at steady state (23), when the height is small enough with very small velocity, then the material stops. It leads to a static thin layer of material, whose height is approximately $h_p$ (see (30b)). It is also observed in Figure 6, where we show a comparison with the lubrication model for the free surface, and the time evolution of the height at point $x = 0$ m. We see that with the multilayer model (10) (similarly model (21)), this height decreases fast until the critical height $h_p$ approximately, and from then it is constant. We also see a convergence toward this threshold for the lubrication model at very long times. In particular, it is consistent with the shape of the height in the first part of the domain predicted by the multilayer model, where a layer remains static and not all the material flows over the slope.

We also observe that the flow with the one-layer model is faster than with the multilayer model, as expected, since multilayer model include vertical viscous effects. We also see that the influence of the multilayer approach is more relevant for large slopes, where there are more differences between one-layer and multilayer results. It is also what was observed (though the physics is quite different) for dry granular flows (see [27]).

The approximation of the front position ($x_f$) is shown in Figure 7 for the one-layer, multilayer and lubrication models. Nevertheless, we shall remark that the front position depends on a lot of different factors (wet/dry treatment, non-hydrostatic effects, and other
<figure><img src="image_9.png" /><figcaption>Figure 8: Time evolution of the front position $x_f$ for $\theta = 15^\circ$ and $25^\circ$ with the multilayer model (10) and using different friction conditions: Coussot’s formula (16) (brown dashed lines), no slip (18) (blue dash-dotted lines) and friction condition (15) for several values of $K_b$ ($K_b = 1.93K$ corresponds to (17), purple squares). Circles are data experiments.</figcaption></figure>

assumptions in the simplified model). From a qualitative analysis of Figure 7, we extract
several conclusions. First, as commented before, multilayer effects are higher for larger
slopes. Second, the front positions with multilayer and one-layer models remain very close
to each other, especially in the inertial phase. We see small differences between them in
the transition to the spreading phase, where the multilayer models predict a shorter front
position. However, it is not enough to recover the experimental results. In particular,
in the first inertial phase of the flow the results of shallow models are not far from the
experimental ones, contrary to the lubrication model. However, the second spreading phase
is not properly reproduced neither by multilayer nor, by one-layer and lubrication models.
It is in agreement with the conclusions in [3]. We also show in Figure 7 the time evolution
of the front velocity ($v_f$). We see an agreement between the experiments and lubrication
model at the last stage of the flow. However, it is far from the experiments at the first
inertial phase. As a consequence, both shallow and multilayer models reproduce better
the global behavior of the experiments, taking into account both phases. In the spreading
phase, looking at the experiments, the front moves with a slightly decreasing velocity. This
decrease of the velocity is smaller with shallow models. In order to reproduce this effect
at long times, it is very important to refine the horizontal discretization, as suggested
by the lines representing the solution with $\Delta x/2$ in the evolution of both the front and
the front velocity. One could conjecture that this effect can also be achieved by adding
a fully non-hydrostatic pressure involving all the contributions of the stress tensor. This
is a difficult task from the numerical point of view, even for the Newtonian case, which
would be interesting to solve in the future. However, this is out of the scope of present
paper. Moreover, friction with lateral walls could also have a significant impact on the
approximation of the front position. Adding this effect in simulations, where a calibration
of friction parameters is needed, should help to reduce discrepancies with the laboratory
experiments, as also explained in [3].

We also show in Figure 8 the evolution of the front position when using different
friction coefficients defined in subsection 3.2, for $\theta = 15^\circ$ and $25^\circ$. We see that in the
initial inertial phase there is almost no influence. However, some differences appear in the
spreading phase. In particular, we see that friction condition (17) produces similar results
to the Coussot’s formula (16), which is expected since we are using a value of $n$ very close
to 0.3, that is the one in the Coussot’s formula. We also see that increasing significantly
the value of $K_b$ ($K_b = 10K$) has not a big influence on the front position. It agrees with
<figure><img src="image_10.png" /><figcaption>Figure 9: Velocity field (left-hand side), and $\log_{10}(|\partial_z u|)$ (right-hand side) for $\theta = 15^\circ$ at times $t = 0.2, 0.5, 1.0, 1.5$ s during the inertial phase of the flow, computed with model (10). We use 32 vertical layers here.</figcaption></figure>

our previous conjecture: adding other 3D effects would be more relevant to improve the approximation of the front position.

Although the front position is not significantly improved, the power of the multilayer approach is the fact that it recovers the vertical information of the flow. Figures 9 and 10 show, for model (10), the velocity field and also the distribution of $||D|| \approx |\partial_z u|$ for the case $\theta = 15^\circ$. Figure 9 focuses on the inertial phase whereas Figure 10 shows a specific time of the spreading phase. We see that the multilayer approach reproduces the expected profile of velocity, which has a bottom sheared layer ($||D|| > 0$) and top unsheared layer with a
<figure><img src="image_11.png" /><figcaption>Figure 10: Velocity field (upper figure), and $\log_{10}(|\partial_z u|)$ (lower figure) for $\theta = 15^\circ$ at time $t = 10$ s during the spreading phase of the flow, computed with model (10). We use 32 vertical layers here.</figcaption></figure>

nearly constant profile ($u(z)$) constant, $||D|| \approx 0$). This divides the flow in a lower sheared and upper unsheared layers, in agreement with conclusions in [4] for these viscoplastic experiments. We remark that this division does not necessarily appear very close to the front position, as it will be discussed later. Actually, a strength of the multilayer model is the ability to approximate this interface between the sheared and pseudoplug layers, that is clearly discerned in these figures when showing $\log_{10}(||D||)$. Remark that, as it has been emphasized in previous subsections, we do not recover machine-precision zero velocity nor zero $||D||$ due to the regularization of the viscosity coefficient. In practice, we use a tolerance $4\delta$ to compute this interface, which is of the order of magnitude of the regularization parameter.

One could suggest that a way to recover the vertical profile of velocity is to reconstruct it from the analytical expression (30d) with (31), which takes into account the term $\partial_x h$. In fact, the interface between the sheared and unsheared layers is also approximated once one has $h_p$ (30b). In Figure 11, we show the velocity profiles at several times and locations obtained by the multilayer model (10) for $\theta = 15^\circ$, together with the profiles obtained by the analytical formula once $h$ (and therefore $\partial_x h$) is known. We remark that the analytical approach is computed from the multilayer height (both $h$ and $\partial_x h$). Actually, we show the analytical profiles with and without the correction corresponding to $\partial_x h$ (31). It can be seen that the analytical profiles differ largely from the multilayer profiles during the inertial phase (upper figure). However, during the spreading phase, the approximations given by the analytical formula are in good agreement with the multilayer profiles. It indicates that in the spreading phase, the flow is very close to a (quasi)uniform regime, as expected. It is interesting to see the influence of the $\partial_x h$ term in the analytical profile. It notably improves the approximation, in particular near the front. Moreover, this correction goes in the good direction with respect to the agreement with the multilayer profiles in both the first and last part of the domain. We also show the analytical interface between the sheared/unsheared zones (including $\partial_x h$). Similar conclusions are obtained, i.e., at initial times the results differ but they are very close for the spreading phase.
<figure><img src="image_12.png" /><figcaption>Figure 11: Profiles of velocity at $x = 0.25, 0.5, ..., 2$ m at times $t = 0.5$ s (inertial phase) and $t = 6, 10$ s (spreading phase). We use 32 vertical layers here. Black lines are the multilayer profiles, cyan (+) lines are the analytical profile (30) whereas circled red profiles are analytical profiles with the influence of $\partial_x h$ (31), both analytical approaches computed with the multilayer height. For a better graphical representation each 0.25 m, all profiles are normalized with $U = 4 \max_{i,\alpha} |u_{\alpha,i}|$ for $i = 1, ..., N_x, \alpha = 1, ..., N$. Blue dashed lines represent the interface between the sheared and unsheared layers, and red dotted lines are its analytical approximations given by $h_c$ with (31).</figcaption></figure>

In Figure 12, we quantify the relative differences made by the analytical and the multilayer velocity fields and interfaces, for $\theta = 15^\circ, 25^\circ$. We see here that these differences are huge in the inertial phase. However, these differences stabilize around 2% for $\theta = 15^\circ$ (5% for $\theta = 25^\circ$) for $L_1$ and $L_2$ norms, and 7% for $\theta = 15^\circ$ (12% for $\theta = 25^\circ$) for $L_\infty$ norm, for the interface approximation. Concerning the velocity field, these differences are greater. They stabilize around 11% for $\theta = 15^\circ$ (9% for $\theta = 25^\circ$) for $L_1, L_2$ norms, and 100% for $\theta = 15^\circ, 25^\circ$ for $L_\infty$ norm. We note that similar results were obtained here when using $\delta = 10^{-6}, 10^{-7}$, except for the initial inertial phase where errors are slightly larger. It is due to some spurious oscillations that appear near the front position when too small values of the regularization parameter are used. Note that the tolerance for recovering the yield surface is also decreased with $\delta$ since we use the value $4\delta$ as tolerance.

Figures 13 and 14 show a qualitative comparison, using model (10), of the velocity
<figure><img src="image_13.png" /><figcaption>Figure 12: Relative differences in semilog scale between the computed velocity field and interface ($h_c$) between the sheared and unsheared layers, and the approximations given by the analytical expressions (30d), (30e) and (31) ($\partial_x h$ computed with the multilayer height). Left: slope $\theta = 15^\circ$. Right: slope $\theta = 25^\circ$. We use 32 vertical layers here.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Vertical profiles of velocity at $x = x_f + \Delta x_f$ for $\theta = 15^\circ$. Blue dash-dotted lines are the multilayer results (model (10)), red circles (resp. green crosses) are the analytical profiles (31) with the gradient $\partial_x h$ computed from the multilayer (resp. experimental) height profile, and gray circles are experimental data of [4]. We use 32 vertical layers here.</figcaption></figure>

profiles measured in the experiments in [4] for $\theta = 15^\circ, 25^\circ$. Let us detail this comparison. In the experiments, the velocity profile is always measured at the observation point $x = 2.55 \text{ m}$. Now, they show different profiles of velocity depending on $\Delta x_f = x - x_f$, that is, the distance from the point where the profile is measured to the front position. Notice that the values of $\Delta x_f$ are small, so they gave the velocity profile very close to the front position. Let us remark that it is a difficult task from the numerical point of view since we need to deal with a dry front, where $\partial_x h$ goes to infinity when approximating in the neighborhood of the front position. Figure 13 shows the velocity profile for $\theta = 15^\circ$ measured when $\Delta x_f = -30.6 \text{ mm}$. We show this profile measured exactly at $x = 2.55 \text{ m}$, and also when we move 1 and 2 cm to the right but keeping the same distance to the front $\Delta x_f$. Notice than each subfigure then corresponds to a different time. We see that the velocity profiles obtained with the multilayer model are in good agreement with the experiments. The analytical solution, where $h$ and $\partial_x h$ are computed from the multilayer results, produces in this case very small (almost zero) velocity, most surely because the shape of the front computed with the multilayer model is not steep enough. We also see in
this figure how the analytical approach with $h$ and $\partial_x h$ corresponding to the experimental results, where $\partial_x h$ is computed by an interpolation technique (see [4] for details), is able to recover the experimental velocity profile.

<figure><img src="image_15.png" /><figcaption>Figure 14: Vertical profiles of velocity at $\Delta x_f = x - x_f$ for $\theta = 25^\circ$ at time $t = 14.4$ s. Blue dash-dotted lines are the multilayer results (model (10)), red circles (resp. green crosses) are the analytical profiles (31) with the gradient $\partial_x h$ computed from the multilayer (resp. experimental) height profile, and gray circles and circle-dotted line are experimental data of [4]. The blue solid line is the multilayer free surface approximation and the blue dashed line is the interface between the pseudoplug layer and the sheared layer. Finally, the gray circle-dotted line is the experimental free surface. We use 32 vertical layers here.</figcaption></figure>

In Figure 14, we show the vertical profiles for a fixed time, and for several distances to the front position $\Delta x_f$. In this case, we cannot compare when the front is close to $x = 2.55$ m because the multilayer model is in a different flow regime than the one in the experiments when the material reaches this location. Let us explain this fact. When the velocity profiles are measured in the experimental case, the flow is in the spreading phase for more than 10 seconds, which is not the case for the multilayer case, where the duration of the inertial phase is much longer. Then, in order to compare the observed velocity profiles, we should place ourselves in the same regime, i.e. when the flow is in the spreading phase since a time interval similar to that of the experiments. Here, we thus show the vertical profiles for several distances $\Delta x_f$ at time $t = 14.4$ s, which allows us to make a qualitative comparison of the obtained vertical profiles. In general, we see a good agreement between the profiles far from the front ($|\Delta x_f| \ge 40$ mm). However, the velocity very close to the front ($|\Delta x_f| \le 20$ mm) is underestimated. For the analytical approach, which is again computed with the multilayer height, we find similar conclusions, except that the analytical approach strongly underestimates the velocity for $|\Delta x_f| \le 72.5$ mm, and the predicted velocity is very small very near the front. It is due to the fact that the height is not well approximated very close to the front position, as we see in this figure. However, when the analytical profile is computed using the experimental height (and also $\partial_x h$), the velocity very close to the front is overestimated (see inset figure in Figure 14 corresponding to a zoom). It is worth mentioning that we do not see a plug zone very close to the front in the experiment but a pseudoplug zone. It is not reproduced neither by model (10) nor the analytical approach.

As explained in subsection 3.2.2, in [8] it was shown that in non-uniform flows we recover a pseudoplug layer, in which the shear rate is $||D(\mathbf{u})|| \sim \varepsilon$. In this section, we have explained how to obtain a model reproducing this pseudoplug layer instead of a plug layer. In order to show how the normal stress contributions allow model (21) to reproduce the pseudoplug layer, we see in Figure 15 the velocity profiles close to the front for the same case as before ($\theta = 25^\circ$). We see how the velocity profiles at $|\Delta x_f| \le 20$ mm have not a plug layer but a pseudoplug layer when they are computed with model (21). Notice that, in this figure, the free surface computed with model (21) presents a small change
of curvature. It is a numerical artifact, which is due to the numerical resolution of the full 2D ($x-z$) linear system close to the front position. We also show in Figure 16 the shear rate $||D_\epsilon||$ (9b), which includes the contribution of the term $\partial_x u$ in this case, where the shear rate is of the order of the shallowness parameter $\epsilon = H/L \approx 0.02/2 = 10^{-2}$ in the pseudoplug layer, as shown in [8]. In these figures, we also observe that the interface separating the pseudoplug and sheared layers changes its convexity near the front, which is in agreement with [19].

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

Figure 15: Comparison of vertical profiles of velocity close to the front position for $\theta = 25^\circ$ computed with multilayer models (10) (blue lines) and (21) (brown lines). Thin dashed lines are the interfaces separating the (pseudo-)plug/sheared layers. In the case of model (21), we use the criterion $||D_\epsilon|| < 10^{-2} \approx \epsilon = H/L$. Gray circles and circle-dotted line are experimental data of [4]. Blue (resp. orange) solid and dashed lines represent the free surface and the interface between the pseudoplug layer and the sheared layer computed with the multilayer model (10) (resp. (21)) and the gray circle-dotted line is the experimental free surface. We use 32 vertical layers here.

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

Figure 16: $\log_{10}(||D_\epsilon||)$ for $\theta = 25^\circ$ during the spreading phase of the flow ($t = 14.4$ s) for model (21). We use 32 vertical layers here.

Let us remark that an added-value of the multilayer approach presented here, with respect to the lubrication theory, is the fact that it is possible to recover the velocity field also in the first inertial phase (where the lubrication theory is not valid) and not only in the last spreading phase.

Figure 17 shows the importance of the discretizations (in time and space) for this dam break problem. In Figure 17a we see the effect of the numerical diffusion on the spreading phase when using the first or second order WB schemes. Recall that the extension to second order in space and time is done with a MUSCL reconstruction and the Heun
method. During the inertial phase, the differences are rather small, while they become
significant in the spreading phase. Actually, those become obvious in the formation of the
static layer (at the back of the flow) whose height is approximately $h_p$. It has a straight
shape with the second order scheme while it is smoothed with the first order scheme, which
is a usual effect of numerical diffusion. Figure 17b shows the results of the proposed scheme
and the scheme without a well-balanced treatment (NoWB). As expected, at initial times
(before $t = 2$ s approximately) there are not significant differences. However, the dynamics
are vastly different at larger times. We see the effect of the numerical diffusion in both the
left part (close to $x = 0$) and the front. At the advancing front, the well-balanced scheme
exhibit an almost straight shape, while an advanced lower layer of material is formed with
the NoWB scheme, which has no physical meaning. In addition, in the left part, we see
the plateau effect with (approximately) height $h_p$ for the well-balanced scheme, whereas
all the material flows with the NoWB scheme, contrary to what is predicted by the steady
state condition (see (23)).

Concerning the computational efforts of multilayer model (10), Table 3 shows the computation time needed to reach the final time $t_f = 15$ s, with the slope $\theta = 25^\circ$, $N_x = 1400$ and $CFL = 0.5$. We also show the ratio with respect to the one-layer case, and with respect to the case of using half of layers. As it is highlighted in subsection 5.1, the computation time increases linearly with the number of layers. In particular, the multilayer system with 16 layers is 8.5 times more expensive than the one-layer model. We remark again that solving model (21) requires a huge computational effort. For instance, for the same case ($\theta = 25^\circ$, $N_x = 1400$ and 16 layers), it takes 24.3 hours, thus, it is 332 times more expensive than model (10).

<figure><img src="image_18.png" /><figcaption>Figure 17: Effect of the numerical diffusion in the spreading phase of the flow (a) with a first/second order time-space scheme (b) without the well-balanced treatment of (26). Model (10) is used here. Inset figures are zooms of the domain.</figcaption></figure>
<table><caption>Table 3: Computation times (s) needed by model (10) to reach $t_f = 15$ s, and speed-ups, for the case $\theta = 25^{\circ}$, $N_x = 1400$ and $CFL = 0.5$.</caption><thead><tr><th>N. layers (N)</th><th>N = 1</th><th>N = 2</th><th>N = 4</th><th>N = 8</th><th>N = 16</th><th>N = 32</th></tr></thead><tbody><tr><td>t<sub>comp</sub>(s)</td><td>30.9</td><td>45.7</td><td>73.8</td><td>137.8</td><td>262.6</td><td>544.4</td></tr><tr><td>t<sub>comp</sub>(N)/t<sub>comp</sub>(N/2)</td><td>-</td><td>1.5</td><td>1.61</td><td>1.87</td><td>1.91</td><td>2.07</td></tr><tr><td>t<sub>comp</sub>(N)/t<sub>comp</sub>(1)</td><td>-</td><td>1.5</td><td>2.4</td><td>4.5</td><td>8.5</td><td>17.6</td></tr></tbody></table>

# 6 Conclusions

Two multilayer Herschel-Bulkley models have been proposed. Their derivation stems from the Navier-Stokes system through an asymptotic analysis and the multilayer approach, where the vertical profile is no more constant along the normal direction but piecewise constant. These models differ in the viscous terms coming from the normal stress components, which allow the model to reproduce a pseudoplug layer instead of a true plug zone. For such multilayer models, a well-balanced finite volume/finite difference scheme has been proposed, allowing us to preserve a relevant family of steady solutions at rest. In addition, we deal with the singularity in the definition of the stress tensor by means of a regularization technique, whose influence has been numerically studied for steady solutions.

In the numerical tests, we have compared the numerical results with an analytical solution for uniform flows, showing the ability of the multilayer model to capture the vertical profile of velocity observed for these fluids. It exhibits a sheared (lower) layer and a (pseudo-)plug (upper) layer. We have shown the ability of the scheme to preserve steady solutions at rest, as well as the dependence of the obtained errors on the regularization parameter ($\delta$). Moreover, we have also shown the importance of the hydrostatic reconstruction of the well-balanced treatment for the numerical scheme. A comparison with viscoplastic dam break experiments of [3, 4] has also been performed. In particular, we show a comparison with the lubrication theory, we measure vertical profiles of velocity, and we demonstrate the ability of the proposed model to reproduce the sheared (bottom) and (pseudo-)plug (top) layers in general flows. For these flows, we distinguish an initial inertial phase and a (late) spreading phase. In addition, the accuracy of the analytic approach during the second quasi-uniform regime has also been studied. Furthermore, the pseudoplug layer has also been recovered for a multilayer model accounting for the normal stress contributions. In practice, both multilayer models, (10) and (21), produce similar results except for the pseudoplug zone and the velocity profiles very close to the front, whereas the computational cost is dramatically increased for model (21). In the future, it would be interesting to design efficient duality methods for multilayer systems, as well as going towards fully non-hydrostatic models including the contributions of the stress tensor to the pressure profile.

# Acknowledgments

We would like to warmly thank Christophe Ancey for kindly providing us with the data of the experimental results from [4, 3].

This research has been partially supported by the Spanish Government and FEDER through the research project RTI2018-096064-B-C22. This research has also been conducted with financial support from the French National Research Agency (ANR) through
the research project VPFlows, ANR-20-CE46-0006, as well as the CNRS through the
research project Visco3plug, of the Infiniti program (2017-2018).

References

[1] C. Acary-Robert, E. D. Fernández-Nieto, G. Narbona-Reina, and P. Vigneaux. A well-balanced finite volume-augmented Lagrangian method for an integrated Herschel-Bulkley model. *Journal of Scientific Computing*, 53(3):608-641, apr 2012. 3

[2] C. Ancey. Plasticity and geophysical flows: A review. *Journal of Non-Newtonian Fluid Mechanics*, 142:4 – 35, 2007. 3

[3] C. Ancey, N. Andreini, and G. Epely-Chauvin. Viscoplastic dambreak waves: Review of simple computational approaches and comparison with experiments. *Advances in Water Resources*, 48:79–91, nov 2012. 4, 6, 10, 16, 21, 24, 32

[4] N. Andreini, G. Epely-Chauvin, and C. Ancey. Internal dynamics of Newtonian and viscoplastic fluid avalanches down a sloping bed. *Physics of Fluids*, 24(5):053101, may 2012. 4, 16, 21, 26, 28, 29, 30, 32

[5] A. Armanini, M. Larcher, E. Nucci, and M. Dumbser. Submerged granular channel flows driven by gravity. *Advances in Water Resources*, 63:1–10, jan 2014. 3

[6] 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, jan-feb 2011. 3, 13

[7] J. L. Baker, T. Barker, and J. M. N. T. Gray. A two-dimensional depth-averaged $\mu(I)$-rheology for dense granular avalanches. *Journal of Fluid Mechanics*, 787:367–395, 1 2016. 21

[8] N. J. Balmforth and R. V. Craster. A consistent thin-layer theory for Bingham plastics. *Journal of Non-Newtonian Fluid Mechanics*, 84(1):65–81, July 1999. 3, 6, 7, 17, 29, 30

[9] N. J. Balmforth, R. V. Craster, A. C. Rust, and R. Sassi. Viscoplastic flow over an inclined surface. *Journal of Non-Newtonian Fluid Mechanics*, 142(1-3):219–243, mar 2007. 3

[10] N. J. Balmforth, I. A. Frigaard, and G. Ovarlez. Yielding to Stress: Recent Developments in Viscoplastic Fluid Mechanics. *Annual Review of Fluid Mechanics*, 46(1):121–146, 2014. 2

[11] M. Bercovier and M. Engelman. A finite-element method for incompressible non-Newtonian flows. *Journal of Computational Physics*, 36(3):313–326, 1980. 3

[12] A. Bermúdez and C. Moreno. Duality methods for solving variational inequalities. *Computers & Mathematics with Applications*, 7(1):43–58, 1981. 5

[13] E. C. Bingham. *Fluidity and plasticity*. Mc Graw-Hill, 1922. 2

[14] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville. Yield stress materials in soft condensed matter. *Reviews of Modern Physics*, 89(3):035005, August 2017. 2

[15] F. Bouchut. *Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws*. Birkhäuser Basel, 2004. 14

[16] 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. 11
[17] M. J. Castro, T. Morales de Luna, and C. Parés. Well-Balanced Schemes and Path-Conservative Numerical Methods. In Rémi Abgrall and Chi-Wang Shu, editor, *Handbook of Numerical Analysis*, volume 18 of *Handbook of Numerical Methods for Hyperbolic Problems Applied and Modern Issues*, pages 131–175. Elsevier, 2017. DOI: 10.1016/bs.hna.2016.10.002. 16

[18] 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

[19] G. Chambon, P. Freydier, M. Naaim, and J.-P. Vila. Asymptotic expansion of the velocity field within the front of viscoplastic surges: comparison with experiments. *Journal of Fluid Mechanics*, 884, dec 2019. 6, 7, 16, 21, 30

[20] 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. 16

[21] GdR Midi CNRS. On dense granular flows. *The European Physical Journal E*, 14(4):341–365, August 2004. 2

[22] P. Coussot. *Mudflow rheology and dynamics*. A.A. Balkema, Rotterdam - Brookfield, 1997. 6, 10, 16

[23] P. Coussot. Bingham's heritage. *Rheologica Acta*, pages 1–14, December 2016. 2

[24] D. C. Drucker and W. Prager. Soil mechanics and plastic analysis of limit design. *Q. Appl. Math.*, 10(2):157–175, 1952. 2

[25] E. D. Fernández-Nieto, J.M. Gallardo, and P. Vigneaux. Efficient numerical schemes for viscoplastic avalanches. Part 1: The 1D case. *Journal of Computational Physics*, 264:55–90, may 2014. 5

[26] E. D. Fernández-Nieto, J.M. Gallardo, and P. Vigneaux. Efficient numerical schemes for viscoplastic avalanches. Part 2: The 2D case. *Journal of Computational Physics*, 353:460–490, jan 2018. 5

[27] 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. 3, 9, 18, 23

[28] 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 multi-layer discretization. *Journal of Computational Physics*, 356:192–219, 2018. 3, 6, 10, 13, 14

[29] 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. 3, 8, 13

[30] E. D. Fernández-Nieto, P. Noble, and J.-P. Vila. Shallow Water equations for Non-Newtonian fluids. *Journal of Non-Newtonian Fluid Mechanics*, 165(13-14):712–732, jul 2010. 3, 6, 7

[31] M. Fortin and R. Glowinski. *Méthodes de Lagrangien augmenté applications à la résolution numérique de problèmes aux limites*. Méthodes mathématiques de l'informatique 9. Dunod, Paris, 1982. 5

[32] P. Freydier, G. Chambon, and M. Naaim. Experimental characterization of velocity fields within the front of viscoplastic surges down an incline. *Journal of Non-Newtonian Fluid Mechanics*, 240:56–69, feb 2017. 16, 21

[33] I. Frigaard. Simple yield stress fluids. *Current Opinion in Colloid & Interface Science*, 43:80–93, October 2019. 2
[34] 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. 16

[35] 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. 12

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

[37] W. H. Herschel and T. Bulkley. Measurement of consistency as applied to rubber-benzene solutions. *Am. Soc. Test Proc.*, 26(2):621-633, 1926. 2

[38] J. M. Hill, D. Bhattacharya, and W. Wu. Bagnold velocity profile for steady-state dense granular chute flow with base slip. *Rheologica Acta*, jan 2022. 17

[39] X. Huang and M. H. García. A Herschel-Bulkley model for mud flow down a slope. *Journal of Fluid Mechanics*, 374:305-333, nov 1998. 3, 16

[40] P. Jop, Y. Forterre, and O. Pouliquen. A constitutive law for dense granular flows. *Nature*, 441(7094):727-730, June 2006. 2

[41] D. Laigle and P. Coussot. Numerical modeling of mudflows. *Journal of Hydraulic Engineering*, 123(7):617-623, jul 1997. 3

[42] 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. 5, 16, 21

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

[44] J.-M. Piau. Flow of a yield stress fluid in a long domain. Application to flow on an inclined plane. *Journal of Rheology*, 40(4):711-723, jul 1996. 3

[45] M. La Rocca, P. Prestininzi, L. Elango, R. Hinkelmann, and A. Montessori. Depth averaged modelling of loose rectangular granular piles collapsing in water. *Advances in Water Resources*, 143:103663, sep 2020. 3

[46] P. Saramito and A. Wachs. Progress in numerical simulation of yield stress fluid flows. *Rheologica Acta*, pages 1-20, January 2017. 3

[47] T. Schwedoff. Recherches expérimentales sur la cohésion des liquides. *J. Phys. Theor. Appl.*, 8(1):341-359, 1889. 2

[48] S. Yavari-Ramshe and B. Ataie-Ashtiani. Numerical modeling of subaerial and sub-marine landslide-generated tsunami waves-recent advances and future challenges. *Landslides*, 13(6):1325-1368, December 2016. 3