A Simplex Multi-Phase Approach for Modelling Debris Flows in Smoothed-Terrain-Following Coordinate System

. 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 the dynamical variation of the composition concentrations. In the present study the component of silt/clay is an individual 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 the viewpoint 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 the interstitial fluid. The grain constituent is treated as a frictional Coulomb-like continuum, and the viscosity of the interstitial depends on the clay concentration. The resultant models are given in a smoothed-terrain-following coordinate system, a compromise between the constraint of shallow curvature for the terrain-fitting 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 and applicability will be presented by back calculation of a historical event.


Introduction
For the sake of hazard assessment or risk analysis, scenario investigations by numerical simulation are often employed for the quantitative evaluation. The present study aims at developing a multi-phase model for debris flows, which arise in mountainous area and sliding down rugged topography. The body of debris flows is generally a mixture of fluid, fine sediment and grains. Because of the complex composition, the rheology of debris flows is certainly different to the conventional Newtonian one, and it is the most challenging task in the modelling work. In the frame work of continuum mechanical approach, there are various 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-phase models [8]. Each type has various pros and cons with respect to the aspect of applications.
In the present study, the flow body is supposed to be composed of water, fine sediment (clay/silt) and grains. The component of silt/clay is treated as an individual phase, and its concentration plays a key role in 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 the momentum 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 a frictional Coulomb-like continuum, and the viscosity of the interstitial depends on the clay concentration.
Probably hindered by the complicated mathematical formulation, 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 number of 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, the definitions of the flow thickness (depth) and the direction 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. In general, the flow velocity is along the basal surface, so that a significant deviation of flow direction may take place for the depth-averaged velocity for a highly rugged topography (cf. ux in Fig. 1). Furthermore, the 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 of the slope should be small or with mild variation. The terrain-following coordinate system is therefore proposed, in which the flow thickness is measured in the normal direction and the resultant mean velocity is parallel 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 the digital elevation model (DEM), a "sub"-topography ζB is introduced to sit on the smoothed surface, and its magnitude fulfils the mass conservation. As shown in Fig. 3, they can be related by B B z x ζ ξ ∆ = ∆ ∆ (from the view point of mass conservation), where ∆zB is the elevation difference due to the smoothing process, ∆x and ∆ξ are the mesh size in the Cartesian and terrainfollowing coordinate system, respectively.

Smoothed-Terrain-Following Coordinates
In the present study, the topographic elevation is smoothed by a 3×3 Gaussian smoothing filter, . (2) The smoothing process by (2) several times until the maximal curvature (minimal curvature radius) meet the shallow curvature requirement.
Once the smoothed topographic surface is determined, a local coordinate system (ξ, η, ζ) is introduce, where the projections of ξand η-axis on the horizontal plane are assigned to coincide with the X-and Y-axis of the Cartesian coordinates [9,10]. The unit normal vector is determined by (3) and the coordinate transformation between the Cartesian and terrain-following coordinates reads with 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, , .

Two-phase approach of three constituents for debris flows
Through the scale analysis with shallowness assumption for the flow body, the depth-averaged are obtained in the non-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 for 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 of momentum conservation read and tan , .

Application to Historical Event
This simplex multi-phase approach is applied to a historical 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 Mm 3 ) is released, flowing over 2.8 km down to the Chishan river. The DEM and initial released mass are identical to the ones used in [12], where the DEM resolution is 10 m. The water level in the river follows the in-situ record in the nearest water level monitoring station. Figure 6 illustrates the mass movement and the merge into the river water, where there is distinct difference of clay content between the landslide body and river water. As shown in the results at time level t = 181.8 sec, two dammed lakes with depth more than 20 m can be recognized due to the deposits.  Fig. 6. The movement of the released mass, and the merge into the river water, where the red outline is the Hsiaolin village.
As shown in Fig. 7, the discrepancy is less than 4.7% between the flow paths computed in the coordinates with smoothing operation of 2 and 40 times, where computations are performed without the water in the river for a thorough investigation. Figure 8 shows that the topography has been significantly smoothed, and the DEM is reproduced by adding sub-topography on the smoothed surface (lower panel). It indicates the applicability of the smoothed-terrain approach.

Concluding Remarks
In the present study a simplex multi-phase approach is proposed for modeling flows, where the merge of flows with various viscosities for the interstitial fluid can be well described. The model equations are given in a smoothed-terrain-following coordinate system, so that the 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 even with 40 times of smoothing operations.  The financial support of the National Science and Technology Council, Taiwan (MOST 110-2221-E-006-045-MY2) for this research is sincerely acknowledged.