E3S Web of Conferences 415, 02022 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502022

A Simplex Multi-Phase Approach for Modelling Debris Flows in
Smoothed-Terrain-Following Coordinate System
Yih-Chin Tai *, Hock-Kiet Wong, and Ching-Yuan Ma
Department of Hydraulic and Ocean Engineering, National Cheng Kung University, 70101 Tainan City, Taiwan
Abstract. Herewith we present a multi-phase model for debris flows, of which the flow body is supposed
to be composed of water, fine sediment (clay/silt) and grains. The rheology of debris flows varies due to thedynamical variation of the composition concentrations. In the present study the component of silt/clay is anindividual phase, and its concentration plays a key role in determining the rheology of the interstitial fluid.
Hence, there are three phases in the mixture, the grain phase, the clay phase and the water phase from theviewpoint of mass conservation. Only the grain phase and fluid phase are considered in the momentum
conservation, since the clay is suspended in the fluid and the relative motion is negligible within theinterstitial fluid. The grain constituent is treated as a frictional Coulomb-like continuum, and the viscosityof the interstitial depends on the clay concentration. The resultant models are given in a smoothed-terrainfollowing coordinate system, a compromise between the constraint of shallow curvature for the terrainfitting coordinate system and retaining the high resolution of the topography. The numerical implementation
is developed with the CUDA-library for GPU-high-performance computations. The feasibility andapplicability will be presented by back calculation of a historical event.

1 Introduction
For the sake of hazard assessment or risk analysis,
scenario investigations by numerical simulation areoften employed for the quantitative evaluation. The
present study aims at developing a multi-phase modelfor debris flows, which arise in mountainous area andsliding down rugged topography. The body of debris
flows is generally a mixture of fluid, fine sediment andgrains. Because of the complex composition, therheology of debris flows is certainly different to theconventional Newtonian one, and it is the mostchallenging task in the modelling work. In the frame
work of continuum mechanical approach, there arevarious types of theoretical models, such as the singlephase models [1,2], quasi-two-phase models [3,4],
solid-fluid two-phase models [5,6,7], and multi-phasemodels [8]. Each type has various pros and cons withrespect to the aspect of applications.
In the present study, the flow body is supposed tobe composed of water, fine sediment (clay/silt) andgrains. The component of silt/clay is treated as anindividual phase, and its concentration plays a key rolein determining the rheology of the interstitial fluid.
Hence, there are three phases, the grain phase, the clay
phase and the water phase for mass conservation. Only
the grain phase and fluid phase are considered in themomentum conservation, because the clay is suspended
in the fluid and the relative motion is negligible within
the interstitial fluid. The grain constituent is treated as africtional Coulomb-like continuum, and the viscosity of
the interstitial depends on the clay concentration.

Fig. 1. The depth-averaged velocity ux (in Cartesian coordinate)
and uξ (in terrain-following coordinate).

