Mul ilaye shallow model o d y g anula lows wi h a weakly
non-hyd os a ic p essu e
C. Escalan e ∗
, J. Ga es-D´ıaz †
, E.D. Fe n´andez-Nie o ‡
, A. Mangeney §
Abs ac
The mul ilaye model p oposed in his pape is a gene aliza ion o he mul ilaye non-hyd os a ic model o shallow
g anula lows Fe n´andez-Nie o e al. [40], he mul ilaye model wi h µ(I) heology Fe n´andez-Nie o e al. [36], and
he monolaye model wi h weakly non-hyd os a ic p essu e o d y g anula lows Ga es-D´ıaz e al.[42]. We show
ha he p oposed model e i ies a dissipa i e ene gy balance. A well-balanced nume ical scheme is p oposed o sol e
he equa ions based on a p ojec ion me hod and a hyd os a ic econs uc ion o he Coulomb ic ion e ms. In o de
o educe he compu a ional cos associa ed wi h sol ing he linea sys em o he p ojec ion me hod, a p ecompu ing
o he ini ial guess o an i e a i e sol e is p oposed. This s a egy allows us o educe he compu a ional ime by
a ound 70% when 20 laye s a e conside ed. In he nume ical es s, we show ha he p oposed model can eco e he
in-dep h eloci y p o iles ypically obse ed in lab expe imen s and cap u e he low/no- low in e ace ha appea s
in g anula a alanches. Du ing he ini ial s age o g anula collapse simula ions, he model is shown o imp o e he
app oxima ion o he mass p o iles compa ed o o he models and o p edic he pa abolic shape o he on eloci y
e olu ion wi h ime, as obse ed in lab expe imen s. In e es ingly, ou nume ical es s show ha he abili y o he
g anula low o o e come obs acles s ongly depends on he model used, which is o s ong in e es o landslide
haza d assessmen .
Keywo ds: Mul ilaye models Non-hyd os a ic p essu e Fini e Volume G anula lows.
1 In oduc ion
Many e o s ha e been de o ed in ecen yea s o he unde s anding o g anula lows, due o hei impo ance in a wide
a ie y o na u al phenomena. Bo h d y and luidized g anula lows ha e been b oadly s udied om a heo e ical poin o
iew, and also many esou ces ha e been employed o pe o m expe imen al s udies, in pa icula a he labo a o y scale
(see e.g. [43,28,48,1] o e iews). The unde s anding and quan i ica ion o he complex beha io o hese lows is s ill a
challenge, since hei dynamics is s ongly condi ioned by many di e en physical p ocesses ( ic ion, dila ancy, luid-g ain
in e ac ions, pa icle seg ega ion, e c.), see o ins ance [65,12,78,6]. In pa icula , he low/no- low local ansi ion ha
con ols e osion/deposi ion p ocesses, in ol es many heo e ical and nume ical di icul ies [55,52,12,63,36].
He e, we deal wi h d y and dense g anula lows, which a e nowadays desc ibed by he µ(I)- heology, in oduced in
[53,54]. This heology is based on a D ucke -P age plas ici y c i e ion, ha de ines he de ia o ic enso as





