Debris flows with pore pressure and intergranular friction on rugged topography

Computers and Fluids 190 (2019) 139–155

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Debris flows with pore pressure and intergranular friction on rugged topography Julian Heß a,∗, Yih-Chin Tai a, Yongqi Wang a TU Darmstadt,Otto-Berndt-Str. 2, Darmstadt, 64287 Germany

a r t i c l e

i n f o

Article history: Received 18 December 2018 Revised 14 May 2019 Accepted 14 June 2019 Available online 15 June 2019 Keywords: Debris flow Hypoplasticity Pore-fluid pressure Terrain-following coordinates Two-phase flow

a b s t r a c t The dynamic behavior of debris flows features the interplay of a non-hydrostatic pore-fluid pressure with the non-linear deformational behavior of the granular skeleton and the internal contact stress between grains. This complex physical background is considered by amending the classical depth-integrated modeling for granular-fluid flows by two additional fields, an extra pore-fluid pressure and a hypoplastic intergranular stress. A scaled and depth-integrated model is developed and transferred into a system of terrain-following coordinates, enabling the application on rugged topography. With this model, numerical investigations are carried out, using a non-oscillatory, shock-capturing central-upwind scheme. Parameter studies show the general impact of the additional fields, completed by comparison to the experimental results of a dam break scenario. Furthermore, application to the landslide event at the village of Hsiaolin in Taiwan, 2009, show the capability of the model to cope with large scale scenarios. The results show that the model and its implementation provide insights in the flow dynamics and the possibility to application on complex topography, considering an enhanced approach to the physics of debris flows. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Both in industrial processes and nature, a variety of granularfluid flows can be observed, in which a mixture of granular solids and fluid components is set into motion. In general, one can distinguish flows that are not affected remarkably by the interstitial fluid, in this case often air, in their flow dynamics, from those that are dominated by the fluid phase. In the latter, granular particles tend to be dispersed and rarely in frictional contact, like in the case of mud flows. In between are debris flows, as collections of solid particles, forming a structure in which the interstitial fluid is of significant influence on the flow and the material behavior. Due to this, with debris flows, one cannot neglect the influence of the pore-fluid, i.e., its importance in the process of generating a mass flow with both solid-like and fluid-like behavior. And distinct from dry granular avalanches as well as mud flows, the increased mobility of debris flows can be linked to this interstitial fluid in its interdependency with the granular material. Occurring in mountainous areas, debris flows are often induced by heavy rainfall and tend to mobilize large masses of soil and rocks. Often tons of debris are accelerated and transported large

∗

Corresponding author. E-mail address: [email protected] (J. Heß).

https://doi.org/10.1016/j.compfluid.2019.06.015 0045-7930/© 2019 Elsevier Ltd. All rights reserved.

distances, sometimes causing destruction to property and casualties on a large scale. In December 1999, a storm hit the northern coast of Venezuela, mostly the state of Vargas, and triggered due to intense rainfall a series of debris flows, burying settlements and killing estimated 30,0 0 0 inhabitants, see [1,2]. Another tragic and more recent example is the landslide that occurred in the mountains of southern Taiwan in August 2009, erasing the village of Hsiaolin, again after heavy rainfall during typhoon Morakot, see [3,4]. All of this emphasizes the need for appropriate hazard management and countermeasures, which themselves rely on sound mathematical modeling, together with effective numerical computation. To face the task of modeling such material behavior and predicting flows accurately, various efforts have been made, particularly in the past 30 years. Many of those attempts were based on the pioneering work of [5], in which a scaled and depth-integrated model was proposed to provide equations describing the movement of a dry, granular flow. This Savage-Hutter (SH) approach has been taken up and extended in various directions, considering further phases, see [6], and various topographies, as in [7,8], forming basically a class of depth-integrated models for granular(-fluid) flows. While the two-phase approach of [6] has been further developed in [9,10], advancement for the treatment of a most general topography has been achieved with the introduction of general coordinates. Their implementation is a response to the difficulties of complex, rugged topography being treated in the framework

140

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

of a depth-integrated model and originates in the connection of the so-called unified coordinate method of Hui et al. see [11] and [12], with the system of terrain-following coordinates of [13]. For a detailed overview, covering the transfer of model equations to topography-following coordinates, see [14]. The complex interaction of the constituents also makes grasping the apparent physical mechanisms a difficult task. In the original SH-model, following the Coulomb theory for dry granular material, a rather simple earth-pressure correlation was adopted, coupling the driving, lateral pressure to the vertical stress via an earth pressure coefficient, see [5]. With an additional fluid phase, momentum interaction effects were proposed in [6,9,10], as well as viscous fluid behavior and sufficient boundary conditions for both components. Nonetheless, the complexity of a saturated granular structure, due to the interdependency of different effects, is not fully grasped in these works. The considerable momentum of the porefluid during onset and movement, resulting from pressure peaks, reduces the internal friction: Following an increased load on the granular structure, a rising pore pressure may, in return, press the granular material apart and thus destabilize the structure itself, increasing its mobility. With this, a stress supporting structure of granular grains may be carried over into a granular-fluid mass flow and already moving mass may accelerate further – and as the examples mentioned above point out, often heavy rainfall, saturating the pore space and inducing higher pore pressure, has been an apparent triggering mechanism. The evolution of the pore-fluid pressure in the context of debris flows is taken up by Iverson and coworkers, see [15–17], who developed a pressure-diffusion equation for the pore-fluid pressure, governed by dilatancy effects. Furthermore, the intergranular friction can be considered more profoundly by referring to the concept of hypoplasticity, introduced by [18], that develops further plasticity models for granular material. Additional fields of research not considered in the following, are, for example, concerned with mechanisms of particle-size segregation and soil erosion, see e.g. [19,20]. This work resumes and continues the derivation of a depthintegrated model for granular-fluid flows that considers pore-fluid pressure evolution and intergranular friction. Based in continuum mechanics, fundamental forms of the material laws for such a mixture have been deduced in a thermodynamically consistent manner and with recourse to mixture theory in [21]. Taking up these results and introducing them into the framework of the SHapproach, a set of scaled and depth-integrated equations has been derived for a 1D debris flow, applied to a curvilinear coordinate system in [22]. This system is now employed and enhanced to a general topography and tested extensively with parametric studies, the comparison to experimental results and its application to large scale events. These modeling efforts also take up a line of development that leads from the work of the Meng-Wang (MW) model, see [10], to an application on general topography in [4], since the current system can be seen as an amendment of the MW-model with additional fields, considering the extra pore-fluid pressure and intergranular contact stress. The structure of the article is as follows: Starting with some brief explanation of the modeling and the derivational process in Section 2, the incorporated physical mechanisms are further explained. We conclude with the numerical aspect of this work in Section 3, presenting the results of the parameter studies, comparison with experimental results and a case study of the Hsiaolin event, before closing with a conclusion in Section 4. 2. Modeling, derivation and physics In this section, an overview is given on the process of mathematical modeling, the applied system of equations and the incor-

Fig. 1. Overview of the derivational process of a scaled and depth-integrated model of the Savage-Hutter type, applied to general coordinates for a rugged topography. Starting with the general model for a system with an unspecified number of constituents α , successively, the scaled and depth-integrated equations for a granularfluid system are derived and applied to general coordinates.

porated physical mechanisms. For reasons of clarity and comprehensibility, we refrain from recapitulating the full actual derivation but state the complete system in a most general, overseeable form before explaining the further steps of the derivation. This differs from the form of presentation as it has been given in [4,22], but since a complete and consistent depiction of the full derivation would go beyond the scope of this work, stretch the patience of the reader and, first of all, has in parts already been presented in the aforementioned works. We believe that it is advisable to shift the focus to the practical use of the model. It is noted that we resort to previous findings on the constitutive functions, derived in [21] with the entropy principle in its formulation by Müller and Liu, see [23,24] and, summarizing, [25]. In the framework of continuum mechanical mixture theory, this allows for the deduction of restrictions on the constitutive functions, based on general principles and the chosen system of balance equations. This not only grounds the derived constitutive functions thermodynamically, here especially the stress and momentum interaction, but also provides a way of looking into how a certain set of balance equations, comprising additional fields, manifests in the material behavior of the respective constituents. In the case of granular-fluid flows, this means that we do not need to apply, for example, the Coulomb theory in order to find a formulation for the solid stress tensor. In contrast, one can for example follow the impact of additional fields like the hypoplastic, intergranular stress on the shape of the solid stress tensor. As a further result worth noting, an additional buoyancy term appears in the momentum interaction due to the derivation. For Savage-Hutter type models, the process of scaling and depth-integration is usually applied to the system of equations to provide closure for the otherwise unknown stress tensors. Due to the derivation in the framework of the entropy principle, the stress tensors already exhibits a given structure and thus is not completely unknown. Nonetheless, the equations are simplified in the process of scaling, the computational efforts are reduced due to depth-integration and the hydrostatic pressure still need a formulation that can be found in conjunction with depth-integration. This process, depicted in Fig. 1, is outlined in the following. It should be noted that with these modeling attempts, to the best of the authors knowledge, for the first time a transfer from entropyprinciple based derivations to – usually rather a priori based – depth-integrated application and numerical simulation has been established in this field. Starting with a general system of equations for α constituents, in a first step, the mixture is confined to two phases, a fluid and

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