Probably hindered by the complicated mathematicalformulation, most of the model equations for hazardous
flows are given in Cartesian coordinates [1,2,3,5].
Moreover, in modelling the mass movement, the flow
thickness is assumed to be shallow, and depthintegration along the flow thickness is utilized for
reducing the complexity of computation and the numberof model equations. The employment of the depthintegration in the Cartesian coordinate system arouses
several crucial concerns in modelling the flows on the
non-trivial topographic surfaces. Namely, thedefinitions of the flow thickness (depth) and thedirection of the flow velocity. In the Cartesian
coordinates, the flow thickness is generally defined in
the vertical z-direction and the depth-averaged flow
velocity is parallel to the horizontal flat xy-plane. Ingeneral, the flow velocity is along the basal surface, sothat a significant deviation of flow direction may takeplace for the depth-averaged velocity for a highly
rugged topography (cf. ux in Fig. 1). Furthermore, the* Corresponding author: yctai@ncku.edu.tw© The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0(https://creativecommons.org/licenses/by/4.0/).

E3S Web of Conferences 415, 02022 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502022

flow depth is generally counted along the normal
direction of the basal surface. Notable deviation of flow
thickness is encountered for deep slopes, cf. Fig. 2.
Consequently, the depth-averaged models are applicable
only for gentle topographies, i.e. the inclination angle ofthe slope should be small or with mild variation.

zismoothed= 0.0751zi −1, j +1 + 0.1238 zi , j +1 + 0.0751zi +1, j +1 +,j0.1238 zi −1, j + 0.2044 zi , j + 0.1238 zi +1, j +. (2)0.0751zi −1, j −1 + 0.1238 zi , j −1 + 0.0751zi +1, j −1

The smoothing process by (2) several times until the
maximal curvature (minimal curvature radius) meet theshallow curvature requirement.
Once the smoothed topographic surface isdetermined, a local coordinate system (ξ, η, ζ) is
introduce, where the projections of ξ- and η-axis on thehorizontal plane are assigned to coincide with the X- andY-axis of the Cartesian coordinates [9,10]. The unit
normal vector is determined by   n =gζ =gξ × gη gξ × gη =nX eˆ X + nY eˆY + nZ eˆZ , (3)Fig. 2. The flow depth h (counted in Cartesian coordinate)versus the hξ (counted in terrain-following coordinate).

and the coordinate transformation between the Cartesianand terrain-following coordinates reads(4)dr = Ω dξThe terrain-following coordinate system is therefore
proposed, in which the flow thickness is measured in thenormal direction and the resultant mean velocity isparallel to the basal surface (see uξ in Fig. 1). However,in addition to the complex mathematical formulation,
remarkable deviations arise in calculating the flow depth,
when the local curvature is large. And several negligible
terms in the model equations become significant. Hence,
the application of the terrain-following coordinate
system is with constraint of shallow curvature.
The present study proposes a strategy of introducing
a general formulation for model equations in an adaptive
terrain-following coordinate system, where highly
rugged topography may exist. That is, the coordinate
axes are defined by the smoothed topographic elevation.
To reproduce the “original” topography as given in thedigital elevation model (DEM), a “sub”-topography ζB
is introduced to sit on the smoothed surface, and its
magnitude fulfils the mass conservation. As shown inFig. 3, they can be related by ζ B ∆ξ =∆z B ∆x (from thewith the transformation matrix Ω. For any scalar a,vector u and symmetric tensor T, the following relations
have been used in derivation of the model equations,Ω −1 = A = ( Aij ) , ∇a = AT ∇ξ a  , J b = det Ω ,J b∇ ⋅ u = ∇ξ ( J b Au ) ,J b∇ ⋅ T = ∇ξ ⋅ ( J b TA ) .

(5)

T

2.2 Two-phase approach of three constituentsfor debris flows
Through the scale analysis with shallowness assumption
for the flow body, the depth-averaged are obtained in thenon-conventional approach (see [7,9,10] for details).
Letting h be the flow depth, v be the depth-averaged
velocity, the leading-order, dimensionless model
equations read∂( J b hϖ ) + ∂∂ξ { J b hϖ vξϖ } + ∂∂η { J b hϖ vηϖ } =J b Eϖ (6)∂tfor mass conservation, where E denotes the entrainment
rate, the super script ϖ = s, w, c stands for the solid,
water and clay constituent, respectively.
The leading-order, dimensionless equations ofmomentum conservation readview point of mass conservation), where ∆zB is theelevation difference due to the smoothing process, ∆xand ∆ξ are the mesh size in the Cartesian and terrainfollowing coordinate system, respectively.

∂∂J b h s v Xs ) + { J b h s v Xs vξs + ε J b hA11 N s }(∂t∂ξ∂+{ J b hs vXs vηs + ε J b hA21 N s }∂η∂=J b p nX − εα ρϕ  A11( J b h p f ) + A21 ∂∂η ( J b h p f )
∂ξ
sv+ α ρϕ sϕ f cD J b h ( v Xf − v Xs ) − Xs,b J b pbs tan δ bvbsb

Fig. 3. The smoothed coordinate (red dashed curve) and therugged topography (blue solid line).

s− ε J b ( N bs + α ρϕ s pbf ) Φ1121

and2 Model equations
2.1 Smoothed-Terrain-Following Coordinates
In the present study, the topographic elevation issmoothed by a 3×3 Gaussian smoothing filter,2∂( J b hs vYs ) + ∂∂ξ { J b hs vYs vξs + ε J b hA12 N s }∂t∂+{ J b hs vYs vηs + ε J b hA22 N s }∂η

(7)

E3S Web of Conferences 415, 02022 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502022(∂∂
J b pbs nY − εα ρϕ s  A12J b h p f ) + A22J b h p f )=((∂η ∂ξs(8)v+ α ρϕ sϕ f cD J b h ( vYf − vYs ) − Y s,b J b pbs tan δ bvbµ f =µ w 4.8α fc + 1 + (α fc ) 2viscosity. They are assumed to be given byφc µ f =cD 0.5cD ,ref 1.0 + c,φref µreff 

for the solid phase, where∂ζ B∂ζ ∂ζ B∂ζ B Φ=+ A21 B  , Φ=+ A2211211222
 A11
 A12.
∂ξ∂η ∂ξ∂η 
(11)(12)0.1N
 α fcϑbf ϑbf,ref 1.0 + R ,ref =.
N R  α reffc
In (11),=cD ,ref 6.0,
=φrefc 0.13125,=
φreff 0.5,=
ϑbf,ref 5.0,

are defined. In (7) and (8), N and p are the meanfsolid and fluid normal stress, and pbs is the solidpressure at bottom. For the fluid phase, the resultant
momentum equations read∂( J b h f vXf ) + ∂∂ξ { J b h f vXf vξf + ε J b hA11 p f }∂t∂+{ J b h f vXf vηf + ε J b hA21 p f }∂η7.54

where µw = 0.001 Pa ⋅ s is the viscosity of the water andαfc stands for the clay concentration in the interstitialfluid. Parameters cD and ϑbf are functions of the− ε J b ( N bs + α ρϕ s pbf ) Φ1222s) ,N R ,ref = 268 are the reference values for reproducing theresults by the two-phase approach in [7].

(9)∂=J b pbf nX + εϕ s  A11( J b h p f ) + A21 ∂∂η ( J b h p f )
 ∂ξfffϑ vεϕ− ϕ sϕ f cD J b h ( v Xf − v Xs ) +FXvis − J b h f b X
ε NRNR− J b h f α Mannn 2 vXf v f(h )f4/3(1 − ϕ ) − ε J ϕ p Φcbffb1121

Fig. 4. The viscosity of the interstitial fluid (cyan curve) andthe measured values (black markers) of in-situ samples against
the concentration, where the red star indicates the viscosity
used in [7]. (reproduced from [11])and∂∂J b h f vYf ) + { J b h f vYf vξf + ε J b hA12 p f }(∂t∂ξ∂+{ J b h f vYf vηf + ε J b hA22 p f }∂η

3 Application to Historical Event(10)∂∂=J b pbf nY + εϕ s  A12J b h p f ) + A22J b h p f )((∂η ∂ξfff
ϑ v
εϕ− ϕ sϕ f cD J b h ( vYf − vYs ) +FYvis − J b h f b Yε NRNR− J b h f α Mannn 2 vYf v f(h f )4/3This simplex multi-phase approach is applied to ahistorical event, the 2009 Hsiaolin event, which took
place during the Typhoon Morakot in 2009 in southern
Taiwan. A huge amount of mass (ca. 2.4±2 Mm3) isreleased, flowing over 2.8 km down to the Chishan river.
The DEM and initial released mass are identical to theones used in [12], where the DEM resolution is 10 m.
The water level in the river follows the in-situ record inthe nearest water level monitoring station.
Figure 6 illustrates the mass movement and themerge into the river water, where there is distinct
difference of clay content between the landslide bodyand river water. As shown in the results at time level t =181.8 sec, two dammed lakes with depth more than 20m can be recognized due to the deposits.

(1 − ϕ ) − ε J ϕ p Φcbffb1222

where NR = ρf ℋ ( 𝑔𝑔ℒ )1/2/µf with 𝑔𝑔 the gravity
acceleration, and ε = ℋ/ℒ << 1 the aspect ratio. In (9)and (10), the viscous drags read
FXvis,Y =2∂ξ ,η  J b h ( A11,12 ∂ξ vξf,η +A21,22 ∂η vξf,η ) + ∂η ,ξ  J b h ( A12 ∂ξ vξf +A22 ∂η vξf + A11∂ξ vηf +A21∂η vηf ) In (7) to (10), the terms in blue represent the contributionof the sub-topography.
2.3 Material parameters
Different to the solid-fluid model in [7], the density andthe viscosity of the interstitial fluid do not remain
constant. They vary along with the variation of the clay
concentration. Referring to various measurement of insitu samples, a clay-concentration-dependent relation
for the viscosity of the interstitial fluid (see the cyan
curve in Fig. 4) is suggested

3E3S Web of Conferences 415, 02022 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502022

Fig. 6. The movement of the released mass, and the merge into
the river water, where the red outline is the Hsiaolin village.
Fig. 8. Upper: Topography after smoothing operation for 40times. Lower: the added sub-topography to recover the DEM.

As shown in Fig. 7, the discrepancy is less than
4.7% between the flow paths computed in thecoordinates with smoothing operation of 2 and 40 times,
where computations are performed without the water inthe river for a thorough investigation. Figure 8 showsthat the topography has been significantly smoothed,and the DEM is reproduced by adding sub-topographyon the smoothed surface (lower panel). It indicates theapplicability of the smoothed-terrain approach.

The financial support of the National Science and TechnologyCouncil, Taiwan (MOST 110-2221-E-006-045-MY2) for thisresearch is sincerely acknowledged.

References
C.L. Chen, J. Hyd. Eng. 114, 237-257 (1988)
K.F. Liu, C.C. Mei, J. Fluid Mech. 207, 505-529(1989)
3. S. Egashira, K. Miyamoto, T. Itoh, Proceedings ofthe 1st International Conference on Debris-flow
Hazards Mitigation, C.L. Chen (Eds.), ASCE:New York: 340-349 (1997)
4. R.M. Iverson, R.P. Denlinger, J. Geophys. Res.,
106, 537-552. (2001)
5. E.B. Pitman, L. Le, Phil. Trans. R. Soc. A, 363,1573-1601.(2005)
6. S.P. Pudasaini, J. Geophys. Res., 117, F03010(2012)
7. Y.C. Tai, J. Heß, Y. Wang, J. Geophys. Res. 124,305-333 (2019).
8. S.P. Pudasaini, M. Mergili, J. Geophys. Res. 124,2920-2942 (2019)
9. Y.C. Tai, C.Y. Kuo, W.H. Hui, Geophys.
Astrophys. Fluid Dyn. 106(6), 596–629 (2012).
10. I. Luca, Y.C. Tai, C.Y. Kuo, Shallow Geophysical
Mass Flows down Arbitrary Topography (SpringerVerlag, 2016)
11. Takahashi, T., Debris Flow: Mechanics, Predictionand Countermeasures (Taylor & Francis, 2007).
12. C.Y. Kuo, Y.C. Tai, C.C. Chen, K.J. Chang, A.Y.
Siau, J.J. Dong, ... & C.T. Lee, J. Geophy. Res.
116, F04007 (2011)1.
2.

4 Concluding Remarks
In the present study a simplex multi-phase approach isproposed for modeling flows, where the merge of flows
with various viscosities for the interstitial fluid can bewell described. The model equations are given in asmoothed-terrain-following coordinate system, so thatthe employment for highly rugged topography without
the constraint of the shallow curvature is possible.
Numerical investigation of the Hsiaolin event reveals
the high potential of engineering application, that only
minor discrepancy is identified for the flow paths evenwith 40 times of smoothing operations.

Fig. 7. The discrepancy of flow paths between viscosity of theinterstitial fluid (cyan curve) and the measured values (blackmarkers) of in-situ samples against4