τ=µ(I)p
∥D∥Di ∥D∥ = 0,
∥τ∥ ≤ µspi ∥D∥= 0,
(1)
whe e pis he p essu e, µsis he angen o he epose angle o he ma e ial, D he s ain- a e enso and µ(I) is a
s ain- a e and p essu e dependen ic ion coe icien , which will be de ined la e . I is known ha he incomp essible
µ(I)- heology is ill-posed o low and high alues o he ine ial numbe I(see [7]). Many wo ks ha e been de o ed o
∗Dp o. Ma em´a ica Aplicada. U. M´alaga, 14071- M´alaga, Spain (escalan[email p o ec ed])
†Dp o. Ma em´a icas. Edi icio Eins ein - U. C´o doba. Campus de Rabanales, 14014-C´o doba, Spain ([email p o ec ed]
‡Dp o. Ma em´a ica Aplicada I. ETS A qui ec u a - U. Se illa. A da. Reina Me cedes S/N, 41012-Se illa, Spain ([email p o ec ed])
§Ins i u de Physique du Globe de Pa is, Seismology eam, U. Pa is-Dide o , So bonne Pa is Ci ´e, 75238, Pa is, F ance. ANGE eam,
CEREMA, INRIA, Lab. J. Louis Lions, 75252, Pa is, F ance. Ins i u Uni e si ai e de F ance (IUF), F ance, ([email p o ec ed])
1
imp o e he well-posed cha ac e o he µ(I)- heology. A egula ized incomp essible µ(I) law, which imp o es he ange
o alues o Iwhe e i is well-posed, was p oposed in [5]. This egula ized exp ession was conside ed in [44,6] among
o he s. A comp essible e sion o he µ(I), called µ(I), ϕ(I)- heology, has also been in es iga ed (see e.g. [17,77]),
conside ing ha he solid olume ac ion ϕis no mo e cons an . Some au ho s sugges ed ha his comp essible e sion
imp o es he well/ill-posed cha ac e o he µ(I) o s eady p oblems (see [8,51]), al hough his is no necessa ily ue
o ansien p oblems (see [80]). As commen ed abo e, when he comp essible heological law is conside ed, he solid
olume ac ion ϕis also a iable (depending on Iin he ae ial case o on he iscous numbe Jin he subme ged case),
and dila ion/con ac ion o he ma e ial occu s (see [77]). Dila ancy e ec s in a shallow dep h-a e aged model o d y
a alanches was ecen ly in es iga ed in [11], whe e he model is ob ained by emo ing he luid phase in he wo-phase
model o deb is lows wi h dila ancy [12]. In addi ion, in he las yea s, some au ho s ha e in oduced non-local e ec s
in he cons i u i e law (see e.g. [76,16]). Howe e , he local µ(I)- heology con inues being e y popula , and ac ually
he mos accep ed law, since many ques ions abou non-local models s ill need o be cla i ied.
In he ollowing, we ocus on he incomp essible µ(I)- heology. This heology has been used in models sol ing he ull
Na ie -S okes equa ions o ee su ace lows. The yield s ess can be deal wi h by de ining a egula ized µ(I)- iscosi y
as in he 2D sol e o [55] o [56]. A di e en app oach was ollowed in [52,63], whe e au ho s used a ini e elemen sol e
combined wi h an Augmen ed Lag angian me hod, also modelling sidewall ic ion e ec s. They showed a quali a i e
good ag eemen wi h g anula collapse labo a o y expe imen s.
Sol ing he ull Na ie -S okes sys em howe e equi es huge compu a ional e o s. Fo his eason, dep h-a e aged
shallow models appea as an in e es ing al e na i e. Many au ho s ha e p oposed such models o g anula lows since he
pionee ing Sa age-Hu e model [79] whe e Moh -Coulomb ic ion e ms wi h cons an ic ion coe icien s a e added o
he classical Sain -Venan sys em. A mo e sophis ica ed ic ion coe icien was p oposed in [74,73,75], which depends on
he hickness and eloci y o he low, e.g. on i s F oude numbe . This a iable ic ion coe icien has been widely used
o model d y g anula lows (see [62,60,72,19] among many o he s). In e es ingly, when hese dep h-a e aged models
a e conside ed in il ed (o local) coo dina es o e a e e ence plane o bo om, he low could s a lowing o slopes
ha could di e (being lowe o highe ) om he angle o epose o he ma e ial. This was shown in [10], whe e au ho s
p oposed a co ec ion o he ic ion coe icien in e ms o a second o de co ec ion o he p essu e. This co ec ion
allows he model o e i y he physical h eshold o low, de ined by he angle o epose.
A igo ous de i a ion o a shallow model wi h he µ(I)- heology unde a dimensional analysis is p esen ed in [46],
whe e au ho s also included a second o de iscous e m in he o m o a ho izon al di usion e m. To de i e he
model, hey needed o assume a Bagnold p o ile o he downslope eloci y e en hough eloci y p o iles a e known o
a y signi ican ly om linea , Bagnold o S-shape p o iles du ing he di e en phases o he low (ini ial accele a ion,
de eloped low, decele a ion) (see [43]). The main d awback o such dep h-a e aged models is he ac ha all he
in-dep h in o ma ion is los due o he shallow hypo hesis and he dep h-in eg a ion p ocedu e. In pa icula , i makes
i impossible o ep oduce he lowe s a ic laye ha exis s in g anula collapse o du ing he s opping phase o g anula
lows. This low/no- low (also called lowing/s a ic) in e ace is shown o depend on he de i a i es o he downslope
eloci y wi h espec o he di ec ion no mal o he opog aphy, an in o ma ion in insically missing in dep h-a e aged
models [57].
As an in e media e s ep be ween ull Na ie -S okes and dep h-a e aged models, mul ilaye (o laye -a e aged) models
[3,39] ha e been shown o be able o eco e impo an in-dep h a ia ions o he low p ope ies such as eloci y p o iles.
In [39] a mul ilaye model is deduced as he se o equa ions e i ied by a pa icula weak solu ion o he ull Na ie -S okes
sys em. The e, he domain is subdi ided in shallow laye s and he eloci y ield is app oxima ed by a laye wise cons an
unc ion. Such a model is shown o ep oduce he di e en shapes o he eloci y p o iles associa ed wi h di e en low
egimes. Fo example, in [36], simula ions pe o med wi h a mul ilaye model o d y g anula lows wi h he µ(I)- heology
we e success ully compa ed o g anula collapse expe imen s in [61], showing an excellen quali a i e ag eemen conce ning
he shape o he deposi s, o he eloci y p o iles, o he in luence o an e odible bed on he a el dis ance ( unou ) o
hese lows. Mo eo e , his mul ilaye model was able o app oxima e he low/no- low in e ace sepa a ing he mobile
and s a ic laye s o ma e ial du ing g anula column collapse. The mul ilaye app oach makes i also possible o p ope ly
accoun o sidewalls ic ion ha has a s ong in luence on he low dynamics [37].
All hese dep h-a e aged models assume ha he p essu e is hyd os a ic. Howe e , du ing he ini ial des abiliza ion,
he mass geome y may signi ican ly di e om a shallow laye , hus inducing s ong non-hyd os a ic e ec s a ising om
he non-negligible accele a ion in he di ec ion no mal o he opog aphy. This could explain pa o he disc epancies
obse ed when compa ing he esul s o shallow dep h-a e aged models wi h g anula collapse expe imen s, in pa icula
du ing he ini ial ins an s e.g.[36]. Going beyond he hyd os a ic assump ion, dispe si e models ha e become a e y
ac i e opic o esea ch in ecen yea s, looking o an imp o emen o nonlinea dispe si e p ope ies.
Two amilies o dispe si e models can be conside ed: Boussinesq ype sys ems, whe e high o de de i a i es o he
a iables a e included in he model (see, e.g., [15,67] among o he s); and non-hyd os a ic sys ems (see, e.g., [23,82]),
whe e he p essu e is spli in o hyd os a ic and non-hyd os a ic coun e pa s, in ol ing new unknowns and es ic ions in
he model. Ac ually, mos classical dispe si e models o wa e lows can be e o mula ed as non-hyd os a ic sys ems, as
2
shown in [32]. The ad an age o he second app oach is ha only i s -o de de i a i es appea in he model, which a e
easie o ea nume ically. In addi ion, he speci ic s uc u e o hese models and hei simila i ies o he Shallow Wa e
Equa ions (SWE) enables he ex ension o nume ous amilia nume ical schemes o SWE o non-hyd os a ic models, as
demons a ed in [30]. One o he d awbacks o non-hyd os a ic models is ha hey do no all in o he class o hype bolic
PDE sys ems. Fu he mo e, he e is a need o add ess implici p oblems o include he non-hyd os a ic p essu e e ms
o a oid he es ic i e ime s ep ha would gi e an explici scheme. Va ious nume ical me hods can be ound in he
li e a u e o app oxima ing hype bolic-ellip ic dispe si e sys ems: he p essu e-co ec ion-based me hods (e.g., [26,81]
in he amewo k o Na ie -S okes equa ions), semi-implici me hods on s agge ed meshes (see, e.g., [4,24,50,82]), o
mo e ecen ends based on ob aining a hype bolic elaxa ion sys em ha can be disc e ized wi h explici me hods (see
[20,27,32,34,47,49]).
Howe e , such non-hyd os a ic models widely de eloped o wa e wa e p opaga ion, ha e been poo ly applied o
g anula lows. The non-hyd os a ic con ibu ion associa ed o he accele a ion in he di ec ion no mal o he slope
desc ibed abo e is di e en om o he non-hyd os a ic e ec s such as hose a ising om he iscous s ess enso in ol ed
in he p essu e [13], which a e always igno ed in g anula low models owing o hei complexi y.
Recen ly, a shallow model o d y g anula lows wi h a weakly non-hyd os a ic p essu e was in oduced in [42]. The
p essu e is called weakly non-hyd os a ic because only he con ibu ion o he no mal accele a ion is conside ed in he
p essu e and no he ones coming om he s ess enso . This is a s ong limi a ion since hese las e ms a e o i s
o de in he dimensional analysis based on he shallow pa ame e εwhile hose associa ed o he no mal eloci y a e o
second o de . This choice elies only on he di icul y o handle he s ess con ibu ion om a nume ical poin o iew.
Despi e his limi a ion, e y p omising esul s we e ob ained. In pa icula , i was possible o eco e he pa abolic shape
o he e olu ion o he on eloci y wi h ime while hyd os a ic models ailed o ep oduce he ini ial on accele a ion.
In ha wo k, he in luence o he coo dina e sys em (local o Ca esian) was s udied, showing ha he non-hyd os a ic
model in Ca esian coo dina es p oduces oughly simila deposi s han he hyd os a ic model in local coo dina es. This
was pa ially sugges ed p e iously by [29], whe e hey conside ed a Ca esian model wi h a co ec ion ela ed o he
e ical accele a ion.
Conce ning non-hyd os a ic models o mul ilaye sys ems, a hie a chy o such models o Eule equa ions was
in oduced in [40], all o hem e i ying a dissipa i e ene gy balance. They a e deno ed as LDNHkmodels, whe e k
is he deg ee o app oxima ion o he e ical eloci y. This wo k p esen s he ex ension o he LDNH0model o g anula
lows wi h he µ(I)- heology, ob aining a mul ilaye model wi h weakly non-hyd os a ic p essu e. I can be seen as he
ex ension o he shallow model in [42] o he mul ilaye amewo k, bu also as he ex ension o he mul ilaye model
[36] o he non-hyd os a ic amewo k. In pa icula , a e he p omising esul s o he one-laye model wi h a weakly
non-hyd os a ic p essu e [42], i is a na u al ques ion whe he i s mul ilaye e sion could imp o e he modelling o
g anula lows, adding he ad an ages o imp o ing he in-dep h desc ip ion o he low shown in [36]. Ou goal is o
p esen his weakly non-hyd os a ic mul ilaye model. I is expec ed ha his model eco e s he good esul s o bo h
p e ious models [36,42], and sol es some o he d awbacks o each o hem, being cheape compu a ionally han a ully
non-hyd os a ic model.
The mul ilaye weakly non-hyd os a ic µ(I)-model p oposed in he pape e i ies a dissipa i e ene gy balance. A
ini e olume disc e iza ion is also in oduced, being well-balanced o s eady solu ions a es wi h non-cons an ee
su ace. As expec ed, he p ize o pay is he inc ease o he compu a ional e o , mainly due o he need o sol e he
non-hyd os a ic p essu e in mul ilaye sys ems, which is made ia an i e a i e linea sys em sol e . Conce ning ha , we
p opose a simple s a egy ha allows us o no ably educe he compu a ional cos .
This pape is o ganized as ollows. In Sec ion 2 he ini ial sys em is s a ed and he dimensional analysis is pe o med.
Sec ion 3is de o ed o he mul ilaye app oach o he sys em, p esen ing he inal model. The nume ical app oxima ion
is desc ibed in Sec ion 4and some nume ical es s a e p esen ed in Sec ion 5. Conc e ely, we s udy he in luence o
he coo dina e sys em in he mul ilaye case, and compa e wi h labo a o y-scale g anula collapse expe imen s. We also
pe o m con e gence and well-balanced es s, and show he e iciency o he p oposed s a egy o speed up he i e a i e
sol e o he linea sys em associa ed o he p essu e unknowns. Finally, some conclusions a e p esen ed in Sec ion 6.
2 De i a ion o a non-hyd os a ic µ(I)-model
In his sec ion we de i e he non-hyd os a ic mul ilaye model ha we p opose. I is based on he dimensional analysis
made in [42], and he non-hyd os a ic ex ension o mul ilaye sys ems in oduced in [40]. Fi s ly, we w i e he go e ning
equa ions and he non-dimensional 2D Na ie -S okes sys em, unde he assump ion o a weakly non-hyd os a ic p essu e.
3
2.1 Go e ning equa ions
Le us b ie ly summa ize he go e ning equa ions desc ibing hese lows. We conside a 2D g anula mass wi h cons an
densi y ρ∈Rand eloci y u∈R2, whose dynamics is desc ibed by he 2D Na ie -S okes sys em
(∇·u= 0,
ρ(∂ u+ (u·∇)u) = ∇·σ+ρg,
whe e σ=−pI+τis he o al s ess enso , wi h p∈R he o al (hyd os a ic and non-hyd os a ic) p essu e and τ he
de ia o ic enso . This enso is de ined by (1), which can be w i en in he egime when he ma e ial lows (∥D∥ = 0),
as
τ=ηD(u),
whe e η∈Ris he iscosi y coe icien and he s ain- a e enso D(u) = 1
2∇u+ (∇u)′. In he µ(I)- heology,
in oduced in [53], he iscosi y is de ined as (see de ails in [55,36] among o he s)
η=µ(I)p
p||D(u)||2+δ2,
wi h δ > 0 a egula iza ion pa ame e (see e.g. [45,66]), ||D|| =√0.5D:Dand µ(I) he ic ion coe icien de ined by
µ(I) = µs+µ2−µs
I0+II,
whe e I0and µ2> µsa e cons an pa ame e s depending on he ma e ial, and
I=2ds||D(u)||
pp/ρs
is he ine ial numbe , wi h dsand ρs he g ain diame e and densi y, espec i ely. The g ain densi y and he appa en
low densi y (ρ) a e ela ed by ρ=φsρswi h φs he solid olume ac ion.
As bounda y condi ions, simila ly o wha is done in [36,42], we conside :
•A he ee su ace: he usual kinema ic condi ion and no mal s ess balance
N +u·nb+h= 0,and p= 0,
wi h N ,nb+h he downwa d ime-space no mal ec o o he ee su ace.
•A he bo om: he non-pene a ion condi ion and a Coulomb ic ion law
u·nb= 0,and σ nb−σ nb·nbnb=−µ(I)pu
|u|,0′
,
wi h nb he downwa d space ec o no mal o he bo om and u he i s componen o u.
2.2 Ini ial sys em and dimensional analysis
Le us s a by es ablishing he no a ion used he e. We conside local (o il ed) coo dina es, as usually done in g anula
low models, ha e e o a ixed inclined plane (a s aigh line in he 1D x−zcase). Conc e ely, we conside a plane gi en
by eb(x)=(xend −x) an θ, wi h (x, z) he Ca esian coo dina es, xend he inal poin o he domain, and θ > 0 a posi i e
slope angle. We adop he e he usual con en ion in geophysical lows, ha conside s posi i e angles θ o nega i e slopes.
Le (X, Z) be he local coo dina es, measu ed in he di ec ion along and no mal o he inclined plane eb, espec i ely,
and u= (u, w) he eloci y ec o , whose componen s a e he downslope and no mal eloci ies. We shall also conside a
local bo om, b(X), as he de ia ion om he e e ence plane, and h(X) he dis ance om b(X) o he ee su ace, bo h
measu ed in he no mal di ec ion (see Figu e 1).
4
Figu e 1: Ske ch o he Ca esian (blue) and local ( ed) e e ence sys em.
The dimensional analysis pe o med he e is he same as in [42] (see also [36,46] among o he s). Le us emind i
o he pu pose o comple eness. De ining H, L and U he cha ac e is ic heigh , leng h and eloci y, we de ine he usual
shallowness pa ame e ε=H/L, and deno ing wi h ildes (e·) he non-dimensional a iables, we ob ain
(X, Z, )=(Le
X, H e
Z, (L/U)e ),
(u, w)=(Ueu, εU ew),
h=Heh, ρ =ρ0eρ, p =ρ0U2ep,
(τXX , τXZ, τZZ ) = ρ0U2(εgτXX ,gτXZ , εgτZZ),
wi h
gτXX =eη∂ e
Xeu, gτXZ =eη
2∂e
Zeu+ε2∂e
Xew,gτZZ =eη∂e
Zew,
and he e o e
e
D(e
u) = U
H
1
2
2ε∂ e
Xeu ∂e
Zeu+ε2∂e
Xew
∂e
Zeu+ε2∂e
Xew2ε∂e
Zew
.
No e also ha ∥e
D(e
u)∥=∂e
Zeu/2 + O(ε). Finally, he F oude numbe is de ined as F =U/√gH cos θ.
Following he asymp o ic expansion pe o med in [42] o he p essu e, and in pa icula o he non-hyd os a ic
coun e pa , we w i e
ep=ρ0
F 2eb+eh−e
Z+εg
pnh
1+ε2g
pnh,(2)
whe e he p essu e is di ided in o a hyd os a ic con ibu ion, and he i s and second o de non-hyd os a ic con ibu ions
g
pnh
1,g
pnh. In ha wo k, looking a he non-dimensional e ical momen um equa ion, hese con ibu ions a e iden i ied o
he s ess enso con ibu ions and he no mal accele a ion, espec i ely. The e, i s o de e ms we e neglec ed whe eas
second o de e ms we e kep . The eason is i s ha dealing wi h he e ms coming om he s ess enso is no ably mo e
di icul , especially om he nume ical poin o iew. Fu he mo e, including he no mal accele a ion (in his sense, he
ob ained model is called weakly non-hyd os a ic) in [42] allowed o ob ain e y p omising esul s o he dep h-a e aged
(e.g. one-laye ) model. I is he e o e in e es ing o ex end his weakly non-hyd os a ic model o i s mul ilaye e sion.
To his aim, we ollow he same app oach he e. Taking in o accoun (2), he non-dimensional 2D Na ie -S okes sys em
eads ( ildes a e d opped o he sake o simplici y)















∂Xu+∂Zw= 0,
ρ0∂ u+u∂Xu+w∂Zu+∂Xp=1
ε
ρ0
F 2 an θ+ε∂XτXX +1
ε∂ZτXZ ,
ρ0ε2∂ w+u ∂Xw+w ∂Zw+ε∂Zpnh
1+ε2∂Zpnh =ε∂XτZX +ε∂ZτZZ .
(3a)
(3b)
(3c)
5

This is he sys em ha will be laye -a e aged using he mul ilaye app oach. Be o e, bounda y condi ions a e w i en
in non-dimensional o m as
∂ h+u|b+h∂X(b+h)−w|b+h= 0, pnh
|b+h= 0,
a he ee su ace, and
u|b∂Xb−w|b= 0,η
2∂Zu|b
=µ(I|b)p|b
u|b
u|b,
a he bo om.
In nex sec ion, we de elop he mul ilaye disc e iza ion and laye -a e aging o sys em (3).
3 Mul ilaye sys em wi h he µ(I)- heology and a weakly non-hyd os a ic
p essu e
In his subsec ion we apply he mul ilaye disc e iza ion in oduced in [39] o sys em (3). We i s ecall he main aspec s
o his app oach and we ocus la e on he e ms ha di e om he p e ious wo k [36], whe e all he de ails o he
hyd os a ic case a e desc ibed.
Figu e 2: Ske ch o he mul ilaye disc e iza ion.
We i s conside (lα)1≤α≤Ns ic ly posi i e coe icien s e i ying PN
α=1 lα= 1. Then, o a posi i e heigh h( , X) we
conside Nshallow laye s whose heigh s a e hα( , X) = lαh( , X) (see Figu e 2). No ice ha hen PN
α=1 hα=hholds.
Conside ing now he le els Zα+1/2=b+Pα
β=1 hβ, we deno e he laye s by
Ωα( , X) = {( , X, Z) : Zα−1/2≤Z≤Zα+1/2}, o α= 1, . . . , N.
No e ha b=Z1/2<··· < Zα+1/2<··· < ZN+1/2=b+h
is a pa i ion o [b, b +h]. As no a ion, we de ine
α=1
hαZZα+1/2
Zα−1/2
( , X, Z)dZ,
he a e aged alue in he laye Ωαo he unc ion . In o de o deal wi h possible discon inui ies o ac oss he in e ace
de ined by Z=Zα+1/2, we de ine
−
α+1/2= lim
Z→Zα+1/2
Z<Zα+1/2
|Ωα( ), +
α+1/2= lim
Z→Zα+1/2
Z>Zα+1/2
|Ωα+1( ),
6
and [[ ]]α+1/2= +
α+1/2− −
α+1/2deno es he jump o ac oss he in e ace. In case o a con inuous unc ion, hen
[[ ]]α+1/2= 0 and he e o e +
α+1/2= −
α+1/2= α+1/2.
In o de o ob ain he mul ilaye disc e iza ion o sys em (3), some assump ions on he unknowns (u, w, q) a e needed.
Then, he mass and momen um equa ions a e in eg a ed on each laye Ωα, leading o 3 laye -in eg a ed equa ions o
each laye .
In his pape , we p esen a gene aliza ion o he model deno ed by LDNH0in [40]. In his model, he wo componen s
o he eloci y uand wa e app oxima ed by cons an unc ions in each laye :
u( , X, Z) =
N
X
α=1
uα( , X)
1
Ωα(Z), w( , X, Z) =
N
X
α=1
wα( , X)
1
Ωα(Z),
and he non-hyd os a ic p essu e pnh( , X, Z), pnh
|Ωα∈P1by a con inuous linea unc ion. No e ha in his case:
pnh
α( , X) = pnh( , X, Zα) = pnh
α+1/2( , X) + pnh
α−1/2( , X)
2(4)
whe e Zα=Zα−1/2+hα/2 is he midpoin o he laye Ωα. Ac ually, o each laye we ha e
pnh(Z) = pnh
α+pnh
α+1/2−pnh
α−1/2
hα
(Z−Zα), o Z∈[Zα−1/2, Zα+1/2].
Unde hese assump ions and by app op ia ely de ining he jump condi ions ac oss he in e aces, he LDNH0model
o he Eule equa ion is ob ained in [40]. He e, he main no el y is o deal wi h he iscous e m ∂ZτXZ , which in ol es
a complex iscosi y, as de ailed in he ollowing sec ions.
3.1 Jump condi ions a he in e aces
F om an asymp o ic analysis, we eco e he ollowing non-dimensional condi ions when imposing he no mal lux jump
condi ion o he mass and he momen um equa ions (all he de ails a e p esen ed in [36]):
The jump condi ion associa ed wi h he mass conse a ion equa ion leads o he de ini ion o he mass ans e ence
e m a he in e ace Z=Zα+1/2, deno ed by Gα+1/2:= G+
α+1/2=G−
α+1/2(co esponding o −Γα+1/2in [40]), whe e
G±
α+1/2=∂ Zα+1/2+u±
α+1/2∂XZα+1/2−w±
α+1/2.
Mo o e , om ha condi ion we eco e
[[u]]α+1/2·nα+1/2= 0,
whe e nα+1/2=∂XZα+1/2,−1is he downwa d no mal ec o a he in e ace.
F om he jump condi ion o he momen um equa ion, i ollows ha
ε2τ±
XX,α+1/2∂XZα+1/2−τ±
XZ,α+1/2=ε2eτXX,α+1/2∂XZα+1/2−eτXZ,α+1/2±1
2ρ0εGα+1/2u+
α+1/2−u−
α+1/2,
τ±
ZX,α+1/2∂XZα+1/2−τ±
ZZ,α+1/2=eτZX,α+1/2∂XZα+1/2−eτZZ,α+1/2±1
2ρ0εGα+1/2w+
α+1/2−w−
α+1/2,
whe e we ecall ha τXX =η∂Xu,τXZ =1
2η∂Zu+ε2∂Xw,τZZ =η∂Zw, and eτis a consis en app oxima ion o τa
he Z=Zα+1/2(see [36]).
3.2 Laye -a e aging p ocedu e
In his subsec ion we de ail he laye -a e aging o sys em (3) in he amewo k o he mul ilaye app oach. We ackle he
mass and he momen um conse a ion equa ion sepa a ely.
7
3.2.1 Mass conse a ion equa ion
The mass conse a ion equa ion (3a) is in eg a ed along he no mal di ec ion in each laye , leading o
0 = ZZα+1/2
Zα−1/2
(∂Xu+∂Zw)dZ =∂X(hαuα)−u−
α+1/2∂XZα+1/2−w−
α+1/2+u+
α−1/2∂XZα−1/2−w+
α−1/2,
and he e o e
lα(∂ h+∂X(huα)) = Gα+1/2−Gα−1/2, o α= 1, . . . , N.
No ice ha combining p e ious equa ions, he explici o mula o Gα+1/2is ound (see [2])
Gα+1/2=
α
X
β=1
lβ∂X(h(uβ−u)) ,wi h u=
N
X
γ=1
lγuγ.(5)
In pa icula , by summing up all he equa ions we ha e he mass conse a ion equa ion
∂ h+∂X(hu)=0,(6)
whe e G1/2=GN+1/2= 0 ha e been assumed as bounda y condi ions.
3.2.2 Ho izon al momen um conse a ion equa ion
We spli in his case equa ion (3b) in con ec i e, p essu e, and iscous e ms. The usual in eg a ion p ocedu e gi es
ZZα+1/2
Zα−1/2
(∂ u+u∂Xu+w∂Zu)dZ =∂ (hαuα) + ∂Xhαu2
α−u−
α+1/2Gα+1/2+u+
α−1/2Gα−1/2.
No e ha in he p esen app oach (LDNH0model) i holds ha u−
α+1/2=u+
α−1/2=uα, since u|Ωα∈P0. P essu e e ms
a e now in eg a ed wi hin he laye in he no mal di ec ion, leading o
ZZα+1/2
Zα−1/2
∂Xρ0
F 2eb+b+h−z+εpnh
1+ε2pnhdz =hα
ρ0
F 2∂Xeb+b+h
+∂Xεhαpnh
1,α +ε2hαpnh
α−εpnh
1+ε2pnh|Zα+1/2
∂XZα+1/2+εpnh
1+ε2pnh|Zα−1/2
∂XZα−1/2,
whe e we ha e used eb(X) = (Xend −X) sin θ. Finally, iscous e ms a e in eg a ed, which yields
ZZα+1/2
Zα−1/2ε∂XτXX +1
ε∂ZτXZ dz
=ε∂X ZZα+1/2
Zα−1/2
η∂Xu dz!+1
εε2τ+
XX,α−1/2∂XZα−1/2−τ+
XZ,α−1/2−1
εε2τ−
XX,α+1/2−τ−
XZ,α+1/2
=ε∂X ZZα+1/2
Zα−1/2
η∂Xu dz!+εeτXX,α−1/2∂XZα+1/2−1
εeτXZ,α−1/2−εeτXX,α+1/2∂XZα+1/2−1
εeτXZ,α+1/2
+1
2ρ0Gα−1/2u+
α−1/2−u−
α−1/2+1
2ρ0Gα+1/2u+
α+1/2−u−
α+1/2.
De ining now
˜
Kα+1/2=εeτXX,α+1/2∂XZα+1/2−1
εeτXZ,α+1/2and UH
Zα+1/2=uα+1 −uα
hα+1/2
,
we ob ain ha
˜
Kα+1/2=Kα+1/2+O(ε),whe e Kα+1/2=−1
ε
ηα+1/2
2UH
Zα+1/2.
8
We can collec inally all he e ms, keep e ms up o i s o de and e ms o o de ε2, go o he dimensional a iables,
and ob ain he ho izon al momen um equa ion
ρ∂ (hαuα) + ∂Xhαu2
α+gcos θhα∂Xeb+b+h+∂Xhαpnh
α=Kα−1/2−Kα+1/2
+pnh
α+1/2∂XZα+1/2−pnh
α−1/2∂XZα−1/2+uα+1/2ρGα+1/2−uα−1/2ρGα−1/2,(7)
o α= 1, . . . , N and
uα+1/2=uα+uα+1
2.
3.2.3 Ve ical momen um conse a ion equa ion
Following he same p ocedu e o equa ion (3c), we ob ain he momen um conse a ion equa ion in he no mal di ec ion
o each laye . Le us jus de ail he iscous e ms. We ha e
ZZα+1/2
Zα−1/2
ε(∂XτZX +∂ZτZZ )dz
=ε∂X ZZα+1/2
Zα−1/2
η∂Zu+ε2∂Xwdz!+ετ+
ZX,α−1/2∂XZα−1/2−τ+
ZZ,α−1/2−ετ−
ZX,α+1/2∂XZα+1/2−τ−
ZZ,α+1/2
=ε∂X ZZα+1/2
Zα−1/2
η∂Zu+ε2∂Xwdz!+εeτZX,α−1/2∂XZα−1/2−eτZZ,α−1/2−εeτZX,α+1/2∂XZα+1/2−eτZZ,α+1/2
+1
2ρ0ε2Gα−1/2w+
α−1/2−w−
α−1/2+1
2ρ0ε2Gα+1/2w+
α+1/2−w−
α+1/2.
I is obse ed ha
ε∂X ZZα+1/2
Zα−1/2
η∂Zu+ε2∂Xwdz!+εeτZX,α−1/2∂XZα−1/2−eτZZ,α−1/2−εeτZX,α+1/2∂XZα+1/2−eτZZ,α+1/2
is de ined only by e ms o o de O(ε) and O(ε3).
Then, collec ing all he e ms and keeping also e ms o o de ε2, he e ical momen um equa ion in dimensional
a iables is
ρ(∂ (hαwα) + ∂X(hαuαwα)) = pnh
α−1/2−pnh
α+1/2+wα+1/2ρGα+1/2−wα−1/2ρGα−1/2,(8)
o α= 1, . . . , N, whe e wα+1/2= (w+
α+1/2+w−
α+1/2)/2. No ice ha he e i s o de e m in εha e been neglec ed whe eas
second o de e ms a e kep , simila ly o he hypo hesis made o he p essu e (2). On he one hand, he disc e iza ion o
he neglec ed i s o de e ms by non-hyd os a ic mul ilaye models has no been ackled ye in he li e a u e. I en ails
impo an di icul ies e en in he case o a hyd os a ic p essu e (see [18]). On he o he hand, including second o de
e ms ga e p omising esul s in [42].
3.3 Final mul ilaye non-hyd os a ic model
No ice ha we ha e 3N+1 unknowns (h, uα, wα, pnh
α) and 2N+1 equa ions ill now. Then, he incomp essibili y condi ion
is used in each laye . Conc e ely, we in eg a e he incomp essibili y equa ion (3a) o Z∈[Zα−1/2, Zα] and, aking in o
9
S ep 3: Non-hyd os a ic p essu e co ec ion
In he las s ep, he non-hyd os a ic e ec s a e added using he momen um equa ions (9b)-(9c) oge he wi h he
incomp essibili y cons ain s (9d). Thus, we implici ly sol e he ollowing sys em:
(∂ U=T(U,P, ∂XU, ∂XP, ∂XZ),
C(U, ∂XU, ∂XZ) = 0,(26)
whe e Tis e alua ed a ime n+1 co esponding o he las e m in (17). As usual, we will use an implici p ojec ion
co ec ion me hod ha will lead us o sol e a linea sys em.
Le us ecall he cons ain a he laye αby
Eα= 0,∀α= 1, . . . , N, whe e Eα=wα−uα∂XZα+
α−1
X
β=1
∂X(lβhuβ) + ∂X(lαhuα)
2.
As i can be seen, he imposi ion o he cons ain o α= 2, . . . , N in ol es a iables om he laye β= 1 o β=α−1.
To simpli y he se o cons ain s (9d), and o imp o e he diagonal dominance o he linea ope a o o be in e ed, we
equi alen ly conside
Nα= 0,∀α= 1, . . . , N, whe e N1=E1,N2=E2−E1, . . . , Nα=EN−EN−1.
Fo cla i y, we will assume in wha ollows ha lα= 1/N. Then, we de ine he ec o Cin (13b) by
C(U, ∂XU, ∂XZ)=2hN (N1, . . . , NN)′,(27)
ha eads
C=






2Nhw1−2Nhu1∂xZ1+h∂x(hu1)
2N(hw2−hw1)−2N(hu2∂xZ2−hu1∂xZ1) + h∂xhu2+h∂xhu1
.
.
.
2N(hwN−hwN−1)−2N(huN∂xZN−huN−1∂xZN−1) + h∂xhuN+h∂xhuN−1







,(28)
whe e we ecall ha Zα=Zα+1/2−1
2Nh. We conside now a p ojec ion me hod, aking in o accoun he luc ua ion
e
Pn+1 de ined in (15). To do ha , we will impose he incomp essibili y cons ain s on he s agge ed mesh xi+1/2
CUn+1
i+1/2,(∂XUn+1)i+1/2,(∂XZ)i+1/2=0,(29)
whe e
Un+1
i+1/2=Un+2/3
i+1/2+ ∆ TUn+1
i+1/2,e
Pn+1
i+1/2,(∂XUn+1)i+1/2,(∂Xe
Pn+1)i+1/2,(∂XZ)i+1/2,(30)
and he conse ed a iables and hei de i a i es on he s agge ed mesh a e app oxima ed om he a e aged alues
Un+2/3
i+1/2=Un+2/3
i+Un+2/3
i+1
2,Un+1
i+1/2=Un+1
i+Un+1
i+1
2,(∂XUn+1)i+1/2=Un+1
i+1 −Un+1
i
∆x.
Simila ly, he componen (∂XZα−1/2) o (∂XZ) is gi en by
(∂XZα−1/2)i+1/2=h+
i+1/2−h−
i+1/2−(hi+1 −hi)
∆x+α−1
N·hi+1 −hi
∆x,
whe e he i s addend is an app oxima ion o ∂xbusing he hyd os a ic econs uc ion (20). Finally, he de i a i es o
he luc ua ion o he non-hyd os a ic p essu es a e app oxima ed by
(∂Xe
Pn+1)i+1/2=e
Pn+1
i+3/2−e
Pn+1
i−1/2
2∆x.
16

Rema k 3 No ice ha he ac o hin he de ini ion o Cin (27)leads o a ew i ing o he incomp essibili y cons ain s
in e ms o solely conse ed a iables (see (28)). Tha gi es a mo e e icien and s able nume ical implemen a ion since
no di ision by hmus be pe o med when compu ing C.
Now, pu ing Eqs. (30) in o (29) yields a linea sys em
A·P=1
∆ R,(31)
whe e Pis a ec o ha con ains a ea angemen o he non-hyd os a ic p essu e luc ua ions unknowns
P=


P1
.
.
.
PN


,Pα=




epnh,n+1
α,1/2
.
.
.
epnh,n+1
α,Nx+1/2





, α = 1, . . . , N,
Ais a block idiagonal ma ix o dimension N·(Nx+ 1) ×N·(Nx+ 1) gi en by
A=






M1,1M1,2
M2,1M2,2M2,3
M3,2
......
......MN−1,N
MN,N−1MN,N







,(32)
being Mi,j idiagonal ma ices o dimension (Nx+1)×(Nx+1) (see de ails in Appendix A). Finally, Ris he igh -hand
side o he linea sys em ha is gi en by
R=


R1
.
.
.
RN


,Rα=






CαUn+2/3
1/2,(∂XUn+2/3)1/2,(∂XZ)1/2
.
.
.
CαUn+2/3
Nx+1/2,(∂XUn+2/3)Nx+1/2,(∂XZ)Nx+1/2







, α = 1, . . . , N, (33)
whe e he sub-index αin Cαdeno es he α− h componen o he ec o (28).
To nume ically app oxima e he solu ion o he linea sys em (31), we conside he ollowing block Gauss-Seidel sol e :
Algo i hm
1: P(0) ←I.C.
2: s←1
3: E ←0
4: do
5: o α= 1,2, . . . , N do
6: Sol e he idiagonal sys em: Mα,αP(s)
α=1
∆ Rα−Mα,α−1P(s−1)
α−1−Mα,α+1P(s−1)
α+1 ,
7: E ←E + ∆x||P(s−1)
α−P(s)
α||p
8: P(s−1)
α←P(s)
α
9: end o
10: while E ≥ϵp
ol.
In line 1 o he algo i hm, I.C. deno es he ini ial condi ion o he i e a i e sol e . In p ac ice, we ake as ini ial condi ion
he luc ua ion a he p e ious ime s ep e
Pn=Pn−Pn−1,al hough o he choices a e possible. As s a ed in line 6 o
17
he algo i hm, a idiagonal linea sys em mus be sol ed o each laye and a each loop o he i e a i e me hod, which
is ca ied ou wi h he e icien Thomas’ algo i hm. The E o a iable E is used o compu e he lpno ms o p= 1,o
p= 2,used as he s opping c i e ia and gi en by
||P(s−1)
α−P(s)
α||p=1
N
Nx
X
i=0
hi+1/2epnh,(s−1)
α,i+1/2−epnh,(s)
α,i+1/2
p, o p= 1,2.
Al e na i e, he l∞no m can be conside ed, and he e o e E in line 7 is upda ed as
E ←max E , max
i=0,...,Nxnepnh,(s−1)
α,i+1/2−epnh,(s)
α,i+1/2o‘
whe eas he s opping condi ion in line 10 is gi en by E ≥ϵp
ol.The in luence o he choice o ha di e en no ms will be
s udied in Subsec ion 5.1.1. Finally, ϵp
ol will be se in he nume ical expe imen s.
Rema k 4 The p oposed linea sol e has shown he abili y o con e ge o a ew i e a ions (less han 62 o 5laye s,
as shown in Table 2) conside ing he s a egy desc ibed in Subsec ion 4.1, ha can be imp o ed.
One could use o he linea sol e s, such as LU ac o iza ions o K ylo g adien -based me hods. The i s one is
disca ded since he ma ix o be in e ed (32)depends on he ime s ep, inc easing he compu a ional e o . K ylo -based
me hods con e ge apidly, especially i a con enien p econdi ione is adop ed and a ma ix- ee e sion o he algo i hm o
main ain low s o age o he inal implemen a ion is conside ed (see, o ins ance [64]). In gene al, ma ix- ee linea sol e s
make he algo i hm mo e con enien o u u e GPU implemen a ions o 3D domains whe e memo y is an impo an
issue.
The me hodology employed he e ollows simila ideas p esen ed in [31,59], whe e he esul ing linea sys ems a e
sol ed using i e a i e Jacobi me hods combined wi h a scheduled elaxa ion pa allelized o 3D domains ha show excellen
compu a ional pe o mance. Some o he au ho s ha e employed simila me hods. Fo ins ance, in ([31,33]), di e en
ad-hoc implemen a ions o he Gauss-Seidel i e a i e p ocedu es we e employed o 2D domains.
Once he non-hyd os a ic p essu e luc ua ions e
Pn+1
i+1/2ha e been compu ed, he discha ges Un+1
ican be upda ed
om
Un+1
i=Un+2/3
i+ ∆ TUn+1
i,e
Pn+1
i,(∂XUn+1)i,(∂Xe
Pn+1)i,(∂XZ)i,
whe e simila ly, a second o de poin alue app oxima ion in he cen e o he cell will be used using (22). Finally, he
non-hyd os a ic unknowns a e upda ed o Pn+1
i+1/2by (16).
Rema k 5 Since non-hyd os a ic e ms appea only in he momen um equa ions, he inal nume ical scheme is
well-balanced o he wa e a es solu ions hanks o he hyd os a ic econs uc ion pe o med in s ep 1 and he de ini ion
(25), as we p o e nume ically in Subsec ion 5.3. I is posi i e p ese ing o he wa e heigh . I is a consequence o he
nonnega i i y o he HLL scheme o a la bo om ([9]) and he nume ical ea men o we /d y a eas [21] (see Rema k
2). Mo eo e , he inal nume ical scheme is linea ly L∞−s able unde he usual CFL condi ion gi en in (14)(see [22]).
Rema k 6 Conce ning bounda y condi ions, he e we use s anda d me hods o wall/ ee-ou low bounda y condi ions as
in [31]. Conc e ely, we use a ghos -cell echnique ha duplica es he heigh and he p essu es a he bounda ies, and ei he
duplica es o akes he opposi e alue o he discha ges, depending on whe he ee-ou low o wall bounda y condi ions a e
conside ed. Tha is, a disc e iza ion o homogeneous Neumann bounda y condi ions a e used o he heigh and p essu es,
and homogeneous Neumann ( ee-ou low case) o Di ichle (wall case) bounda y condi ions a e used o he discha ges.
In addi ion, he ma ix A(32) o he p essu e esolu ion is modi ied acco ding o (37).
4.1 S a egy o speed up non-hyd os a ic mul ilaye simula ions
As commen ed a he beginning o Sec ion 4, he compu a ional cos associa ed wi h non-hyd os a ic mul ilaye sys ems
becomes high, especially when inc easing he numbe o e ical laye s. Mos o his e o is due o he hi d s ep in he
scheme, ha is, he ellip ic p oblem sol ing he non-hyd os a ic con ibu ions. Fo his eason, any s a egy o speed up
his s ep is app ecia ed. We p opose a s aigh o wa d s a egy ha allows us o signi ican ly educe he compu a ional
cos when using many laye s. This s a egy will be analyzed om he nume ical poin o iew in Subsec ion 5.1.
18
The goal is o educe he numbe o i e a ions in he non-hyd os a ic i e a i e sol e , o a leas when conside ing a
high numbe o e ical laye s. We achie e i by imp o ing he ini ial guess P(0), aking he solu ion p o ided by he
linea sys em a he p e ious ime s ep.
The main idea is o sol e he ellip ic p oblem ela ed o dispe si e e ec s in a coa se e ical mesh, wi h a lowe
numbe o laye s Nini guess < N, so ha i s solu ion is a be e ini ial guess o he o iginal sys em wi h Nlaye s. In
his way, we would sol e he dispe si e linea sys em wice: i s ly wi h Nini guess laye s, whose esolu ion is much as e
han he o iginal p oblem; secondly he ull sys em wi h an imp o ed ini ial guess, so ha less i e a ions a e needed o
con e ge o he solu ion Pn+1.
In o de o be able o build he i s linea sys em in he coa se e ical disc e iza ion, we assume a con o mal mesh
as ollows. We de ine Ng =N/Nini guess ∈N ep esen ing he numbe o laye s ha me ge in o a single laye . Then, he
domain is di ided along he no mal di ec ion in o laye s Ωβ, ed, whe e
Ωβ, ed =
Ng
[
j=1
ΩNg (β−1)+j o β= 1,...Nini guess.
Now, we need o de ine auxilia y ec o a iables wi h size Nini guess o his domain subdi ision. Those a e gi en om
he o iginal ec o a iables by simply compu ing app op ia e a e aged alues wi h he da a o he me ged laye s. Fo
ins ance, o he ho izon al eloci ies, we de ine
uβ, ed =1
Ng
Ng
X
j=1
uNg (β−1)+j o β= 1,...Nini guess,
and analogously o all he a iables in ol ed in esol ing he ellip ic p oblem (30). Conc e ely, i is done o de ine U ed,
P ed and e
P ed, which a e he inpu a iables o he i s s age whe e he linea sys em o he educed ellip ic ope a o
is sol ed. Once we ha e he solu ion o he non-hyd os a ic a iable e
P ed, i is used o gene a e he ini ial guess o he
second s age, whe e he ull linea sys em wi h Nlaye s is sol ed. Thus, we de ine, o α= 1, . . . , N,
epnh,(0)
α+1/2=epnh
β+1/2, ed
wi h 1 ≤β≤Nini guess such ha Zβ−1/2≤Zα−1/2< Zβ+1/2, and epnh,(0)
N+1/2= 0. These a iables de ine he i e a i e
algo i hm’s ini ial condi ion P(0). I is expec ed ha he numbe o i e a ions needed in his case o he ull sys em
is lowe han he o iginal one, whe e he solu ion a he p e ious ime s ep is aken. The e o e, he ellip ic p oblem’s
compu a ional ime is also educed, despi e sol ing wo linea sys ems. In Subsec ion 5.1.2 we will see ha his s a egy is
e ec i e o a g anula collapse p oblem and di e en e ical disc e iza ions. Le us ema k ha , al hough his s a egy
is e ec i e o he nume ical expe imen s pe o med he e, o he s a egies could be used o speed up he i e a i e sol e .
5 Nume ical es s
We pe o m he e some nume ical es s wi h he p oposed mul ilaye model and compa e he esul s o hose ob ained
wi h p eceden models, namely he mul ilaye hyd os a ic [36] and he one-laye weakly non-hyd os a ic [42] models. The
goal is o show ha he mul ilaye model wi h weakly non-hyd os a ic p essu e combines some o he s eng hs o hese
p e ious models. In pa icula , we will quan i y he in luence o he mul ilaye and o he non-hyd os a ic app oaches in
se e al p ac ical cases. I will be shown in Subsec ion 5.5. In pa icula , we eco e impo an esul s o bo h models,
al hough he p ize is a highe compu a ional cos . Be o e, conce ning his compu a ional e o , i will be e alua ed in
Subsec ion 5.1, oge he wi h he e iciency o he p oposed nume ical s a egy o speed up he non-hyd os a ic s ep o
he nume ical scheme p esen ed in Sec ion 4. Nex , in Subsec ion 5.4, we will quan i y he in luence o he coo dina e
sys em (Ca esian o local) in he mul ilaye model in o de con i m he esul s ob ained wi h he one-laye model (see
[42]). No e ha his canno be done by simula ing g anula collapse expe imen s since he associa ed ini ial condi ions
canno be de ined in a Ca esian e e ence sys em. We hen mo e o he compa ison wi h labo a o y g anula collapses
expe imen s.
In all he es s, we conside he ma e ial p ope ies gi en in Table 1. In pa icula , he µ(I) coe icien s a e inc eased
by 0.1 o accoun o sidewall ic ion, ollowing he app oach o [52].
Fo he nume ical pa ame e s, unless o he wise is speci ied, we conside ϵp
ol = 10−6as ole ance o he i e a i e sol e
inding he non-hyd os a ic p essu e (see line 10 o he algo i hm in Sec ion 4), and δ= 10−5as egula iza ion pa ame e
in he iscosi y coe icien (equa ion (10)). No e ha his choice o δis in ag eemen wi h o he wo ks o iscoplas ic
luids, as [41,25,58].
19
µsµ2ds(mm) I0φs
an(25.5◦)≃0.48 an(36◦)≃0.73 0.7 0.279 0.62
Table 1: Values o he heological pa ame e s used in he nume ical es s.
5.1 Analysis o he compu a ional e o s o non-hyd os a ic simula ions
We in es iga e he e he compu a ional e o s associa ed o he non-hyd os a ic mul ilaye model. The compu a ion o
he non-hyd os a ic p essu e desc ibed in Sec ion 4is based on an i e a i e sol e o sol e a linea sys em whose size
is N×(Nx+ 1), whe e Nis he numbe o laye s in he no mal di ec ion and Nx he numbe o cells in he ho izon al
disc e iza ion. Then, he compu a ional ime could g ow d ama ically when inc easing he numbe o no mal laye s. This
ac mo i a es he ollowing a emp o educe he compu a ional e o .
We i s s udy he in luence o changing he no m used in he s opping c i e ion o he i e a i e sol e o he
non-hyd os a ic p essu e, and secondly we e alua e he p oposed s a egy in Subsec ion 4.1 o speed up non-hyd os a ic
simula ions.
We conside he e a compu a ional domain [−0.2,3] mand 640 cells o he ho izon al disc e iza ion, and he dam
b eak p oblem gi en by he ini ial condi ions
eb(X) = (3 −X) sin θ, b(X)=0, h0(X) = 0.14 + hii x≤0,
hio he wise,(34)
wi h hi= 1.82 mm and θ= 22◦. We choose a high slope case since he mul ilaye amewo k is mo e app op ia e o la ge
alues o he slopes, as i will be concluded in 5.5. The ma e ial pa ame e s a e in Table 1. A wall bounda y condi ion
is employed ups eam while a ee-ou low condi ion is used downs eam. Finally, we use CFL = 0.7 and = 3 s.
5.1.1 In luence o he choice o he s opping c i e ion in he i e a i e sol e
As explained in Sec ion 4, he compu a ion o he non-hyd os a ic p essu e is based on an i e a i e sol e o ind he
solu ion o a linea sys em o size N×(Nx+ 1). The e o e, he compu a ional ime is sensi i e o he no m used in he
s opping c i e ion o his sol e .
Laye s l2- NºI e l2- comp l∞- NºI e l∞- comp
20 431 779.7 693 1096.7
10 145 84.8 211 197.7
8 101 58.8 142 114.3
5 46 18.1 62 36.4
2 10 4.2 14 5.2
Table 2: Compu a ional imes comp (s) and maximum numbe o i e a ions made in a single ime s ep by he i e a i e
sol e o di e en numbe o no mal laye s Nand no ms l2,l∞.
We ha e measu ed he compu a ional ime and he numbe o i e a ions made by he sol e o di e en numbe o
e ical laye s. Table 2shows he compu a ional imes and he maximum numbe o i e a ions needed in a ime s ep
o N= 20,10,8,5,2 no mal laye s, and Figu e 3shows he i e a ions needed a each ime s ep o N= 5,10,20 laye s
(Figu e 3a) and he compu a ional imes o bo h no ms (Figu e 3b). When using he l2no m he compu a ional e o
is lowe han wi h he l∞no m o all cases, as expec ed. In pa icula , o N= 20, he simula ion needs 46% less ime.
I is due o he ac ha he numbe o i e a ions needed wi h he l∞no m is always g ea e han wi h he l2no m (see
Figu e 3a). In pa icula , we see he di e ence in he numbe o i e a ions needed when using he di e en no ms in he
s opping c i e ion, a e a small ime. No e ha when using he l1no m, he esul s a e simila o he l2no m.
In addi ion, we ha e checked ha he e a e no signi ican di e ences in he esul s when using he di e en no ms in
he i e a i e sol e . Then, he l2-no m is used om now on in he simula ions, also in subsec ions 5.4 and 5.5.
Mos o he compu a ional cos is spen in he S ep 3, ha is, he ellip ic p oblem ela ed o he non-hyd os a ic
p essu e. In Table 3we show he pe cen age o he compu a ional ime ha he algo i hm spends in his s ep and also
o S ep 2 ( iscous e ms) o di e en numbe o laye s. We see ha i is 28% o he one-laye model, bu i inc eases
quickly wi h he numbe o laye s. Fo 5 laye s i is 81%, and i eaches 98% o 20 laye s. In he non-hyd os a ic case,
20
0123
0
50
100
150
200
250
0123
0
200
400
600
800
0 2 5 8 10 15 20
1
2
3
5
13
18
Figu e 3: (a): Numbe o i e a ions needed o N= 5,10,20 laye s wi h he l2and l∞no ms. (b): Compu a ional ime
in minu es o di e en numbe o laye s wi h no ms l2, l∞.
we see ha S ep 2 does no spend a signi ican amoun o esou ces. We also see ha S ep 2 spends app oxima ely 10%
o he compu a ional ime o N= 1,2, ..., 20 in he hyd os a ic case.
Laye s 1 2 5 10 20
comp S ep 3/ comp (NH) 28% 58% 81% 92% 98%
comp S ep 2/ comp (NH) 5% 4% 2% 0.6% 0.2%
comp S ep 2/ comp (H) 9% 11% 11% 10% 9%
Table 3: Pe cen age o compu a ional ime spen by s eps 2 ( iscous e ms esolu ion) and 3 (ellip ic p oblem associa ed
o he p essu e) o he scheme in Sec ion 4in he non-hyd os a ic case (NH). We also show his pe cen age o S ep 2 in
he hyd os a ic case (H).
5.1.2 Speed-up s a egy o non-hyd os a ic simula ions
We in es iga e in his subsec ion he e iciency o he algo i hm o speed up he nume ical scheme o non-hyd os a ic
simula ions ha has been p oposed in Subsec ion 4.1.
(a) N= 8 laye s (b) N= 10 laye s
Nini guess - 4 2
comp (s) 47.3 44.6 77.46
Speed-up - 1.06 0.61
Nini guess - 5 2
comp (s) 84.8 56.5 128.2
Speed-up - 1.5 0.66
(c) N= 20 laye s
Nini guess - 10 5 4 2
comp (s) 594.2 183.2 166.8 215.2 564.6
Speed-up - 3.24 3.56 2.76 1.05
Table 4: Compu a ional imes and speed-ups o non-hyd os a ic simula ions wi h (a) 8 (b) 10, (c) 20 e ical laye s,
educing o Nini guess laye s o he i s s ep o he algo i hm o sol e he non-hyd os a ic p essu e.
21

0 5 10 15 20
0
2
4
6
8
10
0 0.1 0.2 0.3 0.4 0.5
0
100
200
300
400
0.6 1 1.5 2 2.5
0
50
100
Figu e 4: (a): Compu a ional ime in minu es o di e en numbe o laye s used o he p ecompu ing o he ini ial guess
o he p essu e sol e o N= 8,10,20 laye s. (b): Numbe o i e a ions needed by he i e a i e sol e in he case
N= 20 o he simula ion wi hou using a p ecompu ed guess (FS 20 laye s, dashed black line) and when an ini ial guess
is p ecompu ed wi h 10 laye s ( i s s ep wi h educed sys em wi h 10 laye s - RS-20-10 - solid cyan line and second s ep
wi h he ull sys em - FS-20-10 - dashed ed line). Do ed g ay line is he o al o i e a ions (s ep 1 and s ep 2) in he
case o p ecompu ing he ini ial guess.
In Table 4we show he compu a ional imes when educing o Nini guess laye s in he i s s ep o he algo i hm o
compu e a good ini ial guess o he p essu e i e a i e sol e o simula ions wi h 8, 10 and 20 laye s. We see ha he
p oposed algo i hm is e ec i e, especially when using a la ge numbe o no mal laye s. Ac ually, he e is no signi ican
gain in compu a ional e o when using 8 laye s. Howe e , i is use ul o 10 and 20 laye s. When he e ical esolu ion is
inc eased o 20 laye s, he p oposed s a egy allows us o signi ican ly educe he compu a ional cos wi h an app op ia e
choice o Nini guess. In his case, when educing o Nini disp = 10 laye s o he i s s ep o he p essu e compu a ion, he
speed-up is 3.24, which means ha wi h he p oposed s a egy he scheme is as e by 69% han he scheme wi hou using
he p ecompu ed ini ial guess. Fo Nini guess = 5 he speed-up is 3.56 (72% as e ) and o Nini guess = 5 he speed-up
is 2.76 (64% as e ). When using Nini guess = 2 he e is no signi ican gain in compu a ional ime (i is only 5% as e ).
Howe e , his s a egy is no e ec i e o all he choices. Fo ins ance, when using Nini disp = 2 wi h 10 o 8 laye s, he
compu a ional ime inc eases. I is clea ha in he case o 20 laye s he op imal choice is o use Nini disp = 5.
In Figu e 4we s udy he case o 20 laye s and Nini disp = 10 in mo e de ail. Figu e 4a shows he compu a ional
imes needed o di e en choices o Nini disp, while in Figu e 4b we see an analysis when Nini disp = 10 in e ms o he
numbe o i e a ions needed. The mo e expensi e pa o he scheme is he esolu ion o he non-hyd os a ic p essu e as
i is shown in Table 3, whe e he i e a i e sol e is used. We see ha he ac o sol ing a i s sys em wi h a educed
numbe o laye s (10 laye s in his case) gi es us a good ini ial guess o he esolu ion o he ull sys em wi h 20 laye s.
I makes ha he scheme wi h he p oposed s a egy needs much less i e a ions o con e ge han he scheme wi hou he
p ecompu ing o he ini ial guess. Ac ually, by sol ing di ec ly he ull sys em, he maximum numbe o i e a ions is 431,
whe eas wi h he p oposed s a egy, he maximum numbe o i e a ions needed o he ull sys em is 71.
In summa y, he p oposed s a egy is e ec i e, especially when using a la ge numbe o laye s, and i allows us o
signi ican ly educe he compu a ional ime. The key is o educe he numbe o i e a ions needed o sol e he ull linea
sys em, e en making a p e ious esolu ion wi h a smalle linea sys em a any ime s ep.
5.2 Con e gence es
We pe o m he e a con e gence es showing he i s o de accu acy o he p oposed nume ical me hod. We conside
he same es con igu a ion as in [42]. Conc e ely, we conside a g anula mass a es wi h a smoo h ee su ace, which
lows o e a slope, and measu e he o de a an in e media e ime = 0.15 s. The compu a ional domain is [−0.5,1.5] m,
CFL = 0.7 and ee-ou low bounda y condi ions a e assumed. The bo om and ini ial condi ions a e gi en by
eb(X) = (1.5−X) sin 20◦, b(X)=0, h0(X) = 0.02 + 0.4e−5(x−0.25)2.
22
The same ma e ial p ope ies as in p e ious es a e used (see Table 1).
We measu e he e he o de o accu acy o he heigh and he eloci y ield in bo h di ec ions by inc easing Nand
Nxa he same ime. Tha is, we e ine he mesh in bo h di ec ions (x, z). He e, L1absolu e e o s a e compu ed, which
a e gi en by
∥h−h e ∥1= ∆x
M
X
i=1 hi−eh e
i,∥ − e ∥1= ∆x
Nx
X
i=1
hi
N
N
X
α=1  α,i −e e
α,i ,(35)
whe e h e ,e e mus be unde s ood as he p ojec ions o he e e ence solu ions on o he spaces o he local solu ions.
In p e ious equa ion, (z) deno es he ho izon al (u(z)) o e ical eloci y (w(z)). He e he solu ion o he mul ilaye
model wi h Nx= 6400 and N= 128 will be he e e ence solu ion.
(Nx, N) E o hO de hE o u(z) O de u(z) E o w(z) O de w(z)
(100,2) 1.59×10−2- 4.14×10−2- 1.59×10−2-
(200,4) 1.22×10−20.39 2.62×10−20.66 1.22×10−20.68
(400,8) 6.40×10−30.93 1.16×10−21.18 6.40×10−31.29
(800,16) 3.06×10−31.06 5.58×10−31.05 3.06×10−31.22
(1600,32) 1.35×10−31.18 2.65×10−31.07 1.35×10−31.18
Table 5: L1e o s and ela ed o de s o he heigh (h) and he eloci y ields (u(z), w(z)) o he solu ion a ime = 0.15 s.
Table 5shows he L1e o s and nume ical con e gence a es o he heigh (h), ho izon al eloci y (u(z)) and e ical
eloci y (w(z)) a ime = 0.15. We measu e he e o s a a sho ime by wo main easons: i s , solu ions wi h
shocks and o e shoo ing may appea o longe imes; second, he compu a ional cos associa ed o he e e ence solu ion
conside ed he e (Nx, N) = (6400,128) is huge. We see ha he p oposed scheme is i s o de accu a e o all he a iables,
as expec ed.
5.3 Well-balanced es
We pe o m now a es o check he well-balancing o he scheme o non la s eady s a es de ined by (12). Le us
conside he domain [−3,1] mwi h 300 cells and ee-ou low bounda y condi ions. The ini ial condi ions (see igu e 5)
a e eb(X) = 0,
b(X) = 2
X+ 4 −1
2+e−10(X+3/2)2+ 0.15 RAND(X), h0(X)=0.8−µs(X−2) −b(X),
being RAND(X)∈[0,1] a andom alue gene a ed a un ime, and he low is assumed a es (hu)α= (hw)α= 0 o
α= 1, ..., N. We ake he heological pa ame e in Table 1excep o µs= an(20◦). CFL = 0.7 is used he e and he
inal ime simula ion is in = 1 s.
-3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1
0
1
2
3
Figu e 5: Bo om opog aphy and ini ial ee su ace he well-balanced es .
The L1e o s o all he a iables a e measu ed, as in (35), whe e in his case he e e ence solu ion e e is he ini ial
s a e. Table 6shows he maximum along he simula ed ime o hese e o s o di e en alues o he egula iza ion
23
pa ame e δ. We show hese e o s o he scheme desc ibed in Sec ion 4wi h he de ini ion o bhn+2/3
igi en by (25)
(Table 6a) and by conside ing bhn+2/3
i=hn+2/3
i(Table 6b), in o de o show i s in luence. We see he e ha he de ini ion
(25) is essen ial o ge he well-balanced p ope y o he scheme. As commen ed be o e, i is no possible o ob ain ze o
eloci y o null a ia ion o he a iables because o he egula iza ion pa ame e . In ac , he e is a ela ion be ween δ
and he o de o magni ude o he a ia ions (and he esidual eloci ies). Ac ually, hey a e o he same o de . A simila
beha iou was obse ed in [38] o a hyd os a ic mul ilaye model o iscopla ic (He schel-Bulkley) luids.
(a) WB: bhn+2/3
i(25) (b) No WB: bhn+2/3
i=hn+2/3
i
E o hE o u(z) E o w(z)
δ= 10−55.34×10−55.38×10−54.21×10−5
δ= 10−10 5.34×10−10 5.39×10−10 4.21×10−10
δ= 10−15 2.99×10−16 5.39×10−15 4.21×10−15
E o hE o u(z) E o w(z)
δ= 10−55.22×10−41.20×10−37.42×10−4
δ= 10−10 1.76×10−41.93×10−43.48×10−4
δ= 10−15 4.62×10−56.06×10−51.87×10−4
Table 6: Maximum along he ime o L1e o s o he heigh (h) and he eloci y ields (u(z), w(z)) wi h espec o he ini ial s a e by
conside ing he de ini ion o b
hn+2/3
igi en by (a) eq. (25)(well-balanced scheme); (b) hn+2/3
i(no well-balanced scheme).
Le us ema k ha , since in S ep 3 o he scheme we sol e he non-hyd os a ic p essu e, which in his es is o he
o de o δ, we mus use ϵp
ol < δ. Then, we use he e ϵp
ol =δ/100. No ice ha we eco e he expec ed well-balanced
p ope y o he scheme e en in he limi case ∂Xeb+b+h=µs.
In addi ion, i we change he slope o he ini ial condi ion by an(20.01◦) (i.e., 0.01◦g ea e han he angle o epose),
hen we eco e ha he e o be ween he ini ial and inal heigh is o o de 10−5 o any alue o δ= 10−5,10−12,10−15.
Tha is, he ma e ial is no a es in his case.
We mus also ema k ha we a e able o use e y small alues o δin his es because i is a s eady a es solu ion.
In gene al ansien p oblem we use a la ge alue (δ= 10−5) o s abili y and con e gence easons (see [58]).
In he ollowing, we mo e o es s e alua ing he abili y o he model o eco e some impo an physical beha iou s.
We i s s udy he in luence o he coo dina e sys em and hen compa e he ou esul s wi h dam b eak labo a o y
expe imen s.
5.4 In luence o he coo dina e sys em
Le us ecall he di e ences be ween he local and Ca esian models. In hin-laye landslide models, he downslope eloci y
should be compu ed in he di ec ion pa allel o he opog aphy and cu a u e e ec s should be added (see [14,60,70]).
A easonable app oach, wi h he ad an age o being much simple concep ually and om he p ac ical poin o iew, is o
ake a e e ence plane wi h a slope θco esponding o he mean slope o he opog aphy, and hen conside a ia ions o he
bo om wi h espec o his e e ence plane (see e.g. Figu e 1). I is ac ually he local coo dina e sys em conside ed in his
wo k. In his app oach, he downslope eloci y is angen o he e e ence plane eband he equa ions a e dep h-a e aged
in he di ec ion no mal o he plane. In his model, he low is conside ed hin in he no mal di ec ion compa ed o i s
ex ension pa allel o he e e ence plane. On he con a y, in Ca esian models, he eloci y is compu ed in he ho izon al
x-di ec ion, he low is assumed hin in he e ical di ec ion compa ed o he ho izon al one, and dep h-a e aging is
pe o med in he e ical di ec ion. When he slope is no small, i is well known ha hese models a e less accu a e han
local models. Ac ually, Ca esian models o e es ima e he eloci y. Ano he di e ence be ween hese e e ence sys ems in
he pa icula case o non-hyd os a ic models, is he in luence o he bo om in he e ms pnh
α±1/2Zα±1/2in he ho izon al
momen um equa ions (9b). Conc e ely, Zα±1/2is de ined in e ms o he local bo om bo om band no ebin he local
coo dina es case, con a y o he Ca esian case.
5.4.1 Nume ical se -up
To compa e he esul s in Ca esian and local coo dina e sys ems, we conside a es wi h an ini ial g anula mass ha ing
a smoo h shape so ha i s hickness can be de ined in bo h sys ems. We ake he ini ial condi ion as in [42]. I consis s
o a g anula mass lowing in a complex opog aphy wi h a bump loca ed in he middle o he domain.
We conside a Ca esian domain x∈[−3,3] m and 640 nodes. The ma e ial is ini ially a es , and we se a inal
ime = 10 s and CFL = 0.5. In his case we use ee-ou low condi ions a bo h bounda ies. Figu e 6shows he ini ial
condi ions o he bo om and he mass hickness in bo h coo dina es sys em. In he Ca esian case, hey a e
bCa (x)=1− anh(x) + κe−10(x−1)2+ 0.5e−10(x−3)2,
24
and
hCa (x) = max(η(x)−bCa (x),0) + hi,wi h η(x) = y0−1 + e−0.5(x−x0)8,i x≤0,
0,o he wise,
wi h κ= 0.3, hi= 10−4m, x0=−1.53 m and y0= 1.91 m. In ha (Ca esian) case, we ix θ= 0◦and ebCa (x) = 0.
No e also ha a small hickness hihas been added in o de o a oid dealing wi h we /d y a eas and emo e he in luence
o he ea men o such cases. In o de o de ine hese unc ions in he local coo dina e sys em, he e e ence plane
ebloc =−0.7− an θ(x−3), wi h θ= 25◦, is conside ed. Now he local bo om bloc(X) is he dis ance om ebloc(x) o he
bo om measu ed in he di ec ion no mal o he plane ebloc(x). Analogously, hloc(X) is measu ed as he dis ance om
bloc(X) o he ee su ace.
Figu e 6: Ske ch o he ini ial condi ions in bo h Ca esian and local coo dina es. κ= 0.3. The a iables de ined in he
Ca esian sys em a e in blue and hose in he Local sys em a e in ed.
5.4.2 Local/Ca esian and hyd os a ic/non-hyd os a ic model compa ison
The goal he e is o s udy wo main aspec s. We will in es iga e (i) he in luence o he Coo dina e sys em, which
should be smalle o non-hyd os a ic (NH) models han o hyd os a ic (H) ones and (ii) he simila i y be ween he
non-hyd os a ic Ca esian and he hyd os a ic local models. Conce ning (ii), [29] sugges ed ha a Ca esian model wi h
a co ec ion associa ed o he no mal accele a ion p oduces simila esul s han a hyd os a ic local model. In [42] we
ound ha i is only ue o he inal deposi bu no a in e media e imes.
Figu e 7shows he inal deposi s o he local and Ca esian models wi h a hyd os a ic and weakly non-hyd os a ic
p essu e, and 10 e ical laye s. The slope angle o he bo om bCa is also ep esen ed and a ies om −15◦ o 50◦(we
ecall ha posi i e angles co espond o nega i e slopes). We see ha , quali a i ely, he hyd os a ic local model gi es a
deposi simila o he non-hyd os a ic Ca esian model, excep nea he ini ial loca ion o he mass.
In Figu e 8a we show he ime e olu ion o he ela i e e o s be ween he mass hicknesses simula ed wi h hese wo
models. We see ha , o he mul ilaye cases (5, 10 and 20 laye s), he e o is a ound 20% a inal imes whe eas highe
e o s a e ob ained a in e media e imes. These esul s a e in ag eemen wi h [42] and con adic s wha was assumed in
[29]. In addi ion, we see ha he di e ences a e smalle in he mul ilaye case han in he one-laye case.
In Figu e 8b we show he in luence o he coo dina e sys em (Ca esian e sus Local) on he simula ions o bo h
hyd os a ic and non-hyd os a ic models. This compa ison is pe o med o he one-laye and he mul ilaye (5, 10
and 20 laye s) cases. We see ha , as expec ed, non-hyd os a ic models a e less dependen on he coo dina e sys em.
In e es ingly, he in luence o he coo dina e sys em is smalle in he mul ilaye case han in he one-laye case. Mo eo e ,
in he hyd os a ic case, he di e ences s ongly dec ease when inc easing wi h he numbe o laye s, whe eas o he
non-hyd os a ic case he esul s o 5, 10 and 20 laye s a e almos iden ical.
In o de o show how he di e ences in he models could a ec p ac ical applica ions, we conside a bo om wi h a
highe bump in he middle o he domain (a x= 1). Namely, we use κ= 0.45 he e. Figu e 9shows he deposi s o
all models. In e es ingly, he simula ed mass does no o e come he bump o he NH-Local model, whe eas he o he
models p edic ha a signi ican amoun o ma e ial is lowing u he down o he alley. Such si ua ion is usual in eal
applica ions whe e haza d assessmen equi es o es ima e i a gi en landslide will o e come a opog aphy elie such as
his bump and e en ually each a own loca ed in a alley [68,71,69]. The NH-Local model would p edic ha such
25
W(cm) model µsµ2µw
10 NH-µw an(20.9◦)≃0.38 an(32.76◦)≃0.64 an(10.5◦)
20 NH an(23.4◦)≃0.43 an(34.74◦)≃0.69 0
20 NH-µw an(20.9◦)≃0.38 an(32.76◦)≃0.64 an(10.5◦)
Table 8: Modi ied heological pa ame e s conside ed when adding sidewalls ic ion in he simula ions. Da a o NH
model wi h W= 10 m a e in Table 1.
o igh -hand side o equa ion (24). No ice ha , in p ac ice, jus he diagonal elemen o he linea sys em de ined by (24)
is modi ied. The models wi hou he co ec ion in he coe icien s (bu adding he e m (36)) a e deno ed by NH −µw,
and he heological pa ame e s used a e summa ized in Table 8 o he wid hs W= 10,20 cm.
Figu e 18 shows he mass hickness a an in e media e ime and when he ma e ial is a es (deposi )compu ed wi h
he non-hyd os a ic models wi h and wi hou he co ec ion o he ic ion coe icien s, o slopes θ= 0◦,16◦,22◦and
W= 10 cm. We see ha o θ= 0◦ he esul s a e simila in bo h cases. Howe e , o la ge alues o he slope, in
pa icula o θ= 22◦, he inal deposi s a e d ama ically o e es ima ed wi hou he ex a 0.1 in he ic ion coe icien s.
Ac ually, o his slope he low is no a es ye . No e ha in his case he slope is highe han he ic ion angle
a c an(µs) = 20.9◦. As a consequence, o la ge alues o he slope, i is no possible o emo e he co ec ion in he
ic ion coe icien s. In e es ingly, o he in e media e ime he esul s a e simila in all cases. The models NH-µw
sligh ly imp o e he esul s in he ini ial pa o he ho izon al domain. Figu e 19 shows he on eloci y in hese cases
in addi ion o he esul s o a channel wid h W= 20 cm. We see ha o high slopes, ic ion including sidewall e ec s
is no enough o s op he ma e ial, and o θ= 0◦ he e is no signi ican di e ence be ween he esul s wi h o wi hou
he sidewall ic ion e m.
Figu e 15: Flow/no- low in e ace a di e en imes a θ= 9.78◦(le -hand side) and θ= 22◦( igh -hand side) and
hi= 0 mm o he mul ilaye models wi h non-hyd os a ic (dashed g een lines) and hyd os a ic (dashed magen a lines)
p essu e. O ange dash-do ed lines a e esul s in [63].
32

Figu e 16: No mal p o iles o he downslope no malized eloci y (10uα,i/(maxβ,j |uβ,j|),β= 1, ..., N, j = 1, ..., Nx)
ob ained wi h he NH-20 l. model o θ= 22◦and hi= 1.82 mm du ing g anula collapse a di e en posi ions, a (a)
= 0.60 s and (b) = 1.60 s. Dashed g een lines a e low/no- low in e aces.
(a) θ= 0◦,hi= 1 mm (b) θ= 22◦,hi= 1.82 mm
Model E (%) E (%) comp
( ini) ( in) (s)
NH-20 l. - - 1096.7 (19.08 min)
NH-10 l. 0.30 0.33 368.2 (6.1 min)
NH-8 l. 0.45 0.47 197.7 (3.3 min)
NH-5 l. 0.79 0.88 60.1 (1.0 min)
NH-2 l. 1.68 1.93 14.2
NH-1 l. 3.05 2.30 6.4
H-20 l. 10.01 8.32 56.1
H-10 l. 10.41 8.58 28.8
H-1 l. 11.65 11.25 4.7
Model E (%) E (%) comp
( ini) ( in) (s)
NH-20 l. - - 1517.4 (25.3 min)
NH-10 l. 0.92 1.92 544.6 (9.1 min)
NH-8 l. 1.31 2.87 290.6 (4.8 min)
NH-5 l. 2.19 5.58 89.2 (1.5 min)
NH-2 l. 4.70 14.18 21.0
NH-1 l. 6.47 22.80 9.7
H-20 l. 7.36 11.21 79.4
H-10 l. 8.23 12.90 41.0
H-1 l. 11.78 31.93 7.1
Table 9: l2-e o s made in he heigh by he di e en models compa ed o he e e ence solu ion co esponding o he
non-hyd os a ic model wi h 20 laye s a in e media e imes ini = 0.24 s(θ= 0◦),0.32 s(θ= 22◦), and a inal ime in
1.5s(θ= 0◦),3s(θ= 22◦). comp is he compu a ional ime needed o each in.
5.5.3 Quan i ica ion o model di e ences
By compa ing he NH-10 l., NH-1 l. and H-10 l. models, we ha e shown ha he p oposed NH-10 l. model could be seen
as a good comp omise be ween he NH-1 l. and H-10 l. models. Howe e , he p ice o pay is a highe compu a ional
e o . We ha e also concluded ha mul ilaye e ec s ha e mo e in luence o la ge slopes (a leas o cu en models),
while non-hyd os a ic e ec s a e mo e impo an o small slopes. Le us now quan i y he imp o emen made using
he mul ilaye model wi h a di e en numbe o laye s. In Table 9we gi e he di e ence be ween he mass hickness
33
0 5 10 15
0
0.5
1
1.5
0 5 10 15
0
0.5
1
1.5
0 5 10 15
0
0.5
1
1.5
0 5 10 15
0
0.5
1
1.5
0 5 10 15
0
0.5
1
1.5
Figu e 17: Time e olu ion o he no malized eloci y o he on compu ed wi h hyd os a ic (H) and non-hyd os a ic
(NH) models, and expe imen al da a (solid-ci cle blue lines) o hi= 0 mm. He e h0= 0.14 m, 0=√gcos θh0m/s and
τc=ph0/(gcos θ)s.
calcula ed wi h he hyd os a ic and non-hyd os a ic models a an in e media e ime and o he mass deposi o g anula
collapse on θ= 0◦and θ= 22◦, aking as e e ence he solu ion o he non-hyd os a ic model wi h 20 e ical laye s. We
see ha o θ= 0◦, i one uses 2 laye s in he non-hyd os a ic model he e o in he hickness is below 2%, whe eas we
need a leas 10 laye s o θ= 22◦ o ha e he same esul . I is clea hen ha he mul ilaye app oach is needed o la ge
slopes, bu i can be simpli ied o sho alues o θi one does no need addi ional in o ma ion on in-dep h a ia ions
o he low. One should no e ha his ac ( he in luence o mul ilaye e ec s) could a y i he app oxima ion o he
low/no- low in e ace is imp o ed, since i s size is cu en ly o e es ima ed (see Figu e 15). Imp o ing he app oxima ion
o he low/no- low in e ace would make he mul ilaye app oach mo e decisi e o small slope alues, whe e a wide s a ic
laye o ma e ial is expec ed o be eco e ed. Le us ema k ha we ha e also pe o med his quan i a i e analysis in he
case o a igid bed (hi= 0 mm) and wi h a smoo h ini ial condi ion ins ead o he ab up g anula collapse, ob aining
simila esul s o Table 9. No ice also ha hese di e ences shows he laye -independence o he esul s o N= 10 o
la ge and small slopes, ha is he numbe o laye s used in he compa isons wi h expe imen al da a.
Conce ning he compu a ional imes shown in Table 9, he non-hyd os a ic model is mo e expensi e han he
hyd os a ic one, as expec ed. Namely, we pay he p ize o sol ing he ellip ic ope a o associa ed o he non-hyd os a ic
p essu e, which is coupled o all cells and laye s, as seen in Sec ion 4and analyzed in Subsec ion 5.1.
6 Conclusions
In his wo k, we ha e in oduced a mul ilaye model including he µ(I)- heology and a weakly non-hyd os a ic p essu e
o simula e d y g anula lows. I is called a weakly non-hyd os a ic model, in he sense ha a e a dimensional analysis,
non-hyd os a ic i s -o de e ms ela ed o he heology a e neglec ed, whe eas second o de e ms ela ed o no mal
34
Figu e 18: G anula mass o e a plane o slope θ= 0◦(hi= 0 mm), 16◦(hi= 1.4mm), 22◦(hi= 1.82 mm), a
an in e media e ime and when he mass has s opped (deposi ) o he labo a o y expe imen s (solid-ci cle blue line), and
non-hyd os a ic models wi h and wi hou he sidewall ic ion e m o equa ion (36)(see Tables 1and 8). No e ha a
θ= 22◦, he g anula mass i is no a es ye .
accele a ion a e kep .
A well-balanced nume ical disc e iza ion is also p oposed, based on a h ee-s ep spli ing p ocedu e o deal wi h
iscosi y e ec s and non-hyd os a ic p essu e. The scheme is well-balanced o hose s eady s a es a es whe e he slope
o he ee su ace is lowe han he angle o epose o he ma e ial. This is achie ed hanks o a hyd os a ic econs uc ion
echnique in ol ing he Coulomb ic ion e m a he bo om.
Thus, he p oposed model can be seen as an ex ension o p e ious wo ks. Conc e ely, i gene alizes he hyd os a ic
mul ilaye model [36] o he non-hyd os a ic amewo k and he one-laye non-hyd os a ic model [42] o he mul ilaye
amewo k. Then, all he good esul s associa ed wi h he disc e iza ion in he di ec ion no mal o he slope in [36] a e
p esen in his model. In pa icula , we show he abili y o he model o ep oduce changes in he shape o he eloci y
p o iles and o app oxima e he low/no- low in e ace, which is no possible wi h one-laye models. In addi ion o he
mul ilaye app oach, he main achie emen s o he one-laye non-hyd os a ic model o [42] a e also p esen in his model.
In pa icula , he simula ion o he mass hickness a sho imes is imp o ed, and he accele a ion/decele a ion o he
on eloci y is eco e ed, con a y o hyd os a ic models whe e he i s accele a ion phase is no cap u ed.
Conce ning he in luence o he coo dina e sys em, we ha e ein o ced wo o he esul s gi en in [42]: (i)
non-hyd os a ic models a e less dependen on he coo dina e sys em han hyd os a ic models. Mo eo e , his dependence
is e en smalle when a no mal disc e iza ion is conside ed; (ii) he non-hyd os a ic model in Ca esian coo dina es
p oduces deposi s simila o he hyd os a ic model in local coo dina es, bu he simula ion di e a in e media e imes.
Fu he mo e, he di e ences be ween he deposi s a e smalle when he mul ilaye app oach is used.
Thus, on he one hand, he no el model combines he s eng hs o bo h p e ious models, sol ing some o hei
d awbacks. On he o he hand, he compu a ional cos o his model is highe han o he p e ious ones, mainly due
o he esolu ion o he linea sys em associa ed wi h he non-hyd os a ic p essu e, whose size is (M+ 1) ×N, being
M he numbe o cells and N he numbe o no mal laye s. In o de o speed up he scheme, we ha e p oposed a i s
p ecompu ing o he ini ial guess o he i e a i e sol e . This s a egy allows us o educe a ound 70% he compu a ional
ime when N= 20 laye s.
35
In summa y, he p oposed model could be a easonable comp omise be ween e iciency and accu acy, eco e ing he
good esul s o p e ious wo ks a he same ime. Howe e , he esul s he e sugges ha an impo an imp o emen o
hese models would be o conside he iscous con ibu ions o he momen um equa ion in he di ec ion no mal o he
slope, which a e neglec ed he e. I co esponds o sol ing he ull Na ie -S okes sys em in he laye -a e aged (mul ilaye )
amewo k. Tha is a di icul ask ha will be ackled in he u u e.
Acknowledgemen s
This wo k is pa ially suppo ed by g an s RTI2018-096064-B-C21 and RTI2018-096064-B-C22 unded by
MCIN/AEI/10.13039/501100011033 and “ERDF A way o making Eu ope”, and by p ojec s PID2020-114688RB-I00,
PID2022-137637NB-C21 and PID2022-137637NB-C22, he ERC con ac ERC-CG-2013-PE10-617472 SLIDEQUAKES,
he DT-GEO Digi al Twin p ojec and he En Seis Doc o al Ne wo k. C. Escalan e and J. Ga es-D´ıaz ha e been also
pa ially suppo ed by p ojec PID2020-114688RB-I00. J. Ga es-D´ıaz has been also pa ially suppo ed by he Eu opean
Union - Nex Gene a ionEU p og am.
0 5 10 15
0
0.5
1
1.5
0 5 10 15
0
0.5
1
1.5
0 5 10 15
0
0.5
1
1.5
0 5 10 15
0
0.5
1
1.5
0 5 10 15
0
0.5
1
1.5
0 5 10 15
0
0.5
1
1.5
Figu e 19: Time e olu ion o he no malized eloci y o he on compu ed wi h non-hyd os a ic models wi h and wi hou
he co ec ion in he µ(I)coe icien s accoun ing o he sidewall ic ion, and expe imen al da a (solid-ci cle blue lines),
o W= 10,20 cm. He e h0= 0.14 m, 0=√gcos θh0m/s and τc=ph0/(gcos θ)s.
36
A Coe icien s o he linea sys em associa ed o dispe si e e ec s
He e we desc ibe he coe icien s appea ing in (32) whe e we will assume ee-ou low bounda y condi ions o simplici y.
The e o e, he igh -hand side ec o in (33) is modi ied as ollows
R=


R1
.
.
.
RN


,Rα=



















0
0
CαUn+2/3
5/2,(∂XUn+2/3)5/2,(∂XZ)5/2
.
.
.
CαUn+2/3
Nx−5/2,(∂XUn+2/3)Nx−5/2,(∂XZ)Nx−5/2
0
0



















, α = 1, . . . , N.
The idiagonal ma ices Mα,α−1, Mα,α,and Mα,α+1 o α= 2, . . . , N in (32) a e desc ibed, o i, j = 3, . . . , Nx−2, by
(Mα,α−1)i,j =

























h2
i+1/2
∆x2
−hi+1/2
ξα,i+1/2+ξα+1/2,i+1/2
∆x+N2, o j=i−1,
−2
h2
i+1/2
∆x2+ 4ξ2
α,i+1/2+ 2hi+1/2∂xξα,i+1/2+ 2N2, o j=i,
h2
i+1/2
∆x2+hi+1/2
ξα,i+1/2+ξα+1/2,i+1/2
∆x+N2, o j=i+ 1,
(Mα,α)i,j =

























2
h2
i+1/2
∆x2
−hi+1/2
ξα,i+1/2+ξα+3/2,i+1/2−ξα+1,i+1/2−ξα−1/2,i+1/2
∆x
−2N2, o j=i−1,
−4
h2
i+1/2
∆x2
−4ξ2
α,i+1/2+ξ2
α+1,i+1/2−2hi+1/2∂xξα,i+1/2+∂xξα+1,i+1/2−4N2, o j=i,
2
h2
i+1/2
∆x2+hi+1/2
ξα,i+1/2+ξα+3/2,i+1/2−ξα+1,i+1/2−ξα−1/2,i+1/2
∆x
−2N2, o j=i+ 1
(Mα,α+1)i,j =

























h2
i+1/2
∆x2+hi+1/2
ξα+1,i+1/2+ξα+1/2,i+1/2
∆x+N2, o j=i−1,
−2
h2
i+1/2
∆x2+ 4ξ2
α+1,i+1/2
−2hi+1/2∂xξα+1,i+1/2+ 2N2, o j=i,
h2
i+1/2
∆x2
−hi+1/2
ξα+1,i+1/2+ξα+1/2,i+1/2
∆x+N2, o j=i+ 1,
wi h MN,N+1 =0,and
ξα,i+1/2= (α−1/2)hi+1 −hi
∆x+Nh+
i+1/2−h−
i+1/2−(hi+1 −hi)
∆x,
ha is p opo ional o he de i a i e o he midpoin o he laye Zα,and we app oxima e hei second o de de i a i es
as ollows
∂xξα,1+1/2=1
∆xminmod ξα,i+3/2−ξα,i+1/2, ξα,i+1/2−ξα,i−1/2,1
2ξα,i+3/2−ξα,i−1/2,
minmod being he usual slope limi e unc ion.
37

The special cases o M1,1,and M1,2,co esponds o he disc e iza ion a he lowes laye , whe e he cons ain (27) is
simply gi en by 2hNN1= 2hNE1,and eads
(M1,1)i,j =

























h2
i+1/2
∆x2
−hi+1/2
ξα+3/2,i+1/2−ξα+1,i+1/2
∆x
−N2, o j=i−1,
−2
h2
i+1/2
∆x2
−4ξ2
α+1,i+1/2+ 2hi+1/2∂xξα+1,i+1/2−2N2, o j=i,
h2
i+1/2
∆x2+hi+1/2
ξα+3/2,i+1/2−ξα+1,i+1/2
∆x
−N2, o j=i+ 1
(M1,2)i,j =

























h2
i+1/2
∆x2+hi+1/2
ξα+1,i+1/2+ξα+1/2,i+1/2
∆x+N2, o j=i−1,
−2
h2
i+1/2
∆x2+ 2N2, o j=i,
h2
i+1/2
∆x2
−hi+1/2
ξα+1,i+1/2+ξα+1/2,i+1/2
∆x+N2, o j=i+ 1.
Rega ding he ho izon al wall/ ee-ou low bounda y condi ions along he independen a iable x, we conside , o
α= 1,2,...N,
(Mα,α)1,1= (Mα,α)2,2= (Mα,α)Nx−1,Nx−1= (Mα,α)Nx,Nx= 1,
(Mα,α)1,2= (Mα,α)2,3= (Mα,α)Nx−1,Nx−2= (Mα,α)Nx,Nx−1=−1,(37)
and o he es o ma ices, Mα,α±1, hey a e gi en as usually when applying ghos cells echniques o all he a iables.
B P oo o Theo em 1
This appendix gi es de ailed p oo o Theo em 1conce ning he ene gy balance o model (9). Some o he s eps o his
p oo a e simila o he p oo o he ene gy inequali y e i ied by LDNH0model in [40]. Howe e , we gi e he e all he
de ails o he sake o comple eness.
P oo :
Using he ho izon al momen um equa ion (9b) as well as he laye heigh e olu ion equa ion (9a), we ob ain
hα∂ uα+∂xu2
α
2+uα(Gα+1/2−Gα−1/2) + ∂xhαpnh
α
ρ+gcos θhα∂x(eb+b+h)
=1
ρKα−1/2−Kα+1/2+pnh
α+1/2
ρ∂xzα+1/2−pnh
α−1/2
ρ∂xzα−1/2+uα+1/2Gα+1/2−uα−1/2Gα−1/2.
(38)
Now, we sum up (38) mul iplied by uαwi h he equa ion ob ained by mul iplying he mass conse a ion equa ion by
u2
α/2 + gcos θ(eb+b+h), ge ing
∂ hα
u2
α
2+gcos θ(eb+b+h)∂ hα+∂xhαuαu2
α
2+gcos θ(eb+b+h)
−gcos θ(eb+b+h)Gα+1/2−Gα−1/2+1
ρuα∂xhαpnh
α−1
ρuαpnh
α+1/2∂xzα+1/2−pnh
α−1/2∂xzα−1/2
−Gα+1/2u2
α
2−uαuα+1/2+Gα−1/2u2
α
2−uαuα−1/2+1
ρuαKα−1/2−Kα+1/2
(39)
Simila s eps a e ollowed o he equa ions o he e ical eloci y (9c), and using he mass equa ion, we ge
∂ hα
w2
α
2+∂xhαuα
w2
α
2+1
ρwαpnh
α+1/2−pnh
α−1/2=
−Gα+1/2w2
α
2−wαewα+1/2+Gα−1/2w2
α
2−wαewα−1/2.(40)
38
By de ining
Eα=hαu2
α+w2
α
2+gcos θeb+b+h
2,
and summing up equa ions (39) and (40), a e some s aigh o wa d compu a ions we ob ain
∂ Eα+∂xEαuα+gcos θhα
h
2uα+1
ρPNH,α =MTα+1
ρuαKα−1/2−Kα+1/2+gcos θ
2(hα∂ h−h∂ hα),(41)
whe e we ha e used ha
gcos θ(eb+b+h)∂ hα=∂ hαgcos θeb+b+h
2+gcos θ
2(h∂ hα−hα∂ h).
In equa ion (41), PNH,α and MTαcollec he e ms ela ed o non-hyd os a ic p essu e and mass ans e ence, espec i ely,
in he laye Ωα. The non-hyd os a ic con ibu ions a e
PNH,α =uα∂xhαpnh
α+pnh
α+1/2−uα∂xzα+1/2+wα−pnh
α−1/2−uα∂xzα−1/2+wα
=∂xhαuαpnh
α−pnh
α+1/2
α
X
β=1
∂x(hβuβ) + pnh
α−1/2
α−1
X
β=1
∂x(hβuβ)
whe e we ha e used es ic ion (4) and (9d). No ice ha in he pa icula cases α= 1 and α=N, using he
non-pene a ion condi ion and pnh
|b+h= 0, we ha e
PNH,1=∂xh1u1pnh
1−pnh
3/2∂x(h1u1), PNH,N =∂xhNuNpnh
N+pnh
N−1/2
N−1
X
β=1
∂x(hβuβ).
I leads o a conse a i e con ibu ion since he e ms on pnh
α±1/2 anish when summing up om α= 1, . . . , N:
N
X
β=1
PNH,β =
N
X
β=1
∂xhβpnh
βuβ.
Conce ning he mass ans e e ms, hey a e
MTα=gcos θeb+b+hGα+1/2−Gα−1/2+Gα+1/2
2uαuα+1 +wαwα+1−Gα−1/2
2uαuα−1+wαwα−1.
I is easy o see ha , since we ha e se G1/2=GN+1/2= 0, i holds
N
X
β=1
MTβ= 0,and g
2
N
X
β=1
(hα∂ H−H∂ hα)=0,
and hus, summing up (41) in all he laye s, he p oo is inished.
Finally, by using he de ini ions o Kα+1/2 o α= 0, . . . , N, he iscous e ms sa is y
N
X
β=1
uαKβ−1/2−Kβ+1/2=u1K1/2−uNKN+1/2+
N−1
X
β=1
Kβ+1/2(uβ+1 −uβ)
=−µ(I1/2)p1/2
u2
1
|u1|−
N−1
X
β=1
ηβ+1/2
2
(uβ+1 −uβ)2
hβ+1/2
,
ha is hen a dissipa i e con ibu ion, as expec ed.
□
39
Re e ences
[1] B. And eo i, Y. Fo e e, and O. Pouliquen. G anula Media: Be ween Fluid and Solid. Camb idge Uni e si y
P ess, 2013. 1
[2] E. Audusse, M-O. B is eau, and A. Decoene. Nume ical simula ions o 3D ee su ace lows by a mul ilaye
Sain -Venan model. In e na ional Jou nal o Nume ical Me hods in Fluids, 56(3):331–350, 2008. 8
[3] E. Audusse, M.-O. B is eau, B. Pe hame, and J. Sain e-Ma ie. A mul ilaye Sain -Venan sys em wi h mass
exchanges o shallow wa e lows. De i a ion and nume ical alida ion. ESAIM: Ma hema ical Modelling and
Nume ical Analysis, 45(1):169–200, jun 2010. 2,11
[4] N. A¨ıssiouene, M.-O. B is eau, E. Godlewski, A. Mangeney, C. Pa ´es Mad o˜nal, and J. Sain e-Ma ie. A
wo-dimensional me hod o a amily o dispe si e shallow wa e models. The SMAI jou nal o compu a ional
ma hema ics, 6:187–226, Sep embe 2020. 3
[5] T. Ba ke and J. M. N. T. G ay. Pa ial egula isa ion o he incomp essible µ(I)- heology o g anula low. Jou nal
o Fluid Mechanics, 828:5–32, aug 2017. 2
[6] T. Ba ke , M. Rau e , E. S. F. Magui e, C. G. Johnson, and J. M. N. T. G ay. Coupling heology and seg ega ion
in g anula lows. Jou nal o Fluid Mechanics, 909, dec 2020. 1,2
[7] T. Ba ke , D. G. Schae e , P. Boho quez, and J. M. N. T. G ay. Well-posed and ill-posed beha iou o he
µ(I)- heology o g anula low. Jou nal o Fluid Mechanics, 779:794–818, 9 2015. 1
[8] T. Ba ke , D. G. Schae e , M. Shea e , and J. M. N. T. G ay. Well-posed con inuum equa ions o g anula low
wi h comp essibili y and µ(I)- heology. P oceedings o he Royal Socie y A: Ma hema ical, Physical and Enginee ing
Sciences, 473(2201):20160846, may 2017. 2
[9] F. Bouchu . Nonlinea S abili y o Fini e Volume Me hods o Hype bolic Conse a ion Laws. Bi kh¨ause Basel,
2004. 14,18
[10] F. Bouchu , J. M. Delgado-S´anchez, E. D. Fe n´andez-Nie o, A. Mangeney, and G. Na bona-Reina. A bed
p essu e co ec ion o he ic ion e m o dep h-a e aged g anula low models. Applied Ma hema ical Modelling,
106:627–658, jun 2022. 2,28
[11] F. Bouchu , E. D. Fe n´andez-Nie o, E. H. Kon´e, A. Mangeney, and G. Na bona-Reina. Dila ancy in d y g anula
lows wi h a comp essible µ(i) heology. Jou nal o Compu a ional Physics, 429:110013, ma 2021. 2
[12] F. Bouchu , E. D. Fe n´andez-Nie o, A. Mangeney, and G. Na bona-Reina. A wo-phase wo-laye model o luidized
g anula lows wi h dila ancy e ec s. Jou nal o Fluid Mechanics, 801:166–221, Augus 2016. 1,2
[13] F. Bouchu , I. Ionescu, and A. Mangeney. An analy ic app oach o he e olu ion o he s a ic/ lowing in e ace in
iscoplas ic g anula lows. Communica ions in Ma hema ical Sciences, 14(8):2101–2126, 2016. 3
[14] F. Bouchu and M. Wes dickenbe g. G a i y d i en shallow wa e models o a bi a y opog aphy. Communica ions
in Ma hema ical Sciences, 2(3):359–389, 09 2004. 24,26
[15] J. Boussinesq. Th´eo ie des ondes e des emous qui se p opagen le long d’un canal ec angulai e ho izon al, en
communiquan au liquide con enu dans ce canal des i esses sensiblemen pa eilles de la su ace au ond. J. Ma h.
Pu es Appl., pages 55–108, 1872. 2
[16] M. Bouzid, M. T ulsson, P. Claudin, E. Cl´emen , and B. And eo i. Nonlocal heology o g anula lows ac oss yield
condi ions. Physical Re iew Le e s, 111(23), dec 2013. 2
[17] F. Boye , ´
E. Guazzelli, and O. Pouliquen. Uni ying suspension and g anula heology. Physical Re iew Le e s,
107(18), oc 2011. 2
[18] M.-O. B is eau, C. Guicha d, B. di Ma ino, and J. Sain e-Ma ie. Laye -a e aged eule and Na ie -S okes equa ions.
Communica ions in Ma hema ical Sciences, 15(5):1221–1246, 2017. 9
40
[19] M. B une , L. Mo e i, A. Le F ian , A. Mangeney, E. D. Fe n´andez Nie o, and F. Bouchu . Nume ical simula ion o
he 30–45 ka deb is a alanche low o Mon agne Pel´ee olcano, Ma inique: om olcano lank collapse o subma ine
emplacemen . Na u al Haza ds, 87(2):1189–1222, ap 2017. 2
[20] S. Bus o, M. Dumbse , C. Escalan e, N. Fa ie, and S. Ga ilyuk. On high o de ADER Discon inuous Gale kin
schemes o i s o de hype bolic e o mula ions o nonlinea dispe si e sys ems. Jou nal o Scien i ic Compu ing,
87(2), Ma ch 2021. 3
[21] M. J. Cas o, J. M. Gonz´alez-Vida, and C. Pa ´es. Nume ical ea men o we /d y on s in shallow lows wi h a
modi ied Roe scheme. Ma hema ical Models and Me hods in Applied Sciences, 16(06):897–931, jun 2006. 15,18
[22] M. J. Cas o D´ıaz and E. D. Fe n´andez-Nie o. A class o compu a ionally as i s o de ini e olume sol e s: PVM
me hods. SIAM Jou nal on Scien i ic Compu ing, 34(4):A2173–A2196, jan 2012. 14,18
[23] V. Casulli. A semi-implici ini e di e ence me hod o non-hyd os a ic ee-su ace lows. In e na ional Jou nal o
Nume ical Me hods in Fluids, 30(4):425–440, jun 1999. 2
[24] V. Casulli and P. Zanolli. Semi-implici nume ical modeling o nonhyd os a ic ee-su ace lows o en i onmen al
p oblems. Ma hema ical and Compu e Modelling, 36(9-10):1131–1149, dec 2002. 3
[25] J. Chaucha and M. M´edale. A h ee-dimensional nume ical model o dense g anula lows based on he µ(I)- heology.
Jou nal o Compu a ional Physics, 256(0):696–712, 2014. 19
[26] A. J. Cho in. Nume ical solu ion o he Na ie -S okes equa ions. Ma hema ics o Compu a ion, 22(104):745–762,
1968. 3
[27] A. Dedne , F. Kemm, D. K ¨one , C.-D. Munz, T. Schni ze , and M. Wesenbe g. Hype bolic di e gence cleaning o
he MHD equa ions. Jou nal o Compu a ional Physics, 175(2):645–673, Janua y 2002. 3
[28] R Delannay, A Valance, A Mangeney, O Roche, and P Richa d. G anula and pa icle-laden lows: om labo a o y
expe imen s o ield obse a ions. Jou nal o Physics D: Applied Physics, 50(5):053001, 2017. 1
[29] R. P. Denlinge and R. M. I e son. G anula a alanches ac oss i egula h ee-dimensional e ain: 1. Theo y and
compu a ion. Jou nal o Geophysical Resea ch: Ea h Su ace, 109(F1), ma 2004. 3,25
[30] C. Escalan e, T. Mo ales de Luna de Luna, and M. J. Cas o. Non-hyd os a ic p essu e shallow lows: GPU
implemen a ion using ini e olume and ini e di e ence scheme. Applied Ma hema ics and Compu a ion,
338(338):631–659, dec 2018. 3
[31] C. Escalan e, E. D. Fe n´andez-Nie o, T. Mo ales de Luna, and M. J. Cas o. An e icien wo-laye non-hyd os a ic
app oach o dispe si e wa e wa es. Jou nal o Scien i ic Compu ing, 79(1):273–320, oc 2018. 18
[32] C. Escalan e and T. Mo ales de Luna. A gene al non-hyd os a ic hype bolic o mula ion o Boussinesq dispe si e
shallow lows and i s nume ical app oxima ion. Jou nal o Scien i ic Compu ing, 83(3), jun 2020. 3
[33] C. Escalan e S´anchez, E. D. Fe n´andez-Nie o, T. Mo ales de Luna, Y. Penel, and J. Sain e-Ma ie. Nume ical
simula ions o a dispe si e model app oxima ing ee-su ace Eule equa ions. Jou nal o Scien i ic Compu ing,
89(3), oc 2021. 11,18
[34] N. Fa ie and S. Ga ilyuk. A apid nume ical me hod o sol ing se e–g een–naghdi equa ions desc ibing long ee
su ace g a i y wa es. Nonlinea i y, 30(7):2718–2736, May 2017. 3
[35] E. D. Fe n´andez-Nie o, F. Bouchu , D. B esch, M. J. Cas o D´ıaz, and A. Mangeney. A new Sa age-Hu e ype
model o subma ine a alanches and gene a ed sunami. Jou nal o Compu a ional Physics, 227(16):7720–7754, 2008.
26
[36] E. D. Fe n´andez-Nie o, J. Ga es-D´ıaz, A. Mangeney, and G. Na bona-Reina. A mul ilaye shallow model o d y
g anula lows wi h he µ(I)- heology: applica ion o g anula collapse on e odible beds. Jou nal o Fluid Mechanics,
798:643–681, jun 2016. 1,2,3,4,5,6,7,11,19,28,29,35
[37] E. D. Fe n´andez-Nie o, J. Ga es-D´ıaz, A. Mangeney, and G. Na bona-Reina. 2D g anula lows wi h he µ(I)
heology and side walls ic ion: a well-balanced mul ilaye disc e iza ion. Jou nal o Compu a ional Physics,
356:192–219, 2018. 2,14,15,31
41