a granular one, and applied to simple, curvilinear coordinates, as in [10]. With that, the equations are non-dimensionalized via scaling assumptions, introducing characteristic measurements for flow thickness and length, H and L, respectively, and isolating insignificant terms, similar to the treatment that leads to shallow-water equations. These scaled equations are then integrated in the normal direction, assuming a blunt velocity profile and that all quantities are not depending on the normal coordinate. At this time, material laws and further relations are introduced for the treatment of the stress, the hydrostatic fluid pressure, the momentum interaction and the bed friction terms. For a simple, curvilinear 1D system, this has been done in [22]. Here, these equations are then transformed to general coordinates, as it has been done for the Meng-Wang-model in [4], providing a set of scaled and depthintegrated equations for a two-phase debris flow on rugged topography, with evolution equations for the extra pore-fluid pressure and hypoplastic stress. So in this work, the model developed in [22] is extended to 2D and applied to rugged topography, working with general coordinates, where especially the transfer of the evolution equations for pore pressure and hypoplasticity into terrain-following coordinates is a new derivation. Due to the length of a complete illustration of the derivation, and also just the final equations, we refer to [4] and the Appendix A. 2.1. Equations, coordinates and assumptions

the extra pore-fluid pressure the inner solid stress Zisj :

ef

and the evolution equation for

α = s, f ∂ν α ρ α vαi vαj ∂ Tiαj + − − ρ α ν α gi − mαi = 0; ∂xj ∂xj

Balance of momentum for constituent

∂ν α ρ α vα ∂t

i

∂ν α ∂t

+

∂ν α vα ∂ xi

i



   s  s s ρs f f f f =− − 1 h − Ce e − ρ − ρ δZ Zkk ν δi j + ρδZ Zisj , ρf   Ti jf = − hf + Cef ef ν f δi j + 2μ f Difj . (2) Tisj

The appearance of the additional fields, e and Zisj , is in conjunction with two material constants, the pore-fluid pressure paramef ter Ce , and one for the plastic stress, δ Z . Analogously to the fluid viscosity μf , we assume them to be constant material parameters, but they could as well be more sophisticated functions in future f investigations. In the viscosity-term of Ti j , the fluid strain tensor f

Di j = 12 (∂ j vi + ∂i v j ) is in place. The fluid hydrostatic pressure h arises here as well. For momentum interaction, a simple drag term is considered s , and a buoyancy term. Unlike the comwith a drag coefficient cD mon choice for the latter, as for example in [10], an additional term is present here due to the thermodynamically consistent derivation, so it yields f

msi =

f

  s f ρs f f s (ρ − ρ ) − 1 + ( 1 − ν )   h h ρ ρf   + cDs ν s ν f vif − vsi . ∂ν s ∂ xi

f



(3)

The boundary conditions do not take erosion or deposition into account and the flow surface is stress free. At the flow bottom, for the solid, Coulomb friction is applied in conjunction with a basal friction angle δ b , and for the fluid, Navier’s friction law is applied, so



(1a)

Evolution equation for the volume fraction (balance of mass) for constituent

via the gravitational constant gi , the spin tensor of the solid comf f penent, Wisj = 12 (∂ j vi − ∂i v j ), and the fluid diffusion coefficient κ f . The mixture density ρ is composed as follows, with ρ = ρ s ν s + ρfνf. The apparent material behavior is governed to a great extent by the stress tensor for the solid and the fluid, since it captures in which way the forces and inertia acting upon the mixture are resulting in the movement and thus deformation. Both stress tensors exhibit a spherical term, containing the pressure parts, i.e., the hydrostatic pressure and the fluid extra pressure, together with a plastic spherical term for the solid. Further, while the solid stress tensor comprises a plastic term, linked to the hypoplastic intergranular stress, the fluid stress tensor includes a viscous term, both anisotropic. With this, the stress tensors for solid and fluid are

f

After the path of derivation, i.e., the general steps in the treatment of the equations, has been depicted in the previous section, see Fig. 1 for an overview, the actual equations are now introduced together with assumptions and a closer look at the applied coordinates. A saturated two phase system is considered, consisting of a granular/solid (s) and a fluid (f) phase, without any phase transition or chemical reactions. As a further simplification, a constant mixture temperature is assumed, together with constant true densities ρ s , ρ f . The material is described by the balance of mass (for the volume fractions ν s , ν f ), the balance of momentum (for the vef locities of the individual phases vsi , vi ), an evolution equation for

141

α = s, f

= 0;

(1b)

Evolution equation for extra pore-fluid pressure

  ∂ef ∂ef f ∂ κ f ∂ 2 ef + −κf − γf = 0; vi − ∂t ∂ xi ∂ xi ∂ xi ∂ xi

(1c)

Evolution equation for the inner solid stress

 s  ∂ Zisj ∂vsk Zisj ρ − ρ f ∂vsk s s + + Zik Wksj − Wiks Zks j − ν s Z − si j = 0. ∂t ∂ xk ρ ∂ xk i j

(1d)

Here, the material behavior is described by the constituents stress tensor Tiαj , the momentum interaction mα , as well as the i plasticity source term si j and the pressure source term γ . Furthermore, there is the external influence of gravitation, effective f

vsi s T | tan δb , ||vs || n b    Ti jf  n j − Tnf  ni = kl ν f vif  , b b b Tisj  n j − Tns |b ni = b

(4)

where kl is the fluid friction coefficient. Here, the normal stresses s, f s, f Tn |b = Ti j |b ni n j are introduced, together with ni as the normal unit vector of the basal surface. So far, just the basic, dimensional equations have been introduced. After the process of scaling, the inherent physical parameters are transferred into non-dimensional numbers, which will be present in the following, i.e., in the numerical simulations. As in [22], with the introduction of those non-dimensional numbers, the contribution of different influences to the flow dynamics is described, namely of the additionally introduced fields. These are

C f ρ f gH Eu = e , ρ gL

 0  0 0 ρδZ Zxx + Zyy + Zzz NZ = , ρ gL

where, for the dimensionless intergranular friction, a number NZ is introduced with a referential intergranular contact stresses Zii0 . It gives a relation between the hypoplastic frictional stress terms and the inertial forces. The Euler number Eu specifies the influence of the extra pressure and appears in conjunction with the

142

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

Fig. 2. Comparison of a Cartesian coordinate system (X, Y, Z), curvilinear coordinates (x, y, z) following the references plane with a shallow elevation zb of the topography hereon, on the left-hand side and terrain-following coordinates (ξ , η, ζ ) on the right-hand side, depicting a general topographic surface S. While the curvilinear coordinates are suitable for simple chutes, occurring e.g. in an experimental setup, the general coordinates allow for the depiction of large mountainous areas.

non-dimensional field variable for the introduced extra pore-fluid f pressure, e . We also introduce Eu|b to denote the influence of the dynamic pore pressure on the solid bed friction, with

Eu|b =

Cef |b ρ f gH , ρ gL

to pay attention to the fact that the pore pressure elevation may not be uniform in the normal direction, see [17]. Besides, there is a non-dimensional drag number coefficient cD , the viscous number f NR and the scaled fluid bed friction ϑb , as introduced in [10] and [4], with

cD =

cDs

ρf

 g , NR = L

√

ρ f H gL k H2 , ϑbf = l f . f μ μ

As the scaling allows for the identification of physically relevant terms, giving a measurement of their impact on the flow dynamics and in relation to each other via the introduced non-dimensional numbers, the depth-integration reduces the computational efforts remarkably, and thus the complexity of the system of equations. The main assumption is that changes of all quantities in the normal direction can be neglected, so that the originally 3D system reduces to two dimensions, basically with a projection of all quantities on the basal plane. For an inclined plane, the direction normal to the basal surface may not change, so that a Cartesian coordinate system (X, Y, Z) is applicable, see [5]. In contrast, for a chute with transition to a runout, the (position-dependent) normal direction needs to be described by a curvilinear coordinate system (x, y, z), as introduced for the SH-modeling approach in [7] and extended in [8]. In other words, one needs to introduce a slightly more sophisticated coordinate-ansatz in order to describe flows on certain reference surfaces, that resemble, for example, simple laboratory experimental setups. With the introduction of curvilinear coordinates, certain terms enter the equations, accounting for the curvature of the basal surface. Nonetheless, there are some short-comings of curvilinearcoordinates when it comes to depicting the topography of real mountainous areas, since in such cases, there is not only one clearly distinguishable flow direction in place, but a multitude of rutted channels. Here, so-called general coordinates are capable of both describing the rugged topography by maintaining Cartesian coordinates, and, at the same time, a depth-integrated flow in conjunction with terrain-following coordinates (ξ , η, ζ ). This method, applied in [4,20,26], is possible due to the combination of unified coordinates, developed in [11,12], and the approach of [13]. It unites both the advantage of reduced equations with respect to the degrees of freedom, keeping computational costs low, but allows for the easy implementation of topography data, for example

