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

Appl. Nume . Ma h. 215 (2025) 138–156
A ailable online 17 Ap il 2025
0168-9274/© 2025 The Au ho (s). Published by Else ie B.V. on behal o IMACS. This is an open access a icle unde he CC BY-NC-ND license
(h p://c ea i ecommons.o g/licenses/by-nc-nd/4.0/).
Con en s lis s a ailable a ScienceDi ec
Applied Nume ical Ma hema ics
jou nal homepage: www.else ie .com/loca e/apnum
Resea ch Pape
High-o de well-balanced schemes o shallow models o d y
a alanches
M.J. Cas o Díaza, C. Escalan eb,, J. Ga es-Díazc, ,∗, T. Mo ales de Lunaa,
aDp o. Análisis Ma emá ico, Es ad. e I.O. y Ma emá ica Aplicada, Uni e sidad de Málaga, 29071 Málaga, Spain
bDp o. Ma emá ica Aplicada, Uni e sidad de Málaga, 29071 Málaga, Spain
cDp o. Ma emá ica Aplicada II, Uni e sidad de Se illa, 41092 Se illa, Spain
A R T I C L E I N F O A B S T R A C T
Keywo ds:
Fini e olume me hod
Well-balanced schemes
High-o de scheme
Sa age-Hu e model
G anula flows
In his wo k we conside a dep h-a e aged model o g anula flows wi h a Coulomb- ype ic ion
o ce desc ibed by he 𝜇(𝐼) heology. In his model, he so-called lake-a - es s eady s a es a e
o special in e es , whe e eloci y is ze o and he slope is unde a c i ical h eshold defined by
he angle o epose o he g anula ma e ial. I leads o a amily wi h an infini e numbe o lake-
a - es s eady s a es. We desc ibe a well-balanced econs uc ion p ocedu e ha allows o define
well-balanced fini e olume me hods o such p oblem. The echnique is gene alized o high-
o de space/ ime schemes. In pa icula , he second and hi d-o de schemes a e conside ed in
he nume ical es s sec ion. An accu acy es is included showing ha second and hi d-o de a e
achie ed. A well-balanced es is also conside ed. The p oposed scheme is well-balanced o s eady
s a es wi h non-cons an ee su ace, and i is exac ly well-balanced o hose s eady s a es gi en
by a simple cha ac e iza ion.
1. In oduc ion
Geophysical flows ha e been deeply s udied in las decades due o hei impo ance designing ea ly wa ning sys ems o s udies
agains na u al haza ds (floods, landslides, sunamis, snow a alanches, olcanic e up ions, e c.). The ma hema ical modeling and nu-
me ical simula ion o such flows is a e y ac i e esea ch opic, especially since popula iza ion o ad anced echnological esou ces,
such as high-pe o mance compu e s. Ne e heless, many challenges emain: he physical unde s anding and he design o complex
heologies; he de elopmen o sophis ica ed models whose associa ed compu a ional cos is easonable o p ac ical pu poses; he
design o efficien nume ical schemes o sol e complex models including iscous ( iscoplas ic, iscoelas ic,...) e ms. The complexi y
o hese flows equi es p og ess on all hese esea ch lines, whe e diffe en scien is s (geophysicis s, ma hema icians, compu e sci-
en is s,...) may con ibu e wi h aluable ad ances. He e, we ocus on g anula flows, conc e ely ae ial a alanches, and we ackle he
ask o defining efficien nume ical schemes o some dep h-a e aged a alanche models p e iously in oduced in he li e a u e.
Conce ning he modeling o ae ial a alanches, one o he mos accep ed heological laws desc ibing hei dynamics is he so-called
𝜇(𝐼)- heology (see [1]). This heology conside s a a iable ic ion coefficien depending on he eloci y and p essu e. I is based on
a D ucke -P age plas ici y c i e ion defining he de ia o ic enso as
*Co esponding au ho .
E-mail add ess: [email p o ec ed] (J. Ga es-Díaz).
h ps://doi.o g/10.1016/j.apnum.2025.04.008
Recei ed 17 Sep embe 2024; Recei ed in e ised o m 17 Ma ch 2025; Accep ed 15 Ap il 2025
Applied Nume ical Ma hema ics 215 (2025) 138–156
139
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
⎧
⎪
⎨
⎪
⎩
𝝉=𝜇(𝐼)𝑝
‖𝑫‖
𝑫i ‖𝑫‖≠0,
‖𝝉‖≤𝜇𝑠𝑝i ‖𝑫‖=0,
(1)
whe e 𝑫is he s ain a e enso , 𝑝is he p essu e and 𝜇𝑠= an𝜃0, wi h 𝜃0 he angle o epose o he g anula ma e ial. This heology
was conside ed in [2] whe e i was included in a con inuous 2D Na ie -S okes sol e and he g anula ma e ial was conside ed as a
fluid wi h a 𝜇(𝐼)- iscosi y. In o de o a oid he singula i y when he s ain a e anishes, he au ho s conside ed a egula iza ion o
he 𝜇(𝐼)- iscosi y. The egula iza ion echnique was compa ed wi h he Augmen ed Lag angian (AL) me hod in [3] o hese kind o
flows. I was shown ha , o a p ope choice o he egula iza ion pa ame e and wi h a sui able fini e elemen disc e iza ion, ha he
e o associa ed o his egula iza ion is negligible wi h espec o he ime/space disc e iza ion e o s. O he au ho s ha e conside ed
fini e elemen disc e iza ions combined wi h AL me hods o simula e g anula flows and analyze diffe en aspec s associa ed o his
heology, including compa isons wi h labo a o y expe imen s (see e.g. [4,5]). Howe e , he compu a ional cos associa ed o hese
disc e iza ions is e y high. The e o e, simplified models, based on a dep h-a e aging p ocess, ha e become e y popula o hese
flows. This ype o models will be he s a ing poin o his wo k. Le us ema k ha some o he mo e complex models ha e been
p oposed in he las decade o g anula flows. Fo ins ance, in [6–9] one can find models wi h non-local effec s o a comp essible
e sion o he 𝜇(𝐼)- heology, whe e he solid olume ac ion is no mo e cons an .
As p e iously said, dep h-a e aged models ha e been widely used o s udy g anula flows. Since he pionee ing wo k o Sa age-
Hu e [10], whe e he ic ion o ce wi h he bo om is modeled h ough a Coulomb ic ion e m, se e al mo e sophis ica ed dep h-
a e aged models ha e been in oduced. Fo ins ance, G ay and Edwa ds [11] p oposed a shallow model, which includes second-
o de iscous e ms, o d y g anula flows wi h 𝜇(𝐼)- heology (see also [12]). An ex ension o his model, neglec ing second-o de
iscous e ms and using he mul ilaye (o laye -a e aged) amewo k, has been p oposed in [13]. Mo eo e , dispe si e effec s can
be inco po a ed as in [14], whe e a weakly non-hyd os a ic shallow model was p oposed, which akes in o accoun he con ibu ion
o he e ical accele a ion o he flow. This model was also ex ended o he mul ilaye case in [15].
Mos o hese models we e ob ained using a il ed e e ence sys em o e a plane wi h cons an slope and pe o ming an a e age on
he no mal di ec ion. This means ha in hese models he a e age ho izon al componen o he eloci y is measu ed in he downslope
di ec ion, which is an ad an age wi h espec o Ca esian models. Mo eo e , when alida ing such models, some o he labo a o y
expe imen s used can be only ep oduced o hese local models due o he defini ion o hei ini ial condi ions. In e es ingly, in [16]
au ho s showed ha he mo ion c i e ion is w ong in all p e ious local models due o he change o a iables ( o il ed coo dina es),
and hey p oposed a co ec ion o p ope ly cap u e his mo ion c i e ion.
F om he nume ical poin o iew, one o he majo difficul ies o disc e izing hese models is he nume ical ea men o he
mul i-e alua ed s ess enso (1)when he eloci y anishes (case 𝑫=0). On he one hand, in he mo ion phase o he flow (|𝑢|>0),
(semi-)implici disc e iza ions o he iscous e m ha e been conside ed o ic ion e ms in o de o a oid es ic i e ime-s eps (see
e.g. [17–19]). On he o he hand, diffe en ea men s ha e been conside ed o deal wi h he mul i-e alua ed defini ion o he s ess
enso when he flow is a es (𝑢=0). In [13,15] a egula iza ion echnique was conside ed, which makes he eloci y end o ze o in
s op-mo ion si ua ions, al hough null eloci y canno be ob ained, since a esidual eloci y emains depending on he egula iza ion
pa ame e . In [14,20] a diffe en ea men is used based on he in e p e a ion o he ic ion o ce made in [21]. I assumes ha he
ic ion o ce opposes o he mo ion o he flow, bu i canno change he sign o he eloci y. This assump ion means ha he ic ion
o ce should be bounded in p ac ice. Al hough i is difficul o apply his echnique o mul ilaye sys ems, i is a alid app oach
o single-laye models. I has he ad an age ha a ze o eloci y is exac ly eco e ed when he ma e ial s ops and emains a es .
Mo eo e , his leads o an explici disc e iza ion o he model wi h a pos e io i ( ic ion) co ec ion o he eloci y. In pa icula , a
simila app oach will be used in his wo k.
An impo an p ope y conce ning nume ical schemes o geophysical flows is he so-called well-balanced p ope y. Indeed, many
p ac ical applica ions o fluid o g anula flows co espond o small pe u ba ion o s eady s a es. In such si ua ions, he use o a
well-balanced nume ical me hod is manda o y since he ampli udes may be o he o de (o bigge han) he unca ion e o o he
scheme. This is o ins ance he case o a d y a alanche p opaga ing o e a slope which is ini ially s eady o when i con e ges o a
s eady s a e a e a la ge ime e olu ion. Rema k ha i is no always possible o efine he mesh so ha he unca ion e o o he
me hod is lowe han his pe u ba ion. The e o e, i is impo an o p ese e, all o a sub amily o , he s eady s a es o he model.
See o ins ance he wo ks [22–28] o balance laws, in pa icula , shallow wa e sys ems, and [29,20,30–34] o a alanche and
g anula models, among many o he s.
Focusing on g anula flows models, which is he opic o his wo k, he so called lake-a - es a e o special in e es . These equilib ia
co espond o si ua ions whe e he eloci y is ze o and he non-fla ee su ace o he ma e ial has a slope smalle han a c i ical alue
fixed by he epose angle. Rema k ha , in gene al, such s eady s a es a e cha ac e ized by a inequali y ins ead o an equali y. This
ac makes i difficul o design well-balanced nume ical schemes p ese ing such s eady solu ions. In p e ious wo ks ([31,14,15])
fi s -o de well-balanced schemes we e p oposed. The well-balanced p ope y was achie ed by means o a modified hyd os a ic
econs uc ion which includes he Coulomb ic ion e m (see [35]). In [20] a second-o de MUSCL econs uc ion combined wi h
cen al-upwind scheme was p oposed o a one-dimensional Sa age-Hu e ype model o subma ine landslides and gene a ed sunami
wa es. The Coulomb ic ion e m was ea ed using [29] and he p oposed scheme is well-balanced o lake-a - es s eady s a es
Howe e , designing high-o de well-balanced schemes is no a i ial ask and, up o ou knowledge, he e a e no p e ious wo ks
showing a high o de well-balanced scheme o g anula o complex flows including such a ic ion coefficien . Highe -o de schemes
cap u e g adien s and a ia ions in he solu ion mo e accu a ely han low-o de me hods. They educe nume ical dissipa ion and
dispe sion e o s, which a e c ucial o accu a ely cap u ing shocks, discon inui ies, and smoo h solu ion ea u es. Mo eo e , o a
Applied Nume ical Ma hema ics 215 (2025) 138–156
140
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
gi en le el o accu acy, high-o de schemes o en equi e ewe g id poin s compa ed o low-o de me hods. The e o e, high-o de
schemes allow o educe compu a ional cos s, especially in la ge-scale simula ions, which esul s in a be e efficiency.
The goal o his wo k is o p opose a high o de well-balanced disc e iza ion o shallow models o g anula flows including he
mul i-e alua ed defini ion o he s ess enso when he flow is a es . Based on he ideas in [36], we de elop a no el p ocedu e o
ob ain a well-balanced scheme o lake-a - es solu ions his kind o g anula flows. As al eady men ioned, he main difficul y in his
case lies in he exp ession o he s a iona y solu ions, which a e defined in e ms o an inequali y. The main no el y o his pape is o
desc ibe a p ocedu e o ob ain a well-balanced scheme o g anula flows, which can be easily gene alized o he a bi a y high-o de
case.
The pape is o ganized as ollows: Sec ion 2is de o ed o he desc ip ion o he model and o in oduce he basic no a ion used
in he pape . Sec ion 3will desc ibe he nume ical scheme p oposed he e. In pa icula , in Subsec ion 3.1 he gene al p ocedu e
o ob ain a well-balanced scheme based on a econs uc ion ope a o is p esen ed, and Subsec ion 3.2 is de o ed o desc ibe he
econs uc ion p ocedu e o he g anula model. Subsec ion 3.3 deals wi h he ime disc e iza ion. Some nume ical es s will be
p esen ed in Sec ion 4. In pa icula , we will show ha he expec ed o de o accu acy is achie ed o he second and hi d-o de
case, and also ha he scheme is indeed well-balanced. Finally, he conclusions a e p esen ed in Sec ion 5. The pape is comple ed
wi h Appendix A, whe e some pa icula cases o s eady s a es a e s udied.
2. Shallow model o g anula flows
Fo he sake o simplici y and cla i y in he exposi ion, we ocus he e on he hyd os a ic model in oduced in [14] o g anula
flows, which coincides in ac wi h he one in [11] when he second-o de iscous e m is neglec ed. We will conside Ca esian
coo dina es, al hough e e y hing can be i ially ex ended o he il ed coo dina es case. Adding non-hyd os a ic effec s is no majo
p oblem ei he , ollowing [14].
Thus, we conside a bidimensional domain, (𝑥,𝑧)∈ℝ2, and deno e by ℎ(𝑥),𝑢(𝑥)∈ℝ he o al fluid dep h and he dep h-in eg a ed
ho izon al eloci y, espec i ely. Fo simplici y, we conside ha he bo om opog aphy, 𝑏(𝑥), is a con inuous and diffe en iable
unc ion. The hyd os a ic model o g anula flows is hen w i en as
𝜕𝑡ℎ+𝜕𝑥(ℎ𝑢)=0,
𝜕𝑡(ℎ𝑢)+𝜕𝑥(ℎ𝑢2+1
2𝑔ℎ2)+𝑔ℎ𝜕𝑥𝑏=−𝜏𝑥𝑧|𝑏∕𝜌,
whe e 𝜌is he densi y o he ma e ial, assumed o be cons an , 𝜏𝑥𝑧 is he de ia o ic s ess enso , which includes he ic ion wi h he
bo om, defined by
𝜏𝑥𝑧|𝑏=⎧
⎪
⎨
⎪
⎩
𝜌𝑔ℎ𝜇|𝑏
𝑢
|𝑢|i |𝑢|≠0,
|||𝜏𝑥𝑧|𝑏|||≤𝜌𝑔ℎ𝜇𝑠i |𝑢|=0,
(2a)
whe e 𝜇is he a iable ic ion coefficien
𝜇|𝑏=𝜇𝑠+𝜇2−𝜇𝑠
𝐼0+𝐼|𝑏
𝐼|𝑏, wi h 𝐼|𝑏=
𝑑𝑠|||(𝜕𝑧𝑢)|𝑏|||
√𝜑𝑠𝑔ℎ , and (𝜕𝑧𝑢)|𝑏=𝑢
ℎ,(2b)
being 𝑑𝑠 he mean g ain size and 𝜇2,𝜇𝑠,𝐼0,𝜑𝑠cons an heological pa ame e s depending on he conside ed g anula ma e ial.
When looking o he s eady solu ions o sys em (2), we ocus on lake-a - es solu ions, ha is, hose wi h ze o eloci y. Then,
looking a he momen um conse a ion equa ion wi h 𝑢=0, and deno ing by 𝜂=𝑏+ℎ he ee su ace le el, we ge
𝑔ℎ||𝜕𝑥𝜂||=|||𝜏𝑥𝑧|𝑏,𝑢=0 |||
𝜌
≤𝑔ℎ𝜇𝑠.
The e o e, he s eady s a es co esponding o wa e a es , a e gi en by
𝑢=0, and ||𝜕𝑥𝜂||≤𝜇𝑠,(3)
whe e we ecall ha 𝜇𝑠= an𝜃0, wi h 𝜃0 he angle o epose o he ma e ial.
We will now p oceed o design a high o de nume ical scheme o (2) ha p ese es lake-a - es s eady s a es in a sense o be
p ecised.
3. Nume ical scheme
Fi s o all, le us w i e sys em (2)in a compac o m
𝜕𝑡𝑼+𝜕𝑥𝑭(𝑼)+𝑺(𝑼)𝜕𝑥𝑏=𝑻(𝑼),(4)
whe e
Applied Nume ical Ma hema ics 215 (2025) 138–156
141
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
𝑼=(ℎ
ℎ𝑢), 𝑭(𝑼)=(ℎ𝑢
ℎ𝑢2+1
2𝑔ℎ2), 𝑺(𝑼)=(0
𝑔ℎ), 𝑻(𝑼)=(0
−𝜏𝑥𝑧|𝑏∕𝜌).
Recall ha we a e assuming he e a con inuous and diffe en iable opog aphy 𝑏(𝑥). The ex ension o he case o a discon inuous bo om
could be add essed ollowing he ideas o [36,37].
We conside now, as usual in fini e olume schemes, a spa ial disc e iza ion in o con ol olumes o cells 𝑉𝑖=[𝑥𝑖−1∕2,𝑥𝑖+1∕2], o
𝑖∈={1,...,𝑁}. Fo he sake o simplici y, we conside a uni o m disc e iza ion wi h cell sizes Δ𝑥=𝑥𝑖+1∕2 −𝑥𝑖−1∕2. Le us deno e
by 𝑥𝑖=(𝑥𝑖−1∕2 +𝑥𝑖+1∕2)∕2 he cen e o 𝑉𝑖. We keep, o now, he ime a iable con inuous and deno e by 𝑼𝑖(𝑡) he cell a e age o
𝑼(𝑥,𝑡), ha is, o 𝑖∈:
𝑼𝑖(𝑡)= 1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2
𝑼(𝑥,𝑡) 𝑑𝑥. (5)
Following [38,28,37], gi en a sequence o s a es {𝑼𝑖(𝑡)}𝑖∈, we conside , o e e y cell 𝑉𝑖, a 𝑠-o de econs uc ion ope a o
𝑷𝑡
𝑖(𝑥)=𝑷𝑖(𝑥;{𝑼𝑗(𝑡)}𝑗∈𝑖),
whe e 𝑗is he s encil o cells o dependence o he ope a o , sa is ying
1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2
𝑷𝑡
𝑖(𝑥) 𝑑𝑥 =𝑼𝑖(𝑡)=𝑼(𝑥𝑖,𝑡)+(Δ𝑥𝑠).
We hen conside he econs uc ed alues a he in e -cells
𝑼𝑡,+
𝑖−1∕2 = lim
𝑥→𝑥+
𝑖−1∕2
𝑷𝑡
𝑖(𝑥), 𝑼−
𝑖+1∕2 = lim
𝑥→𝑥−
𝑖+1∕2
𝑷𝑡
𝑖(𝑥),
which a e assumed o be app oxima ions o o de 𝑠o he solu ion a he in e aces:
𝑼𝑡,±
𝑖+1∕2 =𝑼(𝑥𝑖+1∕2,𝑡)+(Δ𝑥𝑠).
Then, a semi-disc e e high-o de pa h conse a i e me hod is w i en as
𝑑
𝑑𝑡𝑼𝑖=− 1
Δ𝑥(𝑖+1∕2(𝑡)−𝑖−1∕2(𝑡) + 𝑺𝑖−𝑻𝑖),(6)
wi h 𝑖+1∕2(𝑡)=𝔽(𝑼𝑡,−
𝑖+1∕2,𝑼𝑡,+
𝑖+1∕2), whe e 𝔽(⋅,⋅)is a fi s o de consis en nume ical flux (𝔽(𝑼,𝑼)=𝑭(𝑼)) and
𝑺𝑖−𝑻𝑖≈
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2
𝑺(𝑷𝑡
𝑖(𝑥))𝜕𝑥𝑏 𝑑𝑥 −
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2
𝑻(𝑷𝑡
𝑖(𝑥)) 𝑑𝑥.
In p ac ice, we conside in his wo k a nume ical flux, which may be w i en as (see [39]):
𝔽(𝑼𝑡,−
𝑖+1∕2,𝑼𝑡,+
𝑖+1∕2)=1
2[𝑭(𝑼𝑡,+
𝑖+1∕2)+𝑭(𝑼𝑡,−
𝑖+1∕2)−(𝛽0(𝑼𝑡,+
𝑖+1∕2 −𝑼𝑡,−
𝑖+1∕2)+𝛽1(𝑭(𝑼𝑡,+
𝑖+1∕2)−𝑭(𝑼𝑡,−
𝑖+1∕2)))].(7)
No ice ha he defini ion o he nume ical scheme (6)assumes a con inuous bo om (o a con inuous econs uc ion o he bo om).
O he wise, a con ibu ion o he jumps ela ed o he discon inui ies o he opog aphy a he in e aces should be aken in o accoun .
3.1. Well-balanced scheme
We now adap scheme (6) o p ese e he lake-a - es equilib ia. To do so, le us ecall fi s he defini ions in oduced in [40]: we
shall dis inguish be ween exac ly well-balanced o well-balanced schemes. The o me p ese es he cell-a e ages o he exac s eady
s a e while he la e p ese es he cell-a e ages o a disc e e app oxima ion o he s a iona y solu ions. Mo e p ecisely,
Defini ion 1. Le 𝑼𝑒(𝑥)be any s a iona y solu ion o (4)and he co esponding cell a e ages {𝑼𝑒
𝑖}𝑖∈ob ained using (5)(o any
gi en high o de quad a u e o mula). Then, a nume ical me hod w i en in he o m (6)is said o be:
•exac ly well-balanced o 𝑼𝑒(𝑥)i he cell-a e ages {𝑼𝑒
𝑖}𝑖∈a e an equilib ium o he ODE sys em (6).
•well-balanced wi h o de 𝑟≥𝑠i , o e e y Δ𝑥, he e exis s a disc e e app oxima ion
𝑼𝑒
Δ𝑥,𝑖 =𝑼𝑒
𝑖+(Δ𝑥𝑟)
such ha {𝑼𝑒
Δ𝑥,𝑖}𝑖∈is an equilib ium o he ODE sys em (6).
Applied Nume ical Ma hema ics 215 (2025) 138–156
142
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
In o de o ge well-balanced schemes, ei he exac ly o no , we shall desc ibe a gene al algo i hm based on a econs uc ion
p ocedu e ollowing he ideas p oposed in [36]. Le us s a by assuming ha o any ime 𝑡, gi en he sequence {𝑼𝑖(𝑡)}, we selec
o e e y con ol olume 𝑉𝑖a con inuous s a iona y solu ion 𝑼𝑡,𝑒
𝑖(𝑥)o sys em (4)sa is ying he conse a ion p ope y
1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2
𝑼𝑡,𝑒
𝑖(𝑥) 𝑑𝑥 =𝑼𝑖(𝑡).
In pa icula , i holds
𝜕𝑥(𝑭(𝑼𝑡,𝑒
𝑖))+𝑺(𝑼𝑡,𝑒
𝑖)𝜕𝑥𝑏=𝑻(𝑼𝑡,𝑒
𝑖).(8)
Rema k ha al hough s a iona y solu ions do no depend on ime, he selec ed s a iona y solu ion a each ime may no be
necessa ily he same, hence he no a ion 𝑼𝑡,𝑒(𝑥).
Rema k 1. No ice ha one o he main difficul ies is o gi e an explici exp ession o he e m 𝑻(𝑼𝑡,𝑒
𝑖). Since he e we will ocus on
lake-a - es s eady s a es (3), he ope a o 𝑻is mul i-e alua ed and a p io i we only know ha
|||𝜏𝑥𝑧|𝑏∕𝜌|||≤𝑔ℎ𝜇𝑠. Ne e heless, since
𝑼𝑡,𝑒
𝑖should sa is y (8), his will ac ually fix he defini ion o 𝑻(𝑼𝑡,𝑒
𝑖). Mo e explici ly, he ic ion should balance (8)and he e o e
we ha e:
𝑻(𝑼𝑡,𝑒
𝑖(𝑥))=(0,𝑔ℎ𝑡,𝑒
𝑖(𝑥)𝜕𝑥𝜂𝑡,𝑒
𝑖(𝑥))𝑇.
Then, any semi-disc e e scheme (6)can be ew i en as
𝑑
𝑑𝑡𝑼𝑖=− 1
Δ𝑥(𝑖+1∕2(𝑡)−𝑭(𝑼𝑡,𝑒
𝑖(𝑥𝑖+1∕2)) − 𝑖−1∕2(𝑡)+𝑭(𝑼𝑡,𝑒
𝑖(𝑥𝑖−1∕2)))
−1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2 (𝑺(𝑷𝑡
𝑖(𝑥)) − 𝑺(𝑼𝑡,𝑒
𝑖(𝑥)))𝜕𝑥𝑏 𝑑𝑥 +
1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2 (𝑻(𝑷𝑡
𝑖(𝑥)) − 𝑻(𝑼𝑡,𝑒
𝑖(𝑥)))𝑑𝑥,
whe e we ha e se
𝑺𝑖−𝑻𝑖=−𝑭(𝑼𝑡,𝑒
𝑖(𝑥𝑖+1∕2))+𝑭(𝑼𝑡,𝑒
𝑖(𝑥𝑖−1∕2))+ 1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2 (𝑺(𝑷𝑡
𝑖(𝑥)) − 𝑺(𝑼𝑡,𝑒
𝑖(𝑥)))𝜕𝑥𝑏 𝑑𝑥 −
1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2 (𝑻(𝑷𝑡
𝑖(𝑥)) − 𝑻(𝑼𝑡,𝑒
𝑖(𝑥)))𝑑𝑥.
In p ac ice, he in eg als o he sou ce e ms will be compu ed by means o a quad a u e o mula wi h o de 𝑟≥𝑠, wi h 𝑠is he o de
o he econs uc ion ope a o :
𝑑
𝑑𝑡𝑼𝑖=− 1
Δ𝑥(𝑖+1∕2(𝑡)−𝑭(𝑼𝑡,𝑒
𝑖(𝑥𝑖+1∕2)) − 𝑖−1∕2(𝑡)+𝑭(𝑼𝑡,𝑒
𝑖(𝑥𝑖−1∕2)))
−
𝑘
∑
𝑗=1
𝜔𝑖
𝑗(𝑺(𝑷𝑡
𝑖(𝑥𝑖
𝑗)) − 𝑺(𝑼𝑡,𝑒
𝑖(𝑥𝑖
𝑗)))(𝜕𝑥𝑏)𝑥=𝑥𝑖
𝑗+
𝑘
∑
𝑗=1
𝜔𝑖
𝑗(𝑻(𝑷𝑡
𝑖(𝑥𝑖
𝑗)) − 𝑻(𝑼𝑡,𝑒
𝑖(𝑥𝑖
𝑗))),
(9)
whe e 𝑥𝑖
𝑗and 𝜔𝑖
𝑗, 𝑗=1,…,𝑘a e he quad a u e nodes in 𝑉𝑖and he co esponding weigh s, espec i ely, o he conside ed quad a u e
o mula.
Now, i is easy o p o e he ollowing esul :
Theo em 1. Le 𝑼𝑒(𝑥)be a con inuous s a iona y solu ion o (4), and he co esponding se o cell a e ages {𝑼𝑒
𝑖}𝑖∈defined as in (5).
Assume ha he co esponding selec ed in-cell s eady s a es e i y ei he one o he p ope ies
(i) 𝑼𝑡,𝑒
𝑖(𝑥)=𝑼𝑒(𝑥), o all 𝑥∈𝑉𝑖and o all 𝑖∈,
(ii) 1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2
𝑼𝑡,𝑒
𝑖(𝑥)𝑑𝑥 =𝑼𝑒
𝑖+(Δ𝑥𝑟), o all 𝑥∈𝑉𝑖and o all 𝑖∈and he global solu ion defined piece-wisely by 𝑼𝑡,𝑒
𝑖(𝑥)in 𝑉𝑖is
con inuous.
Suppose ha he econs uc ion ope a o 𝑷𝑡
𝑖(𝑥)is well-balanced o {𝑼𝑡,𝑒
𝑖}𝑖∈in he sense
𝑷𝑡
𝑖(𝑥;{𝑼𝑒
𝑗}𝑗∈𝑖)=𝑼𝑡,𝑒
𝑖(𝑥), ∀𝑥∈𝑉𝑖, ∀𝑖∈.
Then, he nume ical scheme (9)is ei he exac ly well-balanced (i (i) is sa isfied) o well-balanced wi h o de 𝑟(i (ii) is sa isfied) o
𝑼𝑒(𝑥).

Applied Nume ical Ma hema ics 215 (2025) 138–156
143
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
P oo . Thanks o he con inui y o he solu ion 𝑼𝑡,𝑒(𝑥), we au oma ically ge 𝑖+1∕2(𝑡)=𝑭(𝑼𝑡,𝑒
𝑖(𝑥𝑖+1∕2)) and all he e ms on he
igh hand side o (9)cancel, so he esul ollows. □
3.2. Con inuous s a iona y solu ion and econs uc ion ope a o
In his sec ion we will desc ibe an algo i hm ha allows o define he in-cell s eady s a es 𝑼𝑡,𝑒
𝑖(𝑥) o lake-a - es equilib ium (2)
(ei he exac ly o app oxima ed). We will define as well well-balanced second-o de econs uc ion ope a o 𝑷𝑡
𝑖so ha he esul ing
scheme (9)is well-balanced wi h o de 2. Fo he sake o simplici y, we shall desc ibe fi s in Subsec ion 3.2.1 he echnique o he
second-o de scheme, al hough i can be ex ended o a high-o de e sion as i is shown in Subsec ion 3.2.2.
Gi en a collec ion o cell-a e ages {𝑼𝑖(𝑡)}, we shall desc ibe fi s he cons uc ion o he in-cell s eady s a es {𝑼𝑡,𝑒
𝑖(𝑥)} and hen
he econs uc ion ope a o 𝑷𝑡
𝑖(𝑥). In wha ollows, we will no w i e explici ly he dependence on ime in o de o make he no a ion
less cumbe some.
3.2.1. Second o de me hod
We s a by desc ibing he s a iona y solu ion and he econs uc ion p ocedu e in he case o a second-o de well-balanced scheme.
In-cell s eady s a es. Assume a gi en unc ion 𝑼(𝑥)and he co esponding a e age alues {𝑼𝑖}compu ed as in (5)ei he exac ly o
using a high o de quad a u e o mula. We hen define he in-cell s eady s a e unc ions 𝑼𝑒
𝑖(𝑥)=(𝜂𝑒
𝑖(𝑥),0) as ollows:
1. Conside fi s

𝑷𝑖(𝑥)=(𝑝ℎ
𝑖,𝑝𝑞
𝑖)𝑇wi h 𝑝ℎ
𝑖,𝑝𝑞
𝑖∈ℙ1[𝑥] he usual MUSCL econs uc ion ope a o wi h he minmod slope limi e .
As i is usual in shallow wa e amewo k, we use he econs uc ion o he ee su ace alues 𝜂𝑖=ℎ𝑖+𝑏𝑖 o define he heigh
econs uc ion, ha is, we define:
𝑝𝑤
𝑖(𝑥)=𝑤𝑖+(𝛿𝑤)𝑖
Δ𝑥 (𝑥−𝑥𝑖), wi h (𝛿𝑤)𝑖=minmod(𝑤𝑖−𝑤𝑖−1,𝑤𝑖+1 −𝑤𝑖)
o 𝑤∈{𝜂,𝑞}, and 𝑝ℎ
𝑖(𝑥)=𝑝𝜂
𝑖(𝑥)−𝑏(𝑥). Rema k ha s ic ly speaking 𝑝ℎ
𝑖(𝑥)is no a polynomial in gene al, bu his does no
affec o wha is said a e wa ds, since e e y hing is defined in e ms o (𝛿𝜂)𝑖.
2. Define he ope a o 𝑟𝜂
𝑖∈ℙ1[𝑥]by
𝑟𝜂
𝑖(𝑥)=𝜂𝑖+𝛼𝑖(𝑝𝜂
𝑖(𝑥)−𝜂𝑖)wi h 𝛼𝑖=⎧
⎪
⎨
⎪
⎩
1i |(𝛿𝜂)𝑖|≤Δ𝑥𝜇𝑠,
𝜇𝑠Δ𝑥
(𝛿𝜂)𝑖
o he wise.
No ice ha 𝑟𝜂
𝑖(𝑥)is a linea piecewise unc ion wi h slope less o equal o he angen o he angle o epose o he ma e ial 𝜇𝑠.
Mo eo e , 𝑟𝜂
𝑖sa isfies he conse a ion p ope y
1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2
𝑟𝜂
𝑖(𝑥)𝑑𝑥 =𝜂𝑖.
Rema k ha in he case ha

𝑷𝑖(𝑥)defines a lake-a - es s a iona y s a e in he cell 𝑉𝑖, hen 𝑟𝜂
𝑖(𝑥)=𝑝𝜂
𝑖(𝑥). Howe e , one may find
in gene al discon inui ies a he in e aces 𝑥𝑖+1∕2, i.e., he piecewise unc ion 𝑟𝜂(𝑥)gi en by 𝑟𝜂
𝑖(𝑥)in 𝑉𝑖is no a global s a iona y
solu ion.
3. Conside now a each in e ace 𝑥𝑖+1∕2 he le (𝜂−
𝑖+1∕2) and igh (𝜂+
𝑖+1∕2) limi s, as well as he mean in e -cell alue (𝜂𝑖+1∕2) gi en
by
𝜂−
𝑖+1∕2 =𝑟𝜂
𝑖(𝑥𝑖+1∕2), 𝜂+
𝑖+1∕2 =𝑟𝜂
𝑖+1(𝑥𝑖+1∕2), 𝜂𝑖+1∕2 =1
2(𝜂−
𝑖+1∕2 +𝜂+
𝑖+1∕2),(10)
and conside 𝑟𝜂
𝑖(𝑥)∈ℙ2[𝑥] he quad a ic polynomial e i ying
𝑟𝜂
𝑖(𝑥𝑖)=𝜂𝑖, and 𝑟𝜂
𝑖(𝑥𝑖±1∕2)=𝜂𝑖±1∕2.
Rema k ha his quad a ic polynomial sa isfies he conse a ion p ope y up o second o de :
1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2
𝑟𝜂
𝑖(𝑥)𝑑𝑥 =𝜂𝑖+(Δ𝑥2).
Then, conside he maximum slope in he olume 𝑀𝑖=max
𝑥∈𝑉𝑖|||(𝑟𝜂
𝑖)′(𝑥)|||and se
𝜂𝑒
𝑖(𝑥)={𝑟𝜂
𝑖(𝑥)i 𝑀𝑖≤𝜇𝑠,
𝑟𝜂
𝑖(𝑥)o he wise.
Applied Nume ical Ma hema ics 215 (2025) 138–156
144
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
Fig. 1. In-cell s a iona y s eady s a es o he unc ion −𝜇𝑠𝑥2∕2 + 1. Le figu e shows he linea limi ed slope unc ion 𝑟𝑒
𝑖and igh figu e shows he final s eady s a e
𝜂𝑒
𝑖, including a zoom.
Fig. 2. In-cell s a iona y s eady s a es o he unc ion 𝑥4. Le figu e shows he linea limi ed slope unc ion 𝑟𝑒
𝑖and igh figu e shows he final s eady s a e 𝜂𝑒
𝑖. In blue
cells whe e 𝜂𝑒
𝑖=𝑟𝑒
𝑖and in ed whe e 𝜂𝑒
𝑖=𝑟𝑒
𝑖. (Fo in e p e a ion o he colo s in he figu e(s), he eade is e e ed o he web e sion o his a icle.)
Suppose ha he in e nal slopes wi hin he cells, 𝑀𝑖, a e unde he limi ing angen o he epose angle 𝜇𝑠. Then, he piece-wise
unc ion gi en by 𝜂𝑒
𝑖(𝑥)in 𝑉𝑖is a s a iona y solu ion es ic ed o he in e io o each cell 𝑉𝑖and con inuous a he in e aces,
he e o e i is a global con inuous s eady s a e.
To illus a e he cons uc ion o he in-cell s eady s a e, we show in Fig. 1 he case o he ini ial unc ion 𝜂(𝑥, 0) = −𝜇𝑠𝑥2∕2 + 1 in
he in e al [−1,1]. In his case, since 𝑀𝑖≤1, 𝜂𝑒
𝑖co esponds o he pa abola 𝑟𝜂
𝑖. Rema k ha in his case, he global unc ion defined
piece-wise as 𝜂𝑒
𝑖inside each cell coincides (up o second o de ) wi h he ini ial unc ion a disc e e le el, ha is, hey cell-a e aged
alues coincide. In a simila way, Fig. 2shows he case o ini ial unc ion 𝜂(𝑥, 0) = 𝑥4in he in e al [0,1.2]. In his case, he slope
on he igh pa o he unc ion is g ea e han 𝜇𝑠, which means ha we selec he limi ed slope linea unc ion 𝑟𝑖, while on he le
pa 𝜂𝑒
𝑖coincides wi h he pa abola 𝑟𝜂
𝑖.
Recons uc ion ope a o . Using he same no a ion defined o he in-cell s eady s a es, he econs uc ion ope a o 𝑷𝑖(𝑥)inside he
olume 𝑉𝑖will be se as:
𝑷𝑖(𝑥)={(𝜂𝑒
𝑖(𝑥)−𝑏(𝑥), 𝑝𝑞
𝑖(𝑥))𝑇,i 𝛼𝑖=1,
(𝑝𝜂
𝑖(𝑥)−𝑏(𝑥), 𝑝𝑞
𝑖(𝑥))𝑇,o he wise.
No ice ha he s ep 3 in p e ious algo i hm cha ac e ize he solu ions we a e going o p ese e: hose s a iona y solu ion whose
disc e e ep esen a ion by a quad a ic polynomial 𝑟𝜂(𝑥)a e s a iona y (||(𝑟𝜂)′(𝑥)||<𝜇
𝑠). In pa icula , i holds o linea piecewise
solu ions whose slopes a e lowe han 𝜇𝑠in each olume, p o ided ha we a e able o ha e a exac app oxima ion o (𝛿𝜂)𝑖.
Applied Nume ical Ma hema ics 215 (2025) 138–156
145
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
Theo em 2. Le 𝑼(𝑥)=(ℎ(𝑥),0)𝑇such ha 𝜂(𝑥)=ℎ(𝑥)+𝑏(𝑥)is sufficien ly smoo h and sa isfies |𝜂′(𝑥)|≤𝑀<𝜇
𝑠. Then he scheme (9)
is well-balanced o 𝑼(𝑥).
P oo . We will make use o Theo em 1. Fi s , ema k ha o Δ𝑥sufficien ly small, we may assume
|(𝑝𝜂
𝑖)′(𝑥)|=|(𝛿𝜂)𝑖|
Δ𝑥
<𝜇
𝑠.
The e o e we ha e 𝑟𝜂
𝑖(𝑥)=𝑝𝜂
𝑖(𝑥). Mo eo e , assuming he unc ion sufficien ly smoo h, we will ha e 𝜂𝑖+1∕2 =𝜂(𝑥𝑖+1∕2)+(Δ𝑥2).
The quad a ic polynomial 𝑟𝜂
𝑖(𝑥)can be w i en as
𝑟𝜂
𝑖(𝑥)=𝜂𝑖−1∕2 +2
Δ𝑥(𝜂𝑖−𝜂𝑖−1∕2)(𝑥−𝑥𝑖−1∕2)+2
Δ𝑥2(𝜂𝑖+1∕2 −2𝜂𝑖+𝜂𝑖−1∕2)(𝑥−𝑥𝑖−1∕2)(𝑥−𝑥𝑖).
Since (𝑟𝜂
𝑖)′(𝑥)is a linea polynomial, he maximum 𝑀𝑖is eached a one o he in e aces o he olume 𝑉𝑖:
(𝑟𝜂
𝑖)′(𝑥𝑖+1∕2)= 1
Δ𝑥(𝜂𝑖−1∕2 −𝜂𝑖+3(𝜂𝑖+1∕2 −𝜂𝑖))=𝜂′(𝑥𝑖)+(Δ𝑥)
(𝑟𝜂
𝑖)′(𝑥𝑖−1∕2)= 1
Δ𝑥(3(𝜂𝑖−𝜂𝑖−1∕2)+(𝜂𝑖−𝜂𝑖+1∕2))=𝜂′(𝑥𝑖)+(Δ𝑥).
Then o Δ𝑥sufficien ly small we ge 𝑀𝑖<𝜇
𝑠and we ha e ha he in-cell s eady s a e co esponds o he pa abola, 𝜂𝑒
𝑖(𝑥)=𝑟𝜂
𝑖(𝑥),
which is an app oxima ion o 𝜂(𝑥). Mo eo e , 𝑷𝑖(𝑥)=(𝜂𝑒
𝑖(𝑥)−𝑏(𝑥),0)𝑇and he esul ollows. □
Rema k 2. Al hough we ha e p o ed ha he scheme (9)is in gene al well-balanced, we may ha e he exac ly well-balanced p ope y
in some si ua ions. Fo ins ance, i 𝑼(𝑥)=(ℎ(𝑥),0)𝑇is such ha ℎ(𝑥)+𝑏(𝑥)co esponds o a linea polynomial wi h slope smalle
han 𝜇𝑠, hen he scheme is exac ly well-balanced. Indeed, in such cases he slope app oxima ions (𝛿𝜂)𝑖co espond exac ly o 𝜂′(𝑥)
and i is easy o check ha 𝜂𝑖(𝑥)=𝜂(𝑥)=𝑝𝜂
𝑖(𝑥).
O special in e es is he case o wo connec ed linea unc ions (o a con inuous polygonal in gene al) wi h slopes smalle han
𝜇𝑠. This is s udied in Appendix A.
Ano he pa icula case is when he ini ial condi ion p oduces a sequence o cell a e ages {𝜂𝑖}𝑖∈such ha he econs uc ion
p ocedu e esul s in 𝑟𝜂
𝑖inside he cells and hey connec con inuously a he in e ace. In ha case, a disc e e le el, bo h he ini ial
condi ion and he econs uc ed unc ion coincide. This ini ial condi ion is he e o e p ese ed.
3.2.2. Gene aliza ion o highe o de
The algo i hm in oduced p e iously, can be adap ed o ob ain high o de well-balanced schemes as ollows:
In-cell s eady s a es.
1. We conside now

𝑷𝑖(𝑥)=(𝑝𝜂
𝑖−𝑏(𝑥),𝑝𝑞
𝑖)𝑇wi h 𝑝ℎ
𝑖,𝑝𝑞
𝑖∈ℙ𝑠[𝑥]a high o de polynomial econs uc ion ope a o (WENO, CWENO,
e c.). Again, we a he use he econs uc ions on he ee su ace 𝑝𝜂
𝑖.
2. Define he ope a o 𝑟𝜂
𝑖∈ℙ𝑠[𝑥]by
𝑟𝜂
𝑖(𝑥)=𝜂𝑖+𝛼𝑖(𝑝𝜂
𝑖(𝑥)−𝜂𝑖)wi h 𝛼𝑖={1i |(𝑝𝜂
𝑖)′(𝑥)|≤𝜇𝑠,
𝜇𝑠∕𝑀o he wise,
wi h 𝑀=max
𝑥∈𝑉𝑖|(𝑝𝜂
𝑖)′(𝑥)|. No ice ha , since 𝑝𝜂
𝑖(𝑥)is assumed o sa is y he conse a ion p ope y, so does 𝑟𝜂
𝑖(𝑥):
1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2
𝑟𝜂
𝑖(𝑥)𝑑𝑥 =𝜂𝑖.
3. Conside now a each in e ace 𝜂±
𝑖+1∕2,𝜂𝑖+1∕2 as in (10), and a Gauss quad a u e o mula o deg ee 𝑟≥𝑠in 𝑉𝑖wi h nodes 𝑥𝑖
𝑗∈
(𝑥𝑖−1∕2,𝑥𝑖+1∕2)and weigh s 𝜔𝑖
𝑗, 𝑗=1,…,𝑘.
Conside 𝑟𝜂
𝑖(𝑥)∈ℙ𝑘+1[𝑥] he polynomial e i ying
𝑟𝜂
𝑖(𝑥𝑖
𝑗)=𝑝𝜂
𝑖(𝑥𝑖
𝑗), 𝑗=1,…,𝑘, and 𝑟𝜂
𝑖(𝑥𝑖±1∕2)=𝜂𝑖±1∕2.
Rema k ha his polynomial sa isfies he conse a ion p ope y up o o de 𝑟≥𝑠:
1
Δ𝑥
𝑥𝑖+1∕2
∫
𝑥𝑖−1∕2
𝑟𝜂
𝑖(𝑥)𝑑𝑥 =
𝑘
∑
𝑗=1
𝜔𝑖
𝑗𝑟𝜂
𝑖(𝑥𝑖
𝑗)+(Δ𝑥𝑟)=
𝑘
∑
𝑗=1
𝜔𝑖
𝑗𝑝𝜂
𝑖(𝑥𝑖
𝑗)+(Δ𝑥𝑟)=𝜂𝑖+(Δ𝑥𝑟).
Applied Nume ical Ma hema ics 215 (2025) 138–156
146
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
Then, conside he maximum slope in he olume 𝑀𝑖=max
𝑥∈𝑉𝑖|||(𝑟𝜂
𝑖)′(𝑥)|||and se
𝜂𝑒
𝑖(𝑥)={𝑟𝜂
𝑖(𝑥)i 𝑀𝑖≤𝜇𝑠,
𝑟𝜂
𝑖(𝑥)o he wise.
Recons uc ion ope a o . Simila ly as in he second-o de case, he econs uc ion ope a o 𝑷𝑖(𝑥)inside he olume 𝑉𝑖will be se as:
𝑷𝑖(𝑥)={(𝜂𝑒
𝑖(𝑥)−𝑏(𝑥), 𝑝𝑞
𝑖(𝑥))𝑇,i 𝛼𝑖=1,
(𝑝𝜂
𝑖(𝑥)−𝑏(𝑥), 𝑝𝑞
𝑖(𝑥))𝑇,o he wise.
Then he well-balance p ope y o he econs uc ion ope a o ollows analogously.
Theo em 3. Le 𝑼(𝑥)=(ℎ(𝑥),0)𝑇such ha 𝜂(𝑥)=ℎ(𝑥)+𝑏(𝑥)is sufficien ly smoo h and sa isfies |𝜂′(𝑥)|≤𝑀<𝜇
𝑠. Assume a econs uc ion
ope a o 𝑝𝜂
𝑖(𝑥)such ha
𝜂(𝑥)=𝑝𝜂
𝑖(𝑥)+(Δ𝑥𝑝), 𝜂′(𝑥)=(𝑝𝜂
𝑖)′(𝑥)+(Δ𝑥𝑚), ∀𝑥∈𝑉𝑖
o some in ege s 1≤𝑚<𝑝.
Then he scheme (9)is well-balanced o 𝑼(𝑥).
P oo . We will make use o Theo em 1. Fi s , ema k ha o Δ𝑥sufficien ly small, we may assume
|(𝑝𝜂
𝑖)′(𝑥)|≤|𝜂′(𝑥)|+(Δ𝑥𝑚)<𝜇
𝑠.
The e o e, we ha e 𝑟𝜂
𝑖(𝑥)=𝑝𝜂
𝑖(𝑥). Mo eo e , assuming he unc ion sufficien ly smoo h, we will ha e 𝜂𝑖+1∕2 =𝜂(𝑥𝑖+1∕2)+(Δ𝑥𝑝).
Now, conside 𝑟𝜂
𝑖(𝑥)∈ℙ𝑘+1[𝑥] he polynomial e i ying
𝑟𝜂
𝑖(𝑥𝑖
𝑗)=𝑝𝜂
𝑖(𝑥𝑖
𝑗), 𝑗=1,…,𝑘, and 𝑟𝜂
𝑖(𝑥𝑖±1∕2)=𝜂𝑖±1∕2,
and we shall assume 𝑥𝑖−1∕2 <𝑥𝑖
1<…<𝑥𝑖
𝑘<𝑥
𝑖+1∕2. Le us use he no a ion 𝑥𝑖
0=𝑥𝑖−1∕2, 𝑥𝑖
𝑘+1 =𝑥𝑖+1∕2.
Now, le 𝑥∈𝑉𝑖and conside 0≤𝑗≤𝑘such ha 𝑥∈[𝑥𝑖
𝑗,𝑥𝑖
𝑗+1]. Then, we ha e
(𝑟𝜂
𝑖)′(𝑥)=
𝑟𝜂
𝑖(𝑥𝑖
𝑗+1)−𝑟𝜂
𝑖(𝑥𝑖
𝑗)
𝑥𝑖
𝑗+1 −𝑥𝑖
𝑗
+(Δ𝑥)=
𝜂(𝑥𝑖
𝑗+1)−𝜂(𝑥𝑖
𝑗)+(Δ𝑥𝑝)
𝑥𝑖
𝑗+1 −𝑥𝑖
𝑗
+(Δ𝑥)=𝜂′(𝑥)+(Δ𝑥).
The e o e, o Δ𝑥sufficien ly small we ge 𝑀𝑖<𝜇
𝑠and we ha e ha he in-cell s eady s a e co esponds o he polynomial, 𝜂𝑒
𝑖(𝑥)=
𝑟𝜂
𝑖(𝑥), which is an app oxima ion o 𝜂(𝑥). Mo eo e , 𝑷𝑖(𝑥)=(𝜂𝑒
𝑖(𝑥)−𝑏(𝑥),0)𝑇and he esul ollows. □
3.3. Time disc e iza ion
We de ail in his subsec ion he ime disc e iza ion o he semi-disc e e sys em (9). Le us ecall ha he main difficul y he e is
he defini ion o he ic ion e m. I is usual, when dealing wi h ic ion e ms (Coulomb, Manning, Da cy,...), o use a semi-implici
app oach. In pa icula , o he case o single-laye models o g anula flows, he echnique in oduced in [21] can be used. The
basis idea is o assume ha he ic ion o ce is opposed o he mo ion o he flow, bu i canno change he sign o he eloci y ec o
(di ec ion o mo emen ). Mo eo e , i he magni ude o he ic ion o ce is g ea e han he es o o ces, hen he flow mus s op.
Fo he sake o cla i y, le us summa ize his echnique o a e y simple case. Conside he ODE
𝑑𝑢
𝑑𝑡
=−𝑔𝜇ℎ sgn(𝑢),(11)
whe e sgn(⋅)is he sign unc ion and 𝜇>0a ic ion coefficien . Then, applying he semi-implici me hod wi h a linea iza ion o he
p essu e e m, we ge
𝑢𝑛+1 =𝑢𝑛−Δ𝑡𝑔𝜇ℎ𝑛sgn(𝑢𝑛+1),
whe e he supe sc ip s 𝑛,𝑛 +1deno e he solu ions a imes 𝑡𝑛and 𝑡𝑛+1 =𝑡𝑛+Δ𝑡. Since he ic ion e m canno change he sign o
he eloci y, hen sgn(𝑢𝑛+1)=sgn(𝑢𝑛)is assumed, and he ollowing co ec ion s ep on he upda e o 𝑢is pe o med:
𝑢𝑛+1 ={𝑢𝑛−Δ𝑡𝑔𝜇ℎ𝑛sgn(𝑢𝑛)i |𝑢𝑛|>Δ𝑡𝑔ℎ𝑛𝜇
0o he wise. (12)
Thus, looking a p e ious equa ion, we obse e ha in p ac ice (11)is explici ly disc e ized, excep o he co ec ion s ep, which
se s 𝑢𝑛+1 =0when needed. This explici ea men o he eloci y is an ad an age when compa ed o o he echniques such as a
egula iza ion o he sign unc ion. In pa icula , no ice ha null eloci y is ob ained exac ly when he ic ion su passes a c i ical
le el, ep oducing p ope ly he case o ma e ial a es and s opping c i e ion.
Applied Nume ical Ma hema ics 215 (2025) 138–156
153
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
Fig. A.9. Ske ch o wo connec ed linea s eady s a es wi h slopes 𝛽𝐿, 𝛽𝑅.
5. Conclusions
We ha e in oduced a well-balanced econs uc ion p ocedu e o shallow wa e g anula flows model wi h Coulomb- ype ic-
ion e m, which allows o p ese e lake-a - es s eady s a es. This p ocedu e, oge he wi h an app op ia e ime disc e iza ion, is
immedia ely gene alized o an a bi a y o de o accu acy. The esul ing scheme is well-balanced in gene al o such s eady s a es,
al hough i is exac ly well-balanced o some pa icula cases. Conc e ely, hose s eady s a es whose disc e e ep esen a ion, defined
in he o m o cell a e aged alues, co esponds o a global con inuous s a iona y solu ion desc ibed inside he cells by he s a iona y
econs uc ion p ocedu e desc ibed in Subsec ion 3.2 a e exac ly p ese ed. Rema k ha he wo k p esen ed he e is easily ex ended
o he case o il ed coo dina es and weakly non-hyd os a ic amewo k (see [14]).
All his has been shown in he nume ical es s sec ion. In pa icula , i is wo h no icing ha when he scheme is no exac ly
well-balanced o some ini ial s eady s a e, he ma e ial sligh ly flows o a e y sho ime and hen an app oxima ion o he ini ial
condi ion is ob ained, which is hen p ese ed. Tha is, he scheme is well-balanced. Mo eo e , i has been shown ha he expec ed
o de o accu acy is ob ained o he second and hi d-o de schemes. Le us ema k ha , up o ou knowledge, his is he fi s ime
ha a high-o de well-balanced fini e olume me hod has been p oposed o such g anula flow model. Fu u e wo ks could s udy
he ex ension o his echnique o sys ems including he e ical disc e iza ion o he g anula sys em, whe e he e ical effec s a e
essen ial o ob ain ealis ic simula ion in complex iscoplas ic flows.
CRediT au ho ship con ibu ion s a emen
M.J. Cas o Díaz: W i ing – e iew & edi ing, So wa e, Me hodology, Funding acquisi ion, Concep ualiza ion. C. Escalan e:
W i ing – e iew & edi ing, W i ing – o iginal d a , So wa e, Me hodology, Concep ualiza ion. J. Ga es-Díaz: W i ing – e iew
& edi ing, W i ing – o iginal d a , So wa e, Me hodology, Concep ualiza ion. T. Mo ales de Luna: W i ing – e iew & edi ing,
So wa e, Me hodology, Funding acquisi ion, Concep ualiza ion.
Acknowledgemen s
This wo k is suppo ed by p ojec s PID2022-137637NB-C21 and PID2022-137637NB-C22 unded by MICIU/AEI/10.13039/
501100011033/ and FEDER, UE, and p ojec PDC2022-133663-C21 unded by MICIU/AEI/10.13039/501100011033/ and Eu o-
pean Union Nex Gene a ionEU/PRTR. The wo k has also been pa ially suppo ed by p ojec PROYEXCEL_00525 unded by Jun a de
Andalucía.
Appendix A. The case o connec ed s a iona y linea s a es
In his Appendix we conside fi s he case o wo connec ed linea polynomials and hen he mo e gene al case o a polygonal,
as hey may be o pa icula in e es .
Le 𝛽𝐿,𝛽𝑅,𝛾 ∈ℝ, wi h |𝛽𝐿|≤𝜇𝑠, |𝛽𝑅|≤𝜇𝑠and conside
𝜂(𝑥)={𝛽𝐿(𝑥−𝑐)+𝛾, o 𝑥≤𝑐,
𝛽𝑅(𝑥−𝑐)+𝛾, o 𝑥>𝑐. (A.1)
We will conside a space disc e iza ion 𝑉𝑖=[𝑥𝑖−1∕2,𝑥𝑖+1∕2]such ha 𝑥𝑗+1∕2 =𝑐 o some index 𝑗(see Fig. A.9).
Mo eo e , we shall iden i y he cell a e ages wi h he alues a he cen e o he cells, since hey a e second-o de app oxima ions.
Rema k ha he cons uc ions desc ibed in Subsec ion 3.2 a e in a ian by e ical ansla ion in he sense ha i we add a cons an
𝛾 o he cell a e ages, hen he econs uc ions coincide wi h hose o he o iginal alues ansla ed by he cons an 𝛾. The e o e, we

Applied Nume ical Ma hema ics 215 (2025) 138–156
154
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
assume wi hou loss o gene ali y 𝛾=0. Mo eo e , we suppose |𝛽𝐿|>|𝛽𝑅|, o he wise we may conside a symme y wi h espec o
he e ical axis 𝑥=𝑥𝑗+1∕2 and e e y hing ollows analogously.
Rema k ha he econs uc ed slopes (𝛿𝜂)𝑖∕Δ𝑥co espond exac ly o 𝛽𝐿 o 𝑖≤𝑗−1and o 𝛽𝑅 o 𝑖≥𝑗+2. The only pa icula
cases a e hen he cells 𝑗and 𝑗+1, whe e we ge
(𝛿𝜂)𝑗
Δ𝑥
=minmod(𝛽𝐿+𝛽𝑅
2
,𝛽𝐿),
(𝛿𝜂)𝑗+1
Δ𝑥
=minmod(𝛽𝐿+𝛽𝑅
2
,𝛽𝑅).
A sys ema ic s udy o he diffe en possibili ies when |𝛽𝐿|>|𝛽𝑅|le us show ha
minmod(𝛽𝐿+𝛽𝑅
2
,𝛽𝐿)=𝛽𝐿+𝛽𝑅
2
,
minmod(𝛽𝐿+𝛽𝑅
2
,𝛽𝑅)={𝛽𝑅,i 𝛽𝐿𝛽𝑅≥0,
0,o he wise.
I is clea hen ha all slopes a e smalle han he c i ical alue 𝜇𝑠and 𝛼𝑖=1. Now, we easily check ha a in e ace 𝑥𝑗+1∕2 we ha e
𝜂−
𝑗+1∕2 =Δ𝑥
4
(𝛽𝑅−𝛽𝐿), 𝜂+
𝑗+1∕2 =⎧
⎪
⎨
⎪
⎩
0,i 𝛽𝐿𝛽𝑅≥0,
Δ𝑥
2
𝛽𝑅,o he wise, 𝜂𝑗+1∕2 =⎧
⎪
⎨
⎪
⎩
Δ𝑥
8
(𝛽𝑅−𝛽𝐿),i 𝛽𝐿𝛽𝑅≥0,
Δ𝑥
8
(3𝛽𝑅−𝛽𝐿),o he wise,
a he le in e ace, 𝑥𝑗−1∕2,
𝜂−
𝑗−1∕2 =−Δ𝑥𝛽𝐿, 𝜂+
𝑗−1∕2 =−Δ𝑥(3
4𝛽𝐿+1
4𝛽𝑅), 𝜂𝑗−1∕2 =−Δ𝑥(7
8𝛽𝐿+1
8𝛽𝑅),
and a he igh , 𝑥𝑗+3∕2,
𝜂−
𝑗+3∕2 =⎧
⎪
⎨
⎪
⎩
Δ𝑥𝛽𝑅,i 𝛽𝐿𝛽𝑅≥0,
Δ𝑥
2
𝛽𝑅,o he wise, 𝜂+
𝑗+3∕2 =Δ𝑥𝛽𝑅, 𝜂𝑗+3∕2 =⎧
⎪
⎨
⎪
⎩
Δ𝑥𝛽𝑅,i 𝛽𝐿𝛽𝑅≥0,
3Δ𝑥
4
𝛽𝑅,o he wise.
Fo any cell 𝑉𝑖wi h 𝑖≠𝑗−1,𝑗,𝑗+1,𝑗+2, he second-o de polynomial 𝑟𝑖(𝑥)coincides wi h he o iginal linea unc ion inside he
cell and we ha e 𝜂𝑖(𝑥)=𝑝𝜂
𝑖(𝑥)=𝜂(𝑥) o 𝑥∈𝑉𝑖. Now, in o de o conside he pa icula case 𝑘∈{𝑗−1,𝑗,𝑗+1,𝑗+2}, le us w i e
he gene al polynomial 𝑟𝑘(𝑥)∈ℙ2[2] in he o m:
𝑟𝑘(𝑥)=𝜂𝑘−1∕2 +2𝜂𝑘−𝜂𝑘−1∕2
Δ𝑥
(𝑥−𝑥𝑘−1∕2)+2𝜂𝑘+1∕2 −2𝜂𝑘+𝜂𝑘−1∕2
Δ𝑥2(𝑥−𝑥𝑘−1∕2)(𝑥−𝑥𝑘),
which gi es
𝑟′
𝑘(𝑥)= 𝜂𝑘+1∕2 −𝜂𝑘−1∕2
Δ𝑥
+4𝜂𝑘+1∕2 −2𝜂𝑘+𝜂𝑘−1∕2
Δ𝑥2(𝑥−𝑥𝑘).
Since 𝑟′
𝑘(𝑥)is linea , we ge
max
𝑥∈𝑉𝑘|𝑟′
𝑘(𝑥)|=max{|𝑟′
𝑘(𝑥𝑘−1∕2)|,|𝑟′
𝑘(𝑥𝑘+1∕2)|}= 1
Δ𝑥max {|𝜂𝑘−𝜂𝑘+1∕2 +3(𝜂𝑘−𝜂𝑘−1∕2)|,|3(𝜂𝑘+1∕2 −𝜂𝑘)+𝜂𝑘−1∕2 −𝜂𝑘|}.
A e some easy calcula ions we ha e:
•In he cell 𝑉𝑗−1 wo possibili ies a ise:
–I |11𝛽𝐿−3𝛽𝑅|≤8𝜇𝑠, hen 𝑀𝑗−1 ≤𝜇𝑠and we se 𝜂𝑒
𝑗−1(𝑥)=𝑟𝑗−1(𝑥).
–I |11𝛽𝐿−3𝛽𝑅|>8𝜇𝑠, hen 𝑀𝑗−1 >𝜇
𝑠and we se 𝜂𝑒
𝑗−1(𝑥)=𝑟𝑗−1(𝑥)=𝑝𝜂
𝑗−1(𝑥).
•In he cell 𝑉𝑗we ha e 𝑀𝑗≤𝜇𝑠and 𝜂𝑒
𝑗(𝑥)=𝑟𝑗(𝑥). Mo eo e , i 𝛽𝐿𝛽𝑅≥0, hen 𝑟𝑗(𝑥)is linea .
•In he cell 𝑉𝑗+1 we ha e 𝑀𝑗+1 ≤𝜇𝑠and 𝜂𝑒
𝑗+1(𝑥)=𝑟𝑗+1(𝑥).
•In he cell 𝑉𝑗+2 wo possibili ies a ise:
–I 𝛽𝐿𝛽𝑅≥0, hen 𝑟𝑗+2(𝑥)is linea and coincides wi h he o iginal unc ion inside he cell: 𝜂𝑒
𝑗+2(𝑥)=𝑟𝑗+2(𝑥)=𝜂(𝑥).
–I 𝛽𝐿𝛽𝑅<0, hen we only ge 𝑀𝑗+2 <𝜇
𝑠in he case |𝛽𝑅|≤4
7𝜇𝑠.
𝜂𝑒
𝑗+2(𝑥)=⎧
⎪
⎨
⎪
⎩
𝑟𝑗+2(𝑥),i 𝛽𝐿𝛽𝑅<0and |𝛽𝑅|≤4
7𝜇𝑠,
𝑟𝑗+2(𝑥)=𝑝𝜂
𝑗+2(𝑥),i 𝛽𝐿𝛽𝑅<0and |𝛽𝑅|>4
7𝜇𝑠.
Applied Nume ical Ma hema ics 215 (2025) 138–156
155
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
The e o e, we conclude ha when 𝛽𝐿and 𝛽𝑅a e such ha
{|11𝛽𝐿−3𝛽𝑅|≤8𝜇𝑠}and {𝛽𝐿𝛽𝑅≥0o |𝛽𝑅|≤4
7𝜇𝑠},(A.2)
he econs uc ion p ocedu e o 𝜂(𝑥)is such ha i is (exac ly) p ese ed by he scheme.
I he p e ious condi ions a e no sa isfied, hen we po en ially ge wo cells 𝑉𝑗−1 and 𝑉𝑗+2 whe e he unc ion is no necessa ily
p ese ed.
We ecall ha we ha e assumed he case
||𝛽𝐿||>||𝛽𝑅||. O he wise, by symme y easons one ob ains
{|11𝛽𝑅−3𝛽𝐿|≤8𝜇𝑠}and {𝛽𝐿𝛽𝑅≥0o |𝛽𝐿|≤4
7𝜇𝑠}.
Rema k ha his can be easily gene alized. Conside now he gene al case, whe e he ini ial condi ion co esponds o a con inuous
piece-wise linea unc ion wi h slopes unde he c i ical alue 𝜇𝑠. Conc e ely, le us conside he piece-wise linea unc ion connec ing
nodes wi h abscissas 𝑋𝑗 o 𝑗=1,..., 𝑀. Deno ing by
𝛿=min
𝑗=2,...,𝑀(𝑋𝑗+1 −𝑋𝑗),
i Δ𝑥<𝛿∕3, which means ha he e a e a leas h ee consecu i e cells wi h he same slope, hen we locally find wo connec ed
linea polynomials as depic ed in Fig. A.9. The e o e, all ha has been said p e iously emains alid. In pa icula , in he nume ical
simula ions we obse e ha his ini ially unc ion is usually p ese ed wi h he excep ion o he joining cells. In hose cells, he ini ial
condi ion is modified du ing he fi s i e a ions and hen a disc e e s eady s a e is ob ained which is an app oxima ion o he o iginal
unc ions. Tha is, he scheme is well-balanced, bu no exac ly well-balanced. This beha io will be analyzed in he nume ical es s
sec ion (see Subsec ion 4.2).
Re e ences
[1] P. Jop, Y. Fo e e, O. Pouliquen, A cons i u i e law o dense g anula flows, Na u e 441 (2006) 727–730.
[2] P.-Y. Lag ée, L. S a on, S. Popine , The g anula column collapse as a con inuum: alidi y o a wo-dimensional Na ie -S okes wi h a 𝜇(I)- heology, J. Fluid Mech.
686 (2011) 378–408.
[3] C. Lusso, A. E n, F. Bouchu , A. Mangeney, M. Fa in, O. Roche, Two-dimensional simula ion by egula iza ion o ee su ace iscoplas ic flows wi h D ucke -
P age yield s ess and applica ion o g anula collapse, J. Compu . Phys. 333 (2017) 387–408.
[4] I.R. Ionescu, A. Mangeney, F. Bouchu , R. Roche, Viscoplas ic modeling o g anula column collapse wi h p essu e-dependen heology, J. Non-New on. Fluid
Mech. 219 (2015) 1–18.
[5] N. Ma in, I.R. Ionescu, A. Mangeney, F. Bouchu , M. Fa in, Con inuum iscoplas ic simula ion o a g anula column collapse on la ge slopes: 𝜇(I) heology and
la e al wall effec s, Phys. Fluids 29 (2017) 013301.
[6] O. Pouliquen, Y. Fo e e, A non-local heology o dense g anula flows, Philos. T ans. R. Soc. A, Ma h. Phys. Eng. Sci. 367 (2009) 5091–5107.
[7] F. Boye , É. Guazzelli, O. Pouliquen, Uni ying suspension and g anula heology, Phys. Re . Le . 107 (2011).
[8] M. Bouzid, M. T ulsson, P. Claudin, E. Clémen , B. And eo i, Nonlocal heology o g anula flows ac oss yield condi ions, Phys. Re . Le . 111 (2013).
[9] M. Rau e , The comp essible g anula collapse in a fluid as a con inuum: alidi y o a Na ie –S okes model wi h 𝜇(I)- heology, J. Fluid Mech. 915 (2021).
[10] S.B. Sa age, K. Hu e , The mo ion o a fini e mass o g anula ma e ial down a ough incline, J. Fluid Mech. 199 (1989) 177–215.
[11] J.M.N.T. G ay, A.N. Edwa ds, A dep h-a e aged 𝜇(I)- heology o shallow g anula ee-su ace flows, J. Fluid Mech. 755 (2014) 503–534.
[12] A.N. Edwa ds, J.M.N.T. G ay, E osion-deposi ion wa es in shallow g anula ee-su ace flows, J. Fluid Mech. 762 (2015) 35–67.
[13] E.D. Fe nández-Nie o, J. Ga es-Díaz, A. Mangeney, G. Na bona-Reina, A mul ilaye shallow model o d y g anula flows wi h he 𝜇(𝐼)- heology: applica ion o
g anula collapse on e odible beds, J. Fluid Mech. 798 (2016) 643–681.
[14] J. Ga es-Díaz, E.D. Fe nández-Nie o, A. Mangeney, T.M. de Luna, A weakly non-hyd os a ic shallow model o d y g anula flows, J. Sci. Compu . 86 (2021).
[15] C. Escalan e, E.D. Fe nández-Nie o, J. Ga es-Díaz, A. Mangeney, Mul ilaye shallow model o d y g anula flows wi h a weakly non-hyd os a ic p essu e, J. Sci.
Compu . 96 (2023).
[16] F. Bouchu , J.M. Delgado-Sánchez, E.D. Fe nández-Nie o, A. Mangeney, 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 flow models, Appl. Ma h. Model. 106 (2022) 627–658.
[17] E. Audusse, M. B is eau, B. Pe hame, J. Sain e-Ma ie, A mul ilaye Sain -Venan sys em wi h mass exchanges o shallow wa e flows. De i a ion and nume ical
alida ion, ESAIM: Ma h. Model. Nume . Anal. 45 (2010) 169–200.
[18] E.D. Fe nández-Nie o, E.H. Koné, T.C. Rebollo, A mul ilaye me hod o he hyd os a ic Na ie -S okes equa ions: a pa icula weak solu ion, J. Sci. Compu . 60
(2013) 408–437.
[19] S. Bosca ino, G. Russo, M. Semplice, High o de fini e olume schemes o balance laws wi h s iff elaxa ion, Compu . Fluids 169 (2018) 155–168.
[20] A. Ku gano , J. Mille , Cen al-upwind scheme o Sa age–Hu e ype model o subma ine landslides and gene a ed sunami wa es, Compu . Me hods Appl.
Ma h. 14 (2014) 177–201.
[21] A. Mangeney-Cas elnau, J.-P. Vilo e, M.O. B is eau, B. Pe hame, F. Bouchu , C. Simeoni, S. Ye neni, Nume ical modeling o a alanches based on Sain Venan
equa ions using a kine ic scheme, J. Geophys. Res., Solid Ea h 108 (2003), h ps://agupubs.onlinelib a y.wiley.com/doi/pd /10.1029/2002JB002024.
[22] A. Be mudez, M.E. Vazquez, Upwind me hods o hype bolic conse a ion laws wi h sou ce e ms, Compu . Fluids 23 (1994) 1049–1071.
[23] J.M. G eenbe g, A.Y. Le oux, A well-balanced scheme o he nume ical p ocessing o sou ce e ms in hype bolic equa ions, SIAM J. Nume . Anal. 33 (1996)
1–16.
[24] L. Gosse, A well-balanced flux- ec o spli ing scheme designed o hype bolic sys ems o conse a ion laws wi h sou ce e ms, Compu . Ma h. Appl. 39 (2000)
135–159.
[25] E. Audusse, F. Bouchu , M. B is eau, R. Klein, B. Pe hame, A as and s able well-balanced scheme wi h hyd os a ic econs uc ion o shallow wa e flows, SIAM
J. Sci. Compu . 25 (2004) 2050–2065.
[26] M.J. Cas o Díaz, J.A. López-Ga cía, C. Pa és, High o de exac ly well-balanced nume ical me hods o shallow wa e sys ems, J. Compu . Phys. 246 (2013)
242–264.
[27] C. Be hon, C. Chalons, A ully well-balanced, posi i e and en opy-sa is ying Goduno - ype me hod o he shallow-wa e equa ions, Ma h. Compu . 85 (2015)
1281–1307.
Applied Nume ical Ma hema ics 215 (2025) 138–156
156
M.J. Cas o Díaz, C. Escalan e, J. Ga es-Díaz e al.
[28] M.J. Cas o, T. Mo ales de Luna, C. Pa és, Well-balanced schemes and pa h-conse a i e nume ical me hods, in: R. Abg all, C.-W. Shu (Eds.), Handbook o
Nume ical Analysis, in: Handbook o Nume ical Me hods o Hype bolic P oblems Applied and Mode n Issues, ol. 18, Else ie , 2017, pp. 131–175.
[29] E.D. Fe nández-Nie o, F. Bouchu , D. B esch, M.J. Cas o Díaz, A. Mangeney, A new Sa age-Hu e ype model o subma ine a alanches and gene a ed sunami,
J. Compu . Phys. 227 (2008) 7720–7754.
[30] J. Zhai, W. Liu, L. Yuan, Sol ing wo-phase shallow g anula flow equa ions wi h a well-balanced noc scheme on mul iple gpus, Compu . Fluids 134–135 (2016)
90–110.
[31] E.D. Fe nández-Nie o, J. Ga es-Díaz, A. Mangeney, G. Na bona-Reina, 2D g anula flows wi h he 𝜇(𝐼) heology and side walls ic ion: a well-balanced mul ilaye
disc e iza ion, J. Compu . Phys. 356 (2018) 192–219.
[32] F. Campos, M. Sepúl eda, R. Aba ca del Río, D. Issle , Es udio de modelos de a alancha usando esquemas de olúmenes fini os bien balanceados, Ob as y
P oyec os, 2023.
[33] M. Sanz-Ramos, E. Bladé, P. Olle , G. Fu dada, Nume ical modelling o dense snow a alanches wi h a well-balanced scheme based on he 2D shallow wa e
equa ions, J. Glaciol. (2023) 1–17.
[34] E. Fe nández-Nie o, J. Ga es-Díaz, P. Vigneaux, Mul ilaye models o hyd os a ic He schel-Bulkley iscoplas ic flows, Compu . Ma h. Appl. 139 (2023) 99–117.
[35] F. Bouchu , Nonlinea S abili y o Fini e Volume Me hods o Hype bolic Conse a ion Laws, Bi khäuse , Basel, 2004.
[36] M.J. Cas o, C. Pa és, Well-balanced high-o de fini e olume me hods o sys ems o balance laws, J. Sci. Compu . 82 (2020).
[37] M.J. Cas o Diaz, I. Gómez Bueno, C. Pa és, Hype bolic P oblems: Theo y, Nume ics, Applica ions. Volume I, Sp inge , 2022.
[38] C. Pa és, Nume ical me hods o nonconse a i e hype bolic sys ems: a heo e ical amewo k, SIAM J. Nume . Anal. 44 (2006) 300–321.
[39] M.J. Cas o Díaz, E.D. Fe nández-Nie o, A class o compu a ionally as fi s o de fini e olume sol e s: PVM me hods, SIAM J. Sci. Compu . 34 (2012)
A2173–A2196.
[40] I. Gómez-Bueno, M.J.C. Díaz, C. Pa és, G. Russo, Colloca ion me hods o high-o de well-balanced me hods o sys ems o balance laws, Ma hema ics 9 (2021)
1799.
[41] C.A. Kennedy, M.H. Ca pen e , Addi i e Runge–Ku a schemes o con ec ion–diffusion– eac ion equa ions, Appl. Nume . Ma h. 44 (2003) 139–181.
[42] L. Pa eschi, G. Russo, Implici -explici Runge-Ku a schemes and applica ions o hype bolic sys ems wi h elaxa ion, J. Sci. Compu . 25 (2005) 129–155.
[43] S. Bosca ino, F. Filbe , G. Russo, High o de semi-implici schemes o ime dependen pa ial diffe en ial equa ions, J. Sci. Compu . 68 (2016) 975–1001.
[44] S. Go lieb, C.-W. Shu, To al a ia ion diminishing Runge-Ku a schemes, Ma h. Compu . 67 (1996).
[45] X. Zhang, C.-W. Shu, Maximum-p inciple-sa is ying and posi i i y-p ese ing high-o de schemes o conse a ion laws: su ey and new de elopmen s, P oc.
Royal Soc. A, Ma h. Phys. Eng. Sci. 467 (2011) 2752–2776.
[46] I. C a e o, M. Semplice, G. Viscon i, Op imal defini ion o he nonlinea weigh s in mul idimensional cen al WENOZ econs uc ions, SIAM J. Nume . Anal. 57
(2019) 2328–2358.