from geographic information systems (GIS) and their digital elevation models. For a comparison of curvilinear, Cartesian and general coordinates, see Fig. 2. 2.2. Physical mechanisms: pore-fluid pressure and intergranular friction So far, the two additional fields and the development of their balance equations in terms of scaling, depth-integration and the treatment in the general coordinate system have not been discussed. As for scaling and depth-integration, we refer to [22], where this has been done in terms of curvilinear coordinates. The transfer of these equations to general coordinates is further discussed in the Appendix A, since we refrain here from giving explicitly the final equations in their lengthy form. Instead, in the following, the physical implications of these equations are discussed. The principal aim here is to regain information about physical mechanisms that are not depicted in a classical mixture theory model, since the individual particles, and thus their interaction, as well as between pore space and fluid, are not tracked. 2.2.1. Pore-fluid pressure Soil can be described as a granular-fluid mixture in a resting state, a network of pores, built by a granular material and filled – completely or in parts – with fluids, such as air and (more or less muddy) water. The contact between granular particles is of frictional nature and can be altered by the pore-fluid. This can be described in terms of the effective stress principle, marking a transition from load on the granular structure into increased porefluid pressure, which in return destabilizes the granular structure, as visualized in Fig. 3. The granular particles are themselves subject to the pressure of the pore-fluid, which may surround the single particles, as depicted in 3 b. At the same time, these particles act collectively as a structure that induces pressure upon the porefluid, increasing the pressure level of the fluid or even, if possible, squeezing it out of the material, see Fig. 3a. Such mechanisms are discussed in soil mechanics since the early works of Fillunger and Terzaghi on porous media, see [27] and, for an overview, [28]. With this, it is apparent that the pore-fluid pressure is not hydrostatic and should be modeled accordingly, since it is not only important for the mobilization of debris flows but also affecting the dynamic behavior, preventing the material, once in motion, from resting and, with this, prolonging the runout length. Therefore, the impact of the dynamic pore pressure evolution on the velocity of the granular mass and on their material distribution is apparent.

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

Fig. 3. Pore-fluid pressure mechanism during loading cycle. Panel a on the lefthand side illustrates how the granular skeleton transfers external forces into extra fluid pressure and, on the right-hand side in panel b, this pore-fluid pressure presses subsequently apart the granular material.

An approach of Iverson et al. seeks to describe this behavior, developed in a series of publications from early attempts in [15,29], to [16] and finally [17]. We take up the modeling work of [17], where a diffusion equation for the extra pore pressure is formulated in conjunction with a depth-integrated model for a granularfluid mass flow, describing the evolution of the pore-fluid pressure in response to the stress state and the deformational behavior of the granular material. In particular, this is done via the incorpof ration of dilatancy effects in the source term γ in Eq. (1c), thus stating 

γ = f

d σ −hf dt



− αγ˙D tan(ψ ),

(5)

where the first term induces changes in the stress state, and the second term accounts for changes in the pore space due to shearing, i.e., dilatancy, with the shear rate γ˙ and the tangent of the dilatancy angle ψ , which can be modeled as the difference of the solid volume fraction and its equilibrium value, tan(ψ ) = s ), with two coefficients κ κ 1 (ν s − κ 2 νeq ϖ1 and κ ϖ2 . A positive value of the extra pore pressure depicts a state in which the excess pressure unburdens the granular structure and reduces intergranular friction, as the fluid is pushed in between the granular particles. In contrast, a negative pore-fluid extra pressure means that the pore space is absorbing the fluid. So it is apparent, that the pore pressure evolution is connected to the pore space, which is resembled by the fluid volume fraction, respectively also the solid volume fraction in this saturated mixture. This is also captured in Eq. (5): If the actual value of the solid volume fraction s , the dilatancy is smaller than the equilibrium volume fraction νeq term, − tan(ψ ) ≥ 0, becomes positive (with γ˙ > 0, α D > 0). From s correspond to a this relation, it follows that large values of νeq positive extra pore-fluid pressure, since the granular structure allows for a more dense packing until the effects of dilatancy induce a shift to negative extra pore pressure. For the modeling of s =ν s (ν s , γ˙ , μ f , δ ), including values for the average particle νeq ˆ eq P C size δ P and the critical volume fraction νCs , we refer to the detailed derivation in [17]. Besides, the diffusion term in Eq. (1c) contains a diffusivity coefficient κ f that is determined as

κ = f

km

αD μ f

,

where, besides the fluid viscosity μf , the hydraulic permeability km = k0 exp ((νCs − ν s )/0.04 ) and the debris compressibility αD = aD /(ν s σ ) arise. Here, k0 and aD are constants of permeability and compressibility, respectively, νCs is, as before, the critical solid volume fraction, and σ denotes the stress state. This approach has

143

Fig. 4. Intergranular forwarding of load with respect to different micro-structures: Panel a on the left-hand side shows a dense granular skeleton that will redirect forces in a different way than the skeleton depicted on the right-hand side in b, which will likely reduce the pore space. Even more, since the structure in panel b will change under load, e.g. to a structure alike the one depicted in panel a, it will react in different way if the load is changing or reversed. Also see Fig. 3 in [36] and the companying explanation for more details.

been adopted and included in the modeling process of this work. A different pore pressure model is presented in [30]. 2.2.2. Intergranular friction Besides this interdependency of pore-fluid pressure and granular material, the behavior of (also just dry) granular material can be tackled with hypoplasticity. This concept reflects the dependence of the stress on the deformational history of the material, i.e that the material exhibits different paths of deformation during loading and unloading. Going back to [18,31], hypoplasticity enhanced previous approaches of elasto-plasticity and has been adopted in various works, see [32–34] and, for a comparison to elasto-plasticity, also [35]. Hypoplasticity establishes a model for the development of stress and contact forces in simple granular skeletons, see Fig. 4 and also [36]. The idea in Eq. (1d) is that the temporal evolution of a stress tensor, described with an objective, corrotational time-derivative, as introduced by Jaumann, can be balanced by a source term si j that comprises both a linear and a nonlinear part, L and N, respectively. While the first, linear part couples the temporal development of the stress to a product of the stress state itself and the strain rate, the nonlinear term introduces a rate-dependent behavior. Following the formulation in [37,38], the source term is employed as

  si j = L(Zisj , Dsi j ) + N(Zisj )Dsi j   s s s = fs

asZ Dsi j

+

Zi j Zkl Dkl s (Zmm )2

+

 fD asZ

|

Dsi j

|

Zisj

+

s Zmm

dev(Zisj ) s Zmm

 ,

(6)

where, the absolute value of the strain tensor, equivalent here to the strain rate, is |Dsi j | =

Dsik Dsk j δi j and the deviatoric part of the

s δ . As material intergranular stress is given as dev(Zisj ) = Zisj − 13 Zkk ij 1−ν s

parameters, a stiffness factor fs = fs0 fb ( 1−νCs ) pe0 is introduced, toνs

−ν s

gether with a density factor fD = fD0 ( νmax )e0 , see [32,39] and s −ν s max

[36]. Note that

asZ = a0Z

C

√ 3(3 − sin (φint ) ) , √ 2 2 sin (φint )

as in [36], and there are the material constants e0 = 0.1, with a range of 0.0 < e0 < 1.0 given in [32] and restricted to 0.1 < e0 < 0.3 in [39], while the other exponent pe0 should slightly exceed the value of 1; here, [39] suggest 1.0 < pe0 < 1.1, so we apply pe0 = 1.1. Furthermore, there is φ int as the internal friction angle and the

144

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

maximum packing fraction νmax . The stiffness factor fs includes a barotropy function fb , given as



fb =

Zii (asZ )2 −

 √ 3asZ

pn

,

which is slightly simplified compared to [32] and where pn , exhibiting a possible range of 0.3 < pn < 0.5, see [39], is set to pn = 0.3. There is a wide literature on the formulation of these parameter functions and the possible calibration of the material constants, see for instance [36,40] and, as mentioned before, [32,39]. It is apparent, that with introduction of several limiting values for the s and ν solid volume fraction, such as νCs , νeq max , one seeks to capture the developments of the microstructure of the granular material and its influence on the pore space, thus the fluid pressure and the friction between the grains. As with the formulation for the extra pore-fluid pressure, this granular material behavior is linked to Terzaghi’s principle of effective stress, see [36].

Fig. 5. A chute with curved bed, applied for parametric studies in Section 3.2. The actual topographic surface, indicated by the dashed lines, is superposed over a simple curved plane, denoted by the dotted lines.

3. Numerical studies After the derivation of a set of depth-integrated equations for the applied general coordinates, this system needs to be solved numerically. With this, it is the aim to look into the behavior of the mass flow, especially in conjunction with the additional fields of pore-fluid pressure and hypoplasticity. 3.1. System and numerical treatment The difficulties in the numerical treatment of depth-integrated equations for granular-fluid flows are well known, arising due to steep gradients and phase coupling, see [41,42]. A main requirement here is a robust numerical solver that can handle shocks. Therefore, in this work, a non-oscillatory central (NOC) scheme for conservative laws is employed, derived by Kurganov and Tadmor [43], together with total variation diminishing through a minmodlimiter. The time stepping of the semi-discrete scheme is handled with a two-step Runge-Kutta scheme. As for insights into the numerical implementation for a (simpler) system of partial differential equations, its vectorial form and the discretization, we refer to [4]. The Courant-Friedrichs-Lewy (CFL) condition is set to hold during the computation, determining the time step size via estimates of the wave speeds and the grid size. The associated CFL number is set to = 0.1 during the computations. The local speeds of propagation, i.e., the wave numbers, are left as in [4], that is for a system without the additional fields. Here, they do not depict the actual eigenvalues of the Jacobians of the flux-term, but are rather bounds to the maximum wave speeds that are evaluated via reasonable estimations, see [44,45] for further details. It might appeal to the reader that the actual performance of the numerical code is quite satisfactory, especially in comparison to discrete particle studies. An average computation for studies of the laboratory scale with about 105 (here: 1.2E5) degrees of freedom (DOF), as in Sections 3.2 and 3.3, takes about 5 min, a large scale study, as in Section 3.4, with 106 DOF (here: 8.2E5) from 5 up to 8 h, all on a simple personal computer (intel i5 CPU @3.10 GHz with 4 cores and 8 GB memory, ubuntu OS). 3.2. Parameter studies It is the aim of this first series of investigations to explore the impact of the additional fields of extra pore-fluid pressure and hypoplasticity, and, with this, provide its comparison to a “regular model”, as presented in [10] and transferred to general coordinates in [4]. While results for such a comparison for a 1D case, treated in a system of curvilinear coordinates, have already presented in

[22], we now also investigate contour plots, i.e., the distribution in two dimensions, defined at levels in relation to the increments hst = 0.0125. Moreover, the results presented here shall provide a first overlook, dealing with the general behavior of the mass flow in conjunction with the influence of the additional fields, before in the following sections, special cases and events on a larger scale are tackled. In the following, a combination of Cartesian coordinates (X, Y, Z) and terrain-following coordinates (ξ , η, ζ ) are applied, with respect to the general coordinate ansatz, see Fig. 5 and, for more information, Tai et al. [4]. The computational domain is defined as (X, Y) ∈ ([ 0, 80 ], [ 0, 36 ]), or (ξ , η) ∈ ([ 0, 88.34 ], [ 0, 39.43 ]), respectively. For the parameter studies concerning the role of the extra porefluid pressure and the intergranular stress, the case of a curved, inclined chute was considered, running out into the horizontal plane, see Fig. 5. The transition between the ramp and the horizontal plane is smooth and is described by the inclination angle ϕ , as

⎧ ⎨ϕs , ϕ = ϕs (1 − (X − 24 )/10 ), ⎩ ◦ 0 ,

0 ≤ X ≤ 24, 24 < X < 34, X ≥ 34.

The chute has an inclined part in the range X ≤ 24, a horizontal part in the range for X ≥ 34, and a transitional part in between. The curved character is achieved with an additional elevation zb , such that

⎧ 2 ⎪ ⎨Yd / cos ϕ ,  zb = Y 2 / cos ϕ 1 − (X − 24 )/10 , d ⎪ ⎩ 0,

0 ≤ X ≤ 24, 24 < X < 34, X ≥ 34,

where Yd = |η − η0 |/36 and η0 = 18 (center line). A pile, consisting of a saturated granular-fluid mixture, is released out of rest at the initial time t = 0. This initial mass has a parabolic profile, defined with respect to the terrain-following coordinates as





h(t = 0, ξ , η ) = h0 1 − (ξ − ξ0 )2 /16 − (η − η0 )2 /16 , with h0 = 2/ cos ϕs and ξ0 = 8. Besides the density ratio αρ = ρ f /ρ s , there is a set of four ba-

sic parameters (δb , NR , ϑb , cD ), that govern the flow dynamics of a simple, granular-fluid mixture as in [10], given for these simulations in Table 1. First, there is an angle of basal solid friction δ b , determining the resistance of the granular material at the base. This parameter has a precise physical meaning and can be measured by the inclination angle, at which granular material begins to slide on an inclined plane. The basal friction governs the main drag f

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155 Table 1 Numerical and physical parameters for the studies in Section 3.2. Num. Parameter

Value

Description

(nX , nY ) X Y CFL

ϕs ν s,f |0 vs,X f |0 , vYs, f |0

(161,73) ∈ [ 0, 80 ] ∈ [ 0, 36 ] 0.1 40◦ 0.5 0

Points in X- and Y-direction Domain in X-direction Domain in Y-direction CFL number Maximum chute inclination angle Initial volume fractions Initial velocities

Phys. Parameter

Value

Description

αρ = ρ f /ρ s δb ϑbf

1.06/2.65 23◦ 20 6 200

Density ratio Angle of basal friction Navier fluid friction coefficient Drag coefficient Viscous number

cD NR

on the granular material, so with an increasing angle of basal friction, steeper gradients of height arise, the mixture mass is moving less far and zones of solid concentration develop more distinctly. For determination of values in experimental investigations, we refer to [7,46], where δ b -values between 19.5◦ and 25◦ have been detected. For the fluid, the viscous behavior is adjusted with NR . Neglecting the complexity of the viscous fluid behavior in debris flows, in which the viscosity parameter may depend on the shear-rate, on the presence of clay and silt particles in the pore-water, or on other influences, an idealized case is assumed with a constant and dimensionless viscous number NR = 200. This refers to rather low viscosity of 0.5247P · s, a value fairly close to the range presented in [47]. The corresponding value of μf can be determined with √ NR = (ρ f H gL )/μ f , where L = H = 0.1 and the fluid density is f ρ = 1060kg/m3 . In conjunction with viscosity, the fluid friction f coefficient ϑb determines the basal resistance of the fluid. Here, as well as with the drag coefficient cD for the drag between the two phases, we rely on values given and tested in [10] and also [4], not backed by experimental values. For further discussion of these basic parameters, please see both aforementioned works, as well as [9]. As for the effects of the extra pore-fluid pressure, three different parameter settings have been compared in Fig. 6: A flow without extra pore-fluid pressure and two configurations, in which the extra pore-fluid pressure takes influence with different intensity, via variations of the Euler number Eu. The respective terms act upon the velocities in solid basal friction and in the pressure terms as a driving force. As a main result here, it is apparent that the dynamic extra pressure drives the bulk mass further. When the material is pushed together and compressed during deceleration, the extra pressure is increased, counters settling and reduces the intergranular friction, so the runout length is increased. This matches expectations, since the effects of the dynamic pore-fluid pressure are particularly present during the last, decelerating phase of the flow, when the material reaches the horizontal runout. As shown in [22], an increased initial value of the pore pressure brings this effect to the initial phase, representing a flow onset due to liquef faction, e.g. after intense rainfall. Here, with e |0 = 0.001, a rather small value is chosen, implying little dynamic pressure at the initiation. The parameters for the pore pressure evolution are chosen with reference to [17], with a permeability constant of k0 = 1.0 · 10−13 and the constant of compressibility aD = 1.0, as well as a critical solid volume fraction νCs = 0.62. The average particle size is chosen as δP = 0.0 0 01 and the solid density is ρ s =2650kg/m3 . Further model parameters in conjunction with dilatancy are κ 1 = 1.1 and κ 2 = 15. A possible range for these values is given in [17] via √ the bounds 10−10 ≤ (km L/g)/(αD · μ f · H2 ) ≤ 104 , depending on

145

√ the scale of the event, while we have (k0 L/g)/(aD · μ f · H2 ) ≈ 10−10 . For further explanation of these parameters and their range, see [17], and for the application to the current model, [22]. The f initial pore-fluid pressure is considered as e |0 = 0.001 Further, the influence of hypoplasticity was evaluated in Fig. 7. It provides an inner stabilization, preventing the bulk mass from dissolving and enhances the concentration of solid material. So due to the introduction of the intergranular stress, the material piles up with a steeper front. This is in particular apparent if the solid volume fraction is observed, see Fig. 7(c). Since the solid material stays in an assembled form, with a rather peaked than copped top, the fluid tends to leak out at the front during settling at later time steps, here of course depending on the other material parameters f such as ϑb and cD . This depicts a situation in which the granular skeleton is stable after the bulk mass has decelerated, while the fluid runs out. The influence of the hypoplastic stress can be varied together with the dimensionless number NZ , giving the relation of the intergranular stress to inertial forces. The value of the initial intergranular stress is chosen as Ziis |0 = 0.1, and, as further parameters, an internal friction angle φint = 23◦ is applied, see [7,46] for experimental determination, together with the coefficients fD0 = 1.0 and fs0 = 0.01, as well as a0Z = 0.01. Literature values can be found in [32,39] and there are numerous works on the calibration of the model parameters. A possibly present draw-back is the appearance of unrealistically high values of the solid volume fraction, see [22], since in single-layer models, the maximum packing of the granular particles may be exceeded. As pointed out there and in [9], this calls for a multi-layer approach as in [45], in which real phase separation can be depicted. Nonetheless, the introduction of an actual field in order to depict the contact forces of the granular skeleton provides a further approximation to realistic debris flow behavior. As for a study of the combined influence, both of extra porefluid pressure and hypoplastic stress, we refer to the Appendix Appendix B. There, two more configurations are presented and discussed to display the full model. 3.3. Comparison with experiment Now to test the model with data, a comparison with an experiment was performed. The setup is a classic dam break scenario, conducted by Y.-C. Tai at National Cheng Kung University (NCKU), Taiwan. The results have not been published so far. In this experiment, a block of glass beads, saturated with water, is released from rest and flows to one side while on the other side there is a wall supporting the material. There are two transparent smooth sidewalls, so the configuration is approximating a quasi-2D case, with the (X, Z)-plane being observed. Fig. 8 shows the course of the experiment, captured in a photo, together with a sketch of the experimental setup on the right-hand side (white lines). The flow front is investigated by taking pictures from the side at two time steps after the initiation, t=0.2s and t = 0.4s. The main aim of this comparison is to capture the influence of hypoplasticity on the flow front. For further dam break experiments and numerical investigation, we refer to [48–50]. Note that the role of the initial packing, as examined in [51], has not been examined here. In Fig. 9, the experimental results are compared to numerical results with and without the additional effects at two time steps, i.e., with the parameter configurations Eu = 0.0, Eu|b = 0.0, NZ = 0.0 and Eu = 0.00125, Eu|b = 0.125, NZ = 1.5. For further parameter values, please see Table 2 in Appendix C; values not appearing here are set as before. Panel (a) shows the developing flow front at t=0.2s, together with the initial profile (grey dashed line). It is apparent that the results for the numerical simulation with the additional effects are roughly compatible to the

146

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

Fig. 6. Study on the extra pore-fluid pressure: The height profile is compared at t = 10 and t = 14, together with the volume fraction, for a varying influence of the extra pore-fluid pressure, denoted by changes in Eub . In panel (a), the dashed line denotes no pore pressure, the solid line gives the results for Eub = 0.05 and the dotted line plots the results for Eub = 0.1. Note that, while in panel (a), a cutline at Y = 1.8 (center line) shows the height development, panel (b) shows a height contour plot (scaled with increments 0.25 · hst for t = 10 and 0.1 · hst for t = 14), and panel (c) gives the contour of the solid volume fraction (increments 0.1 · hst for t = 10 and 0.05 · hst for t = 14). For the contour plots, only two configurations are presented, the results for Eub = 0.0 (black line) and Eub = 0.1 (gray line).

experimental results, maintaining a rear part with the initial height and, followed by a rather steep decrease near the front while, without the additional effects, the whole bulk mass fluidizes. The same different behavior is visible at the next time step, t = 0.4s, in panel (b) where a cliff-like socket near the front can be depicted only with the additional fields. This means, the intergranular friction provides a mechanism that hinders the material from dissolving like a fluid, preserving its structure. In contrast, the simulation results without hypolasticity show a rather smoothed front without remnants of the originally sharp edge of the initial block. To provide further qualitative insights, the flow front at t=0.4s is depicted more closely in Fig. 10, together with the height of the

solid phase, hν s . This zoom-in compares the distribution of the solid material with and without an influence of the intergranular stress. This emphasizes the role of the intergranular stress, since it enables the material to reproduce a state in which there is a more distinct, sharp solid cliff together with a water flowout. It can be seen that the steep flow front apparent in the experimental results can be reproduced well by the numerical results by taking hypoplasticity into consideration. All of this highlights the possibility of describing the granular material also as a solid-like structure that maintains its form and refrains from dissolving like a fluid. The apparent limitation of this comparison is due to the application of a system of depth-integrated equations for an actual 2D

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

147

Fig. 7. Height development and volume fraction distribution for a study on the influence of hypoplasticity, given at two different times, t = 10 and t = 14 in the two rows, respectively. The influence is investigated by varying NZ , giving three different configurations in panel (a), denoted by a dashed line for NZ = 0.0, a solid line for NZ = 0.5 and a dotted line for NZ = 0.8. While panel (a) presents a plot at Y = 18 (center line), panel (b) gives the height contour (scaled with increments 0.25 · hst for t = 10 and 0.14 · hst for t = 14) and panel (c) the volume fraction contour plot (with increments 0.1 · hst for both t = 10 and t = 14). For the contour plots, only two configurations are presented, the results for NZ = 0.0 (black line) and NZ = 0.8 (gray line).

problem. Nonetheless the results are usable and satisfying, showing that distinct features in height development are roughly captured by incorporating the additional effects. 3.4. Large scale application: Hsiaolin landslide In August 2009, typhoon Morakot hit Taiwan, causing severe flooding and setting new records with respect to precipitation. Accompanied by heavy rainfall, it triggered a series of catastrophic landslides and mudflows. In the mountains of southern Taiwan, the village of Hsiaolin was destroyed by a debris avalanche, resulting in 474 casualties, see [52]. The events were particularly devastating

since the transported debris reached the Qishanxi river and formed a temporary dam there that, when it broke within an hour, flooded the previously spared parts of the village, see [3,53]. Afterwards, digital terrain models (DTM) were created in order to retrace the course of events and with this, it was estimated that the major body of the landslide had an extent of 57 · 104 m2 and a volume of roughly 24 · 106 m3 , see [52]. For the geomorphological features and a more detailed description of the terrain, also see [53]. In the following, we seek to retrace the disastrous event, and present a brief comparison of the results with respect to the additionally incorporated fields, presenting a further aspect of the practical use of the new model. It is especially interesting to inves-

148

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

Fig. 8. Background and setup of the dam break experiment. The photo depicts the height development of the mass (at t = 0.2s) as it was evaluated during the experiment. On the right-hand side, in the top corner, a sketch shows the abstract setup, including the initial block of glass beads and water (dashed white line) that is supported on the left side by a wall and released at t = 0, as well as the developing flow height (solid white line) at a later time.

Fig. 9. Dam break Experiment: Comparison of the observable height development at two different time steps with numerical simulations. For the numerical results, two different configurations are applied to investigate on the impact of the additional fields with respect to the description of the height development. The black dashed line depicts the experimental results for the height, the gray dashed line gives the initial height profile, the dotted black line gives the configuration without additional fields, Eub = 0.0 and NZ = 0.0, and the black solid line depicts the results for Eub = 0.125 and NZ = 1.5. Note that, while the results shown in the photograph in Fig. 8 display an equal axis scaling, the aspect ratio in the plot is roughly 1: 3.

Fig. 10. Dam break experiment: Zoom-in at the flow front, for t=0.4s. Besides the height of the mixture for the experiment (black dashed line), the height from the numerical studies are given (black solid and dotted lines) and the height of the solid part is depicted (red solid and dotted lines). The dotted lines account for Eub = 0.0, NZ = 0.0 and the solid lines for Eub = 0.125, NZ = 1.5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

tigate the final extents of the height, since here, a comparison to gathered DTM data and aerial photography is possible, see Fig. 12c and d. Besides, a discussion of the volume fraction distribution may provide insights that go beyond the findings from the DTM data. Nonetheless, it is not our aim here to present an extensive reassessment of the events, since this would demand a separate work of its own. A broader investigation together with discussion of parameters and results has been conducted in [4]. f The physical parameters are set to (δb , ϑb , cD , NR , αρ ) = ◦ (16 , 8, 6, 268, 1.06/2.65 ). These chosen basic parameters are discussed in more detail in [4]. It should be noted that they are not the result of a process of optimized fitting but resemble reasonable values and are rather picked to show the applicability of the model to a large scale event. As for the numerical setting, the grid is given with (nX , nY ) = (371, 222 ) on a domain of X ∈ [ 0, 3700m ], Y ∈ [ 0, 2210m ], and CF L = 0.1. The non-dimensional numbers for the extra pore-fluid pressure and the hypoplastic stress are set to Eu|b = 0.075, NZ = 0.035. With respect to the scale of the event, a larger particle diameter value is chosen, δP = .01, together with a slightly changed critical volume fraction νCs = 0.65. Furthermore, for the hypoplastic equations, the initial condition is changed to Ziis |0 = 0.0, accounting for a fully liquefied state in which the inter-

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

149

Fig. 11. Numerical simulation of the Hsiaolin landslide: Height development and deposition depth at 4 different times (panels (a)-(d)), as well as volume fraction development for two time steps (panels (e)-(f)). The basal friction angle is set to δb = 16◦ , the fluid friction coefficient ϑbf = 5, the drag coefficient cD = 6 and the viscous number NR = 268.

granular contact forces are ignored. Besides, the values are left as before. The height development is shown in panels (a)-(d) of Fig. 11, together with the distribution of the solid volume fraction at two different times in panels(e)-(f). With this, it can be retraced how the flow advances, splits in different branches, and settles in the valley of the river, forming three deposits. To allow for an evaluation of the results, Fig. 12 depicts a comparison of the deposition at the last time step, t = 181.83s, both for a configuration with the additional fields (Eu|b = 0.075, NZ = 0.035, panel (a)) in Fig. 12 and without (Eu|b = 0.0, NZ = 0.0, in panel (b)). Furthermore, the results are confronted with an aerial photograph (panel (c)), and some features are added in the figure with the final deposition height (panel (a)): the village Hsiaolin is marked by a thick red line, together with the briefly formated dam

and the flow direction of the river. The formation of the temporary dam can be retraced, while the subsequent failure and the flooding of the remaining parts is not part of the simulation, since here, the river acting upon the dam is not considered. For more details on the extents of this dam, see [3,54]. While the overall results for the height distribution are not differing massively for both studies (with and without additional fields), it is still apparent that due to the influence of the dynamic pore pressure, more material is driven into the deposition area. With this, the results of the DTM study are met, see panel (d) in Fig. 12, without artificially reducing parameters of solid or fluid friction. Please note that in our simulations, no erosion is depicted, so only the green and blue colors in panel (d) resemble deposition also visible in our results. Even more, it can be reasoned that with the incorporation of an intergranular stress, the apparent formation of the (temporary) dam can

150

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

Fig. 12. Comparison of the results for the Hsiaolin landslide: In panel (a), the final deposition height is displayed and the extents of the Hsiaolin village, the flow direction of the river and the briefly formated dam have been included, together with the results without the impact of additional fields in panel (b). An aerial photograph in panel (c) shows the extents of the event, followed by the results of a DTM study in panel (d).

not only to investigate on a laboratory scale but to look into events on real terrain. It is apparent that the results match qualitatively the investigations that have been conducted after the events, see [53]. It also should be noted that with a large scale study like this, small parameter alterations may have a greater effect. Due to the incorporation of the extra pore fluid pressure, a realistic bed friction angle δb = 16◦ can be applied, so instead of mimicking effects with parameter adjustment (as it is sometimes suggested in the literature), the actual physical mechanisms provide for accurate results. Even more, the hypoplastic stress supports the building of a temporary dam. As for the differences to a study without these additional effects, we refer to [4], where the respective results have been presented in greater detail. Fig. 13. Comparison of the solid volume fraction distribution for t = 181.83s, highlighting the influence of both the dynamic pore-fluid pressure and the intergranular stress. A flow configuration with the additional effects (panel (a)) is compared to the distribution without the influence of extra pore-fluid pressure and without hypoplasticity (panel (b)), showing that more solid material is transported further.

be expected in a favorable way, since a solid body formates more distinctly, visible in the increased share of the solid volume fraction at deposition areas at the flow front, see Fig. 13. In this section, the applicability to large scale events, with GIS/DTM data or alike, has been demonstrated, making it possible

4. Conclusion In this work, a two-phase model for debris flows, incorporating pore-pressure evolution and hypoplastic intergranular stress, has been applied to a rugged, general topography and tested extensively with parameter studies, comparison to experimental results and a study of the Hsiaolin event. Taking into account the additional fields of extra pressure and stress state, interaction mechanisms and complex material behavior are captured in greater detail. With this, it is the aim to come closer to a realistic modeling of the mechanisms that govern the dynamics of debris flows. Basi-

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

cally, the fluid phase is treated as a pore-fluid, so with the evolution of the pore pressure, some of the structural information that are absent in classical mixture theory are regained, while the ability of the computational treatment of large flows is kept. A new contribution, both as the continuation of previous work by the authors and in comparison to other approaches in this field, is in particular the combined introduction of pore pressure evolution and hypoplasticity for depth-integrated flows on a general topography. The numerical implementation of this model allows for the efficient numerical treatment of problems, taking about 5 min computational time on a simple personal computer for a small scale laboratory experiment and 5–8 h for large scale natural hazard in mountainous terrain. Several numerical studies are conducted. Starting with the investigation of the additional fields for a simple, curved slope, it is apparent that the extra pore-fluid pressure increases the mobility of the flow, pushing the material further, as it hinders the material from settling. In contrast, the intergranular stress provides a stabilizing effect, preventing the granular structure from dissolving fluid like. Both results are in accordance with the expectations and previous studies, see [22]. The model has also been tested in comparison with the experimental results of a dam break study. Here, the additional fields of the dynamic pore pressure and the intergranular stress also show their justification, providing qualitatively more realistic results. Furthermore, a study of the Hsiaolin landslide has been conducted, to test the model with a real case event. This exhibits the possible application to real-scale events with good agreement between our computed results and gathered data of the real events. As a general advantage here, a most general model is derived in a consistent process. It includes the description of internal mechanisms by additional fields, thus depicting the physical processes in greater detail, compared to classical SH models. Of course, this brings the efforts of additional calibration. However, as has been pointed out, we follow an approach of greater generality and by this, more realistic description of the process. This enables the investigation not only of height deposition, but also of internal mechanisms, the pressure evolution and the stress state. Without this, unrealistically adjusted parameters may hide those internal mechanisms, like, for example, if it is argued that the basal frictional angle should be lowered due to the apparent pore pressure. It can be argued that, without the additional effects, one may still achieve similar results by respective parameter fitting, but mechanisms behind the dynamics are rather concealed. Taken together, the model and its numerical implementation provides some kind of tool box, open for further investigations and, with the additional fields, also the possibility of additional modeling work on the parameters. From here, the derivation of simpler models is possible, making assumptions and constrictions on generality transparent. However, the model remains open to further amendments. Processes of erosion and deposition, i.e., a non-material basal surface, are an important field of research, as well as effects of segregation. Even more, the interdependency of such effects with both the pore pressure and the intergranular friction would leave room for extensive amendments and improvements – as would the mutual interference of pore pressure and hypoplastic stress in the framework of this model as well. Also, as a point of possible connection, there are further approaches that apply a pressure dependence to the rheology and thus may provide ideas for such interdependencies, see for example [55] and [56], as well as [49], integrating the Drucker-Prager plasticity in the framework of the μ(I)-rheology. Further determination of parameters and calibration is worth looking into. As stated above, the presented approach rather offers the basis for further modeling work than delivering an ultimate instrument.

151

Acknowledgement The authors want to acknowledge the financial support of the Deutsche Forschungsgemeinschaft (DFG) for this work within the project number WA 2610/3-1, and Y.C. Tai would like to thank the Ministry of Science and Technology, Taiwan, for the financial support (Project: MOST 107-2221-E-006-032).

Appendix A. Derivation of pore pressure and hypoplasticity equations Since the full derivational process is quite lengthy, we rather point at the way in which the equations can be transformed in terrain-following coordinates. For the balance of (solid) mass, the transfer of the equation in curvilinear coordinates

∂ ∂ ∂  s f hν vy = 0, (hν s ) + (hν s vsx ) + ∂t ∂x ∂y

(7)

to final equation yields

 ∂   ∂ ∂  J hν s vsξ + J hν s vηf = 0, (J hν s ) + ∂t b ∂ξ b ∂η b

(8)

that is easy to follow. The introduced Jb denotes the determinate of the basal Jacobian. Also note that c = nz is the z component of the unit normal vector n = (nx , ny , c ) of the topographic surface S. The equation for ν f is derived accordingly. For the momentum, 4 balances arise, with X- and Y-direction and solid/fluid constituent.  s and the respective fluid We introduce the solid pressure term   f as pressure term 

  ρ f s Jb h2 c J νs f h s = 1 −  − Eu b s e , ν s ρ 2 ρ

f = νf 

h2 c ν f ef h + Eu . 2 ρf

(9)

For the centrifugal forces terms κ s , κ f we refer to [4]. The solid, X momentum equation states

 ∂   ∂ ∂  J hν s vsξ vsx + J hν s vsη vsx (J hν s vsx ) + ∂t b ∂ξ b ∂η b ∂ ∂  s A11 ) −   s A21 ) = Jb hν s (c −  μ κ s )nx −  ( ( ∂ξ ∂η       (i) gravitation  cD ν ν Jb h f vx − vsx + ρs    s

(ii) pressure gradients

f

(iii) drag

     s s ∂ Jb hν s Zξs ξ + Zηη ∂ Jb hν s Zξs ξ + Zηη ρf + NZ 1 − s A11 + A21 ρ ∂ξ ∂η   s ∂ Jb ρ hZxξ ∂ Jb ρ hZxsη  + s NZ + ρ ∂ξ ∂η   f 2 ρ s ∂ Jb h c ∂ Jb h2 c − ν A11 + A21 2ρ s ∂ξ ∂η     f ν f ρs − ρ f 2 h cJb ∂ν s ∂ν s ρ + s A11 + A21 ρ ρ 2 ∂ξ ∂η    f ρ f − sgn(vsx ) tan(δb )hν s c 1 − s − Eu|b es −  μ κ s , (10) ρ ρ    

(iv) basal friction

and the Y momentum equation is given as

152

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

 ∂   ∂   ∂ J hν s vsy + J hν s vsξ vsy + J hν s vsη vsy ∂t b ∂ξ b ∂η b ∂ ∂  s A12 ) −   s A22 ) = Jb hν s (c −  μ κ s )ny −  ( ( ∂ξ ∂η  cD ν s ν f Jb h  f vy − vsy + ρs       s s ∂ Jb hν s Zξs ξ + Zηη ∂ Jb hν s Zξs ξ + Zηη ρf − NZ 1 − s A12 + A22 ρ ∂ξ ∂η   s ∂ Jb ρ hZyξ ∂ Jb ρ hZys η  + s NZ + ρ ∂ξ ∂η   ∂ Jb h2 c ∂ Jb h2 c ρf s − ν A12 + A22 2ρ s ∂ξ ∂η     f s f ∂ν s ∂ν s ρ f ν ρ − ρ h2 cJb + s A12 + A22 ρ ρ 2 ∂ξ ∂η    f   ρ f −sgn vsy tan(δb )hν s c 1 − s − Eu|b es −  μ κ s , (11) ρ ρ where, please, be reminded that s μ 12 s Zxsξ = 11 b Zξ ξ +  b Zηξ ,

s 12 s Zxsη =  μ 11 b Zξ η + b Zηη .

(12)

s 22 s Zys η =  μ 21 b Zξ η + b Zηη .

(13)

and s μ 22 s Zys ξ = 21 b Zξ ξ +  b Zηξ ,

ij

The components b of the basal Jacobian transformation matrix are applied, as well as the components of the inverse transformation matrix, Aij . Besides the driving gravitational term (i), the hindering basal pressure (iv) and the pressure gradients (ii), the solid momentum equations are comprised by the momentum interaction between both phases via drag (iii) and buoyancy (lines 5–6), as well as the influence of the intergranular friction (lines 3–4). For the fluid, the respective X momentum balance equation yields

∂   ∂   ∂ Jb hν f vxf + Jb hν f vξf vxf + J hν f vηf vxf ∂t ∂ξ ∂η b   ∂  f  ∂  f   A11 −  Jb  A21 = Jb hν f c −  μ κ f nx −  Jb   ∂ξ ∂η       (i) gravitation

(ii) pressure gradients

 cD ν ν Jb h  f − vx − vsx f ρ    s

f

∂   ∂   ∂ Jb hν f vyf + Jb hν f vξf vyf + J hν f vηf vyf ∂t ∂ξ ∂η b   ∂  f  ∂  f   A12 −  Jb  A22 = Jb hν f c −  μ κ f ny −  Jb   ∂ξ ∂η  cD ν s ν f Jb h  f − vy − vsy ρf   ϑ f ν f Jb ∂ Jb vηf ∂ Jb vηf  ∂ + hA12 + hA22 − b f vyf (15) NR ∂η ∂ξ ∂η ρ   ∂ Jb vξf ∂ Jb vξf ∂ Jb vηf ∂ Jb vηf  1+μ ∂ + h A12 + A22 + A11 + hA21 NR ∂ξ ∂ξ ∂η ∂ξ ∂η   h2 cJb ∂ν s ∂ν s − A12 + A22 2 ∂ξ ∂η  s    f f 2 ν ρ − ρ h cJb ∂ν s ∂ν s − A12 + A22 . (15) ρ 2 ∂ξ ∂η The fluid momentum equations comprise, additionally to the already introduced terms (i)-(iii), Navier basal friction (iv), buoyancy (line 5) and the viscous terms in conjunction with NR , see lines 3 and 4. The evolution equation for the extra pore-fluid pressure is given as

∂   ∂   ∂ km β f Jb hef + Jb hef vξf + Jb hef vηf − 2 Jb ef ∂t ∂ξ ∂η hμ f αD    f f s s ∂ Jb vξf ∂ Jb vηf σe ∂ Jb hν (vξ − vξ ) ∂ Jb hν s (vsη − vη ) = + +h + 2 ∂ξ ∂η ∂ξ ∂η   ∂ Jb hν s (vsξ − vξf ) ∂ Jb hν s (vsη − vηf ) +ef + − Jb hγ˙ tan(ψ ). (16) ∂ξ ∂η In the following the inner stress is defined with respect to the tangential vector space (not its reciprocal bases), i.e.,



For intergranular hypoplastic stress, we only state the Xξ equation, since the others follow accordingly and, as above with the momentum equations, for reasons of length. It states:

(iii) drag

  ∂ Jb vξf ∂ Jb vξf  ∂ + hA11 + hA21 − NR ∂ξ ∂ξ ∂η

(iv) Navier friction

 ∂ Jb vξ ∂ Jb vξ ∂ Jb vηf ∂ Jb vηf ∂ + h A12 + A22 + A11 + hA21 NR ∂η ∂ξ ∂η ∂ξ ∂η     ν f ρ s − ρ f h2 cJb h2 cJb ∂ν s ∂ν s − A11 + A21 − 2 ∂ξ ∂η ρ 2   ∂ν s ∂ν s × A11 + A21 . (14) ∂ξ ∂η 

1+μ



ϑbf ν f Jb f v ρf x   

f

f

and the Y momentum equation states



Zi, j = Zi, j (g j  g j ).



∂ Jb hZxsξ ∂ Jb hvsξ Zxsξ ∂ Jb hvsη Zxsξ + + ∂t ∂ξ ∂η   ∂ Jb hvsξ ∂ Jb hvsξ s ∂ Jb hvsη ∂ Jb hvsη + A11 + A21 − A12 − A22 Z xη ∂ξ ∂η ∂ξ ∂η    ρs − ρ f s ∂ J hν s (vsx − vxf ) ∂ J hν s (vsx − vxf ) − 1 + ν s Zxξ A11 b + A21 b ρ ∂ξ ∂η  ∂ Jb hν s (vsy − vyf ) ∂ Jb hν s (vsy − vyf ) +A12 + A22 ∂ξ ∂η −Zxsζ Jb hvsx (1 − λκ ) =  sxξ Jb h.

(17)

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155



The source is given as





+Zys η  h A12

∂ Jb vsξ ∂ Jb vsξ sxξ Jb h = fs  asz h A11 + A21 ∂ξ ∂η 



∂ Jb vsξ ∂ Jb vsξ Zxsξ  h A11 + A21 ∂ξ ∂η

  ∂ Jb vsη ∂ Jb vsη + A22 +Z ∂ξ ∂η  (∗ )





153

+



hZxsξ Zxsξ + Zys η + Zzsζ



∂ Jb vsη + Zxsη  h A11 ∂ξ  ∂ Jb vsξ ∂ Jb vsξ ∂ Jb vsη +A21 + A12 + A22 ∂η ∂ξ ∂η

2

    ∂ Jb vsξ ∂ Jb vsξ  Zxsξ h  + fD asZ h A11 + A21 s  ∂ξ ∂η Zxξ + Zys η + Zzsζ  2 s  1 s 1 s +

h

Z

3 xξ Zxsξ

− 3 Zyη − 3 Zzζ + Zys η + Zzsζ

.

(18)

Appendix B. Parameter studies: combined influence

Fig. 14. Study on the influence of the pore-fluid pressure and hypoplasticity acting together, with three different configurations. Two time slices are depicted, t = 10 and t = 14. The influence is investigated by varying NZ and Eu|b , giving three different configurations in panel (a), denoted by a dashed line for no influence, a solid line for Eu|b = 0.1, NZ = 0.5 and a dotted line for Eu|b = 0.05, NZ = 0.1. While panel (a) presents a plot at Y = 18 (center line), panel (b) gives the height contour (scaled with increments 0.25 · hst for t = 10 and 0.14 · hst for t = 14) and panel (c) the volume fraction contour plot (increments 0.1 · hst for t = 10 and 0.075 · hst for t = 14). For the contour plots, only two configurations are presented, the results for Eu|b = 0.0, NZ = 0.0 (black line) and Eu|b = 0.1, NZ = 0.5 (gray line).

154

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155

Appendix C. Comparison with experimental results: parameters

Table 2 Parameters for dam break simulations in Section 3.3. Parameter

Value

Description

(nX , nY ) X Y CFL αρ = ρ f /ρ s ν s |0 vs,X f |0 , vYs, f |0

(161,41) ∈ [ 0, 8 ] ∈ [ 0, 2 ] 0.1 1.0/2.4 0.6 0 0.001 0.1 0.62 0.75 10−13 1.0 0.01 36◦ 40 12 300 23◦ 1.0 0.01

Points in X- and Y-direction Domain in X-direction Domain in Y-direction CFL number Density ratio Initial solid volume fraction Initial velocities Initial Pore-Fluid Pressure Initial Hypoplastic Stress Critical volume fraction Maximum packing fraction Hydraulic permeability coefficient Compressibility factor Fluid viscosity factor Angle of basal friction Navier fluid friction coefficient Drag coefficient Viscous number Angle of internal friction Density coefficient Stiffness coefficient

ef |0

Ziis |0

νCs s νmax k0 aD

μf δb ϑbf

cD NR

φ int f D0 f s0

Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.compfluid.2019.06.015

References [1] Larsen MC, Wieczorek GF, Eaton L, Torres-Sierra H. Natural hazards on alluvial fans: the debris flow and flash flood disaster of december 1999, Vargas state, Venezuela. In: Proceedings of the sixth caribbean Islands water resources congress, Puerto Rico, 965; 2001. p. 1–7. [2] Wieczorek G, Larsen M, Eaton L, Morgan B, Blair J. Debris-flow and flooding hazards associated with the december 1999 storm in coastal Venezuela and strategies for mitigation. Tech. Rep.; 2001. [3] Dong J-J, Li Y-S, Kuo C-Y, Sung R-T, Li M-H, Lee C-T, et al. The formation and breach of a short-lived landslide dam at Hsiaolin village, Taiwan – part I: post-event reconstruction of dam geometry. Eng Geol 2011;123(1–2):40–59. [4] Tai Y-C, Heß J, Wang Y. Modeling two-phase debris flows with grain- fluid separation over rugged topography: application to the 2009 Hsiaolin event, Taiwan. J Geophys Res 2018;124(2):305–33. [5] Savage SB, Hutter K. The motion of a finite mass of granular material down a rough incline. J Fluid Mech 1989;199:177–215. [6] Pitman EB, Le L. A two-fluid model for avalanche and debris flows. Philos Trans R Soc A 2005;363(1832):1573–601. [7] Savage SB, Hutter K. The dynamics of avalanches of granular materials from initiation to runout. Part I: analysis. Acta Mech 1991;86:201–23. [8] Gray J, Wieland M, Hutter K. Gravity-driven free surface flow of granular avalanches over complex basal topography. In: Proceedings of the royal society of London A: mathematical, physical and engineering sciences, 455; 1999. p. 1841–74. [9] Pudasaini SP. A general two-phase debris flow model. J Geophys Res 2012;117. [10] Meng X, Wang Y. Modelling and numerical simulation of two-phase debris flows. Acta Geotech 2016;11(5):1027–45. [11] Hui W, Li P, Li Z. A unified coordinate system for solving the two-dimensional Euler equations. J Comput Phys 1999;153(2):596–637. [12] Hui WH. A unified coordinates approach to computational fluid dynamics. J Comput Appl Math 2004;163(1):15–28. [13] Bouchut F, Westdickenberg M, et al. Gravity driven shallow water models for arbitrary topography. Commun Math Sci 2004;2(3):359–89. [14] Luca I, Tai Y-C, Kuo C-Y, et al. Shallow geophysical mass flows down arbitrary topography. Springer; 2016. [15] Savage S, Iverson RM. Surge dynamics coupled to pore-pressure evolution in debris flows, 1. Millpress; 2003. p. 503–14. [16] George DL, Iverson RM. A two-phase debris-flow model that includes coupled evolution of volume fractions, granular dilatancy, and pore-fluid pressure. Italian J Eng GeolEnviron 2011;10. 2011–03

[17] Iverson RM, George DL. A depth-averaged debris-flow model that includes the effects of evolving dilatancy. I. Physical basis. In: Proceedings of the royal society A: mathematical, physical and engineering sciences, 470. The Royal Society; 2014. [18] Kolymbas D. A rate-dependent constitutive equation for soils. Mech Res Commun 1977;4(6):367–72. [19] Gray J, Thornton A. A theory for particle size segregation in shallow granular free-surface flows. In: Proceedings of the royal society of London A: mathematical, physical and engineering sciences, 461. The Royal Society; 2005. p. 1447–73. [20] Tai Y-C, Lin Y-C. A focused view of the behavior of granular flows down a confined inclined chute into the horizontal run-out zone. Phys Fluids 2008;20(12):123302. [21] Heß J, Wang Y, Hutter K. Thermodynamically consistent modeling of granular-fluid mixtures incorporating pore pressure evolution and hypoplastic behavior. Continuum Mech Thermodyn 2017;29(1):311–43. [22] Heß J, Wang Y. On the role of pore-fluid pressure evolution and hypoplasticity in debris flows. Eur J Mech-B/Fluids 2019;74:363–79. [23] Müller I. A thermodynamic theory of mixtures of fluids. Arch Ration Mech Anal 1968;28(1):1–39. [24] Liu I-S. Method of lagrange multipliers for exploitation of the entropy principle. Arch Ration Mech Anal 1972;46(2):131–48. [25] Müller I, Liu I-S. Thermodynamics of mixtures of fluids. Rational thermodynamics. Truesdell C, editor. Springer; 1984. [26] Tai Y, Kuo C. Modelling shallow debris flows of the coulomb-mixture type over temporally varying topography. Nat Hazards Earth Syst Sci 2012;12(2):269–80. [27] Terzaghi Kv. The shearing resistance of saturated soils and the angle between the planes of shear. In: Proceedings of the 1st international conference on soil mechanics and foundation engineering, 1; 1936. p. 54–6. [28] de Boer R, Ehlers W. The development of the concept of effective stresses. Acta Mech 1990;83(1):77–92. [29] Iverson RM, Denlinger RP. Flow of variably fluidized granular masses across three-dimensional terrain 1. coulomb mixture theory. J Geophys Res 2001;106:537–52. [30] Goren L, Aharonov E, Sparks D, Toussaint R. Pore pressure evolution in deforming granular material: a general formulation and the infinitely stiff approximation. J Geophys Res 2010;115(B9). [31] Kolymbas D. An outline of hypoplasticity. Arch Appl Mech 1991;61(3):143–51. [32] Bauer E. Calibration of a comprehensive hypoplastic model for granular materials.. Solid Found 1996;36(1):13–26. [33] Wu W, Bauer E, Kolymbas D. Hypoplastic constitutive model with critical state for granular materials. Mech Mater 1996;23(1):45–69. [34] Arnold M, Herle I. Hypoplastic description of the frictional behaviour of contacts. In: Numerical methods in geotechnical engineering; 2006. p. 101–6. [35] Niemunis A. Hypoplasticity vs. elastoplasticity, selected topics. Modern ApproachPlastic 1993:277–307. [36] Herle I, Gudehus G. Determination of parameters of a hypoplastic constitutive model from properties of grain assemblies. Mech Cohesive-Frict Mater 1999;4(5):461–86. [37] Teufel A. Simple flow configurations in hypoplastic abrasive materials. Technische Universität Darmstadt; 2001. [38] Fang C, Wang Y, Hutter K. Shearing flows of a dry granular material-hypoplastic constitutive theory and numerical simulations. Int J Numer Anal Methods Geomech 2006;30(14):1409–37. [39] Gudehus G. A comprehensive constitutive equation for granular materials. SoilsFound 1996;36(1):1–12. [40] Mašín D. A hypoplastic constitutive model for clays. Int J Numer Anal Methods Geomech 2005;29(4):311–36. [41] Wang Y, Hutter K. Comparisons of numerical methods with respect to convectively dominated problems. Int J Numer Methods Fluids 2001;37(6):721–45. [42] Wang Y, Hutter K, Pudasaini SP. The Savage-Hutter theory: a system of partial differential equations for avalanche flows of snow, debris, and mud. ZAMM Zeitschrift für Angewandte Mathematik und Mechanik 2004;84(8):507–27. [43] Kurganov A, Tadmor E. New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations. J Comput Phys 20 0 0;160(1):241–82. [44] Kurganov A, Miller J. Central-upwind scheme for Savage–Hutter type model of submarine landslides and generated tsunami waves. Comput Methods Appl Math 2014;14(2):177–201. [45] Meng X, Wang Y, Wang C, Fischer J-T. Modeling of unsaturated granular flows by a two-layer approach. Acta Geotech 2017;12(3):677–701. [46] Hutter K, Koch T, Plüss C, Savage S. The dynamics of avalanches of granular materials from initiation to runout. Part II. experiments. Acta Mech 1995;109(1–4):127–65. [47] O’Brien JS, Julien PY. Laboratory analysis of mudflow properties. J Hydraul Eng 1988;114(8):877–87. [48] Balmforth N, Kerswell R. Granular collapse in two dimensions. J Fluid Mech 2005;538:399–428. [49] Ionescu IR, Mangeney A, Bouchut F, Roche O. Viscoplastic modeling of granular column collapse with pressure-dependent rheology. J Nonnewton Fluid Mech 2015;219:1–18. [50] Wang C, Wang Y, Peng C, Meng X. Dilatancy and compaction effects on the submerged granular column collapse. Phys Fluids 2017;29(10):103307. [51] Pailha M, Pouliquen O. A two-phase flow description of the initiation of underwater granular avalanches. J Fluid Mech 2009;633:115–35.

J. Heß, Y.-C. Tai and Y. Wang / Computers and Fluids 190 (2019) 139–155 [52] Kuo C, Tai Y, Chen C, Chang K, Siau A, Dong J, et al. The landslide stage of the Hsiaolin catastrophe: simulation and validation. J Geophys Res 2011;116(F4). [53] Lo C-M, Lin M-L, Tang C-L, Hu J-C. A kinematic model of the hsiaolin landslide calibrated to the morphology of the landslide deposit. Eng Geol 2011;123(1–2):22–39. [54] Wu C-H, Chen S-C, Feng Z-Y. Formation, failure, and consequences of the Xiaolin landslide dam, triggered by extreme rainfall from typhoon Morakot, Taiwan. Landslides 2014;11(3):357–67.

155

[55] Gray J, Edwards A. A depth-averaged μ(i)-rheology for shallow granular free– surface flows. J Fluid Mech 2014;755:503. [56] Edwards A, Gray J. Erosion–deposition waves in shallow granular free-surface flows. J Fluid Mech 2015;762:35–67.