Acta Polytechnica CTU Proceedings DOI:10.14311/APP.2020.29.0030 Acta Polytechnica CTU Proceedings 29:30–34, 2020 © Czech Technical University in Prague, 2020 available online at http://ojs.cvut.cz/ojs/index.php/app RECENT ADVANCES IN SIMULATING LANDSLIDE AND DEBRIS FLOW Wei Wu Universität für Bodenkultur Wien, Institut für Geotechnik, Department für Bautechnik und Naturgefahren, Feistmantelstraße 4 , 1180 Wien, Austria correspondence: wei.wu@boku.ac.at Abstract. This paper summarizes the outcome of a research project funded by the European Commission on the simulation of landslide and debris flow. The numerical techniques adopted include the mesh-based method (FEM), the point-based method (SPH) as well as a hybrid method (Lattice-Boltzmann Method and DEM). Centrifuge model tests were carried out for the purpose of validation. Keywords: Debris flow, DEM, FEM, landslide. 1. Introduction Landslides are caused by slope failures. The stability of slopes is a classical topic in soil mechanics. Debris flows are often triggered by landslides either partially or completely mobilizing into debris flows. Globally, landslides cause billions of dollars in damage and thousands of deaths and injuries each year. Methods based on limit equilibriums have been most popular for slope stability analysis, as in the slice method of Bishop. However, the limit equilibrium ap- proach cannot capture the initiation of hydrologically driven slope failure, which is characterized mainly by transient processes. The modelling of infiltration, evaporation, propagation of the wetting front and re- duction of matric suction involves different time scales governing both short-term and long-term behaviour. Moreover, the initiation and propagation of the failure surface cannot be explained by the limit equilibrium approach [1]. A consistent analysis of slope stability requires a fully coupled analysis of the hydrological and mechanical processes in multiphase media, and which considers spatial and temporal multiscale prop- erties. Whereas slope stability has been modelled mainly in the realm of soil mechanics, the modelling of debris flows has been developed mainly within the terms of fluid mechanics, in particular the shallow water model [2]. There are two main distinct approaches to com- putational modelling of debris flows: the continuum and discrete methods. In most continuum models, the flow is considered as a one-phase medium with constant density. Much effort has been devoted to the debris viscosity. Several debris flow runout models have been developed. Contrasting the one-phase ap- proach with constant viscosity are field observations demonstrating that such an approach is unable to de- scribe unsteady (spatial and temporal) flow behaviour properly. However, some important phenomena, such as segregation of coarse particles and the transition be- tween solid-like to fluid-like behaviour are still elusive to such models. Consistent continuum modelling of debris flows requires an advanced constitutive model capable of describing the local inhomogeneity and solid-like/fluid-like behaviour and a robust numerical scheme, e.g. SPH (Smooth Particle Hydrodynam- ics), MPM (Material Particle Method) and PFEM (Particle Finite Element Method). The project MUMOLADE (Multiscale Modelling of Landslides and Debris Flows) was funded within the framework of the Marie Curie ITN (Initial Training Networks) by the European Commission in order to address needs in research and training in the numerical simulation of landslides and debris flows. The funding duration was three years between 2012 and 2015 and involved a budget of about 4.3 xmillion euros. This paper gives a brief account of the results achieved in MUMOLADE. 2. Strain localization and triggering of slope failure by FEM Soil often fails through localized deformation in the form of shear bands, which give rise to progressive fail- ure in slopes with over-consolidated soil. The precur- sor of failure is the localisation of deformation within thin shear bands. The bifurcation theory and the re- lated numerical schemes provide a rational approach to this multiscale phenomenon. The multiphase FEM analysis provides a realistic description of the fail- ure process from initiation and aftermath both under drained and partially undrained conditions. Once a shear band is formed, the underlying partial differential equation (PDE) becomes rather poorly expressed and requires regularization. The elliptical PDF for quasi-static loading will change to parabolic. This change brings about some difficulties for the nu- merical calculations using the Finite Element Method (FEM). The numerical solutions show a spurious de- pendence on the mesh alignment and mesh size. More- 30 http://dx.doi.org/10.14311/APP.2020.29.0030 http://ojs.cvut.cz/ojs/index.php/app vol. 29/2020 Recent advances in simulating landslide and debris flow Figure 1. Equivalent viscoplastic strain profiles (a) and contours for (b) coarse mesh (5 × 17), (c) intermediate mesh (10 × 34) and (d) fine mesh (20 × 68) in the case viscosity η = 30 s. [3] over, the thickness of the shear bands and the spacing between the shear bands require a theory with a char- acteristic length. It should be noted that conventional continuum mechanics does not possess a length scale. There are several techniques to regularize the PDF, e.g. nonlocal theory, micropolar theory, strain gradi- ent theory and viscosity. In most of these theories a characteristic length is introduced. The team led by Prof. Sanavia at Padova, Italy studied the regularization by nonlocal theory. The nonlocal theory assumes that the response of the ma- terial at some specific point is determined not only by the state at that point but also by the state of its neighbouring points. The variables are then substi- tuted with their non-local counterparts obtained by weighted averaging over a spatial neighbourhood of each point under consideration. Obviously, the size of the neighbourhood introduces a characteristic length scale. The following examples shows a plane strain com- pression test on water saturated sand. An elastic- viscoplastic constitutive model is used. The finite ele- ment calculations are carried out using three meshes with an 85, 340 and 1340 element respectively. “The finite element mesh consists of eight-node quadrilat- eral isoparametric elements with a reduced Gaussian (2 × 2) integration scheme to eliminate shear locking. The bottom of the sample was assumed to be fixed and rough whereas the boundaries are impervious and adi- abatic. Quasi-static loading conditions are assumed, and gravity acceleration is taken into account. Axial compression is applied to the specimen by imposing vertical velocity (1.2mm/s) on the top [3]. As can be seen from Fig. 1, the thickness of the shear bands remains fairly independent of the mesh size. This gives us some confidence when using the FE code that the problem of slope stability, where shear band propagation and progressive failure are concerned, can indeed be solved. 3. Stability of vegetated slope by FEM Most shallow landslides are triggered by rainfall. Veg- etation is known to play an important role in reduc- ing landslide hazard. We carried out FE simulations to assess the impact of rainfall and vegetation on slope stability. The vegetated slope is considered as a combined mechanical–hydrological system. A sophis- ticated elastoplastic constitutive model for partially saturated soil was used to describe the mechanical behaviour of the soil. The reinforcement by plant roots is regarded as a hardening mechanism in the constitutive model. The hydrological effect of plant by root uptake is considered in the underlying Richards equation. The architecture of plant roots is also taken into consideration. Fig. 2 shows the displacements of a slope with and without vegetation after a rainfall of up to 120 hours with an intensity of 18mm/h. The slope without vegetation loses stability before the vegetated slope by about 10 hours. The slope has an inclination of 45◦ and a height of 10metres. The root architecture is assumed to be half spherical with a diameter of 1 m. The numerical model has some advantages over the traditional slice method based on limit state analy- sis: “The model is capable of capturing and coupling the main features that should be considered when considering the stability of vegetated slopes. The me- chanical root reinforcement is considered based on a macroscopic approach, which is a significant advantage when applying it to large-scale problems. The root architecture can be approximated with geometrical representations. The soil parameters within the root zones are enhanced. However, the ‘constant additional cohesion’ assumption is avoided, one which is widely used in many existing models. The soil strength is enhanced by the plant roots, whose strength is not constant, but hardens and softens depending on strain. Importantly, the numerical model can simulate the 31 Wei Wu Acta Polytechnica CTU Proceedings Figure 2. Contours of displacements, obtained at the end of calculations, for the two extreme cases: (a) non-vegetated slope; (b) vegetated slope with both mechanical root reinforcement and evapotranspiration considered. Different shapes of failure surfaces can be observed. Note that the circles are a projection of root zones (values of displacements are given in metres) [4] root water uptake of plant roots, which leads to the desaturation of the soil mass and gives rise to further enhancement of soil strength [4]. 4. Discrete Element Simulation Following the work by Cundall, the Discrete Element Method (DEM) has become quite commonly used for granular materials. Debris flows involve a large number of granules. Therefore, a working programme led by Prof. Nicot at INRAE, Grenoble, France is devoted to the numerical simulation of flexible barriers against debris flow by DEM. Check dams are often used as protection struc- tures against debris flow. Recently, flexible barriers (steel meshes) are widely used as protection structures against debris flow in mountainous regions. Such flexible structures absorb the impact energy better than the rigid check dams. The complex interaction between debris flow and flexible structures presents some challenge to the continuum simulation. The DEM offers an attractive alternative to the contin- uum approach. Moreover, the DEM offers possible links between the microscopic properties of granular materials and their macroscopic behaviour. The grains of debris flow are usually modelled by spheres. The normal and tangential forces of the particles are obtained from interactions with the help of the normal and tangential aspects of stiffness. The tangential force can be limited by some failure criteria, e.g. Mohr-Coulomb. Usually the model parameters are calibrated by considering some simple and well- defined problems. Fig. 3 shows a snapshot of the interaction between the debris flow and the flexible barrier. The latter is modelled by net elements. The calculations provide also the forces in the flexible barrier and in the anchor elements. The numerical calculations are carried out with the DEM code YADE developed by colleagues in Grenoble. Figure 3. Snapshots of the granular flow impact on the flexible barrier [5] Figure 4. Set-up of numerical simulation (a) Geome- try of the configuration before the simulation starts. (b) Undeformed shape of the barrier. The snapshots on the right show the geometry of the flow before (c), during (d), and after (e) the impact on the barrier [6] 5. A Hybrid Method combining LBM with DEM Although DEM may provide some insight into the complex interaction between debris flow and the flex- ible barrier, large scale problems require continuum modelling. The group led by Prof. Herrmann at ETH, Switzerland works on the numerical simulation with the hybrid method by combining LBM (Lattice Boltz- mann Method) and DEM. The LBM is well established in computational Fluid Dynamics (CFM). The fluid (water together with fines) is simulated by LBM while the large particles are modelled by DEM. The fluid is 32 vol. 29/2020 Recent advances in simulating landslide and debris flow considered as a viscous fluid with yield limit by the Bingham model. The interaction between the large particles and the surrounding fluid is considered by an interaction term. The flexible net is modelled by finite elements. This hybrid method benefits from the strength of both methods and offers an optimal way to model the complex interaction between debris flow and a flexible barrier. Fig. 4 shows the set-up of the numerical simulation. A fluid saturated granular mass flows down a chute with an inclination of 15◦. The seepage of fluid and the holding back of large particles upon impact with the flexible net can be clearly observed in Fig. 3 (e). Figure 5. Development of failure surface depicted by total equivalent strain with reduction factor f = 1 : 9 [7] 6. Smooth Particle Hydrodynamics For problems with large deformation the FEM often suffers from accuracy and convergence problems as a result of skewed mesh. In recent years, various meshless methods have been proposed to remedy this flaw. Among these, the first and most widely used meshless method is the SPH. Originally proposed to handle hydrodynamic problems in astrophysics, SPH has become very common in fluid dynamics and more recently in solid mechanics. SPH is a pure Lagrangian method, which makes use of fictive particles with assigned properties like mass, density, velocity, stress and strain. The approximation of the physical properties of the particles is based on an integration interpolation over a support domain. The size of this domain defines the interaction among the particles within its neighbourhood. Through an updated Lagrangian scheme, SPH can handle large deformation and complex material interface with some ease. Like FEM, a constitutive model for the granular material is still required. Until now, slope failure and debris flow are usually treated in different domains with different methods. SPH offers a promising approach to handle both prob- lems within one single consistent framework and with one numerical method. Fig. 5. shows the development of the failure surface in a slope over the time. The slope is 5m high and 10m wide. A hypoplastic constitutive model was used for the soil. Some 9000 particles were used for the simulation. The progressive failure process is well captured by the SPH. In the SPH calculations, the thickness of shear band is mainly dictated by the particle spacing and size of the support domain. The thickness of the shear band is known to depend on the mean grain diameter of soil. This needs to be improved. 6.1. Centrifuge tests on a scaled slope model There are only a few cases of landslides with well doc- umented data and instrumentation in the real world. Sophisticated numerical models require more data for calibration and validation. Geotechnical centrifuge of- fers a cost-effective alternative to study slope stability under well-defined initial and boundary conditions on a scaled model. In the scaled model, the size of the prototype slope is reduced by a factor N to maintain the geometrical similarity. The mechanical similarity is maintained by using N-fold of the earth accelera- tion. Centrifuge tests on a scaled slope model are carried out to study the stability of soil slopes under rainfall. The setup and the technical specifications of the geotechnical centrifuge at BOKU can be seen from Fig. 6 and Table 1. A fine sand with d50 = 0.232 mm with fines was used. The grading was as follows: 83% sand, 13.4% silt and 2.8% clay. We carried out triaxial tests on unsaturated soil. The friction angle is about 35◦ but depends on the initial density and the confining pressure. The cohesion depends mainly on the suction i.e. water content. The slope height is 20 cm and the slope inclination is 45◦. The slope was prepared by moist tamping. The displacement was measured through a glass win- dow and evaluated by PIV. Rainfall was simulated by purpose designed nozzles. Fig. 7 shows the shear strain from PIV measurements. The acceleration is 30g, which corresponds to a prototype slope of about 6 m. The test was carried out with initial dry density 1.4 g/cm3, the initial water content 11.5 %, and the rainfall intensity 11.8 mm/h. Some 20 secs after rainfall is initiated, a shear band can be observed. This shear band propagates and becomes highly localized at t = 30 s. By varying the parameters, we can establish a relationship between rainfall intensity and duration. This compares well with the available data in the literature. Our dataset can be used for validation of the numerical simulations. 7. Conclusion Landslide and debris flow remain an active research area. The research is highly rewarding for actual prac- tice but also highly challenging. Better understanding of soil behaviour and better numerical techniques are needed. In order to model slope failure and debris flow within one consistent framework, we need a con- stitutive model for solid-like and fluid-like behaviour. 33 Wei Wu Acta Polytechnica CTU Proceedings Figure 6. Geotechnical centrifuge at BOKU Diameter of Centrifuge [mm] 3, 000 Radius to base of swinging platform, rmax [mm] 1, 308.5 Radial acceleration [g] 0 to 200 Angular velocity [rpm] 0 to 400 Deviation in angular velocity [%] ± 0.1 Max. load capacity [g− kg] 9, 070 Max. payload [kg] 90.7 Max. model dimensions w × d × h [mm] 538 × 556 × 540 Table 1. Technical specifications of the centrifuge at BOKU 6 500 1000 1500 2000 200 400 600 800 1000 1200 1400 1600 y´ [px] z´ [p x] 0 5 10 15 20 25 30 35 40 45 (a) t= 10 s 500 1000 1500 2000 200 400 600 800 1000 1200 1400 1600 y´ [px] z´ [p x] 0 5 10 15 20 25 30 35 40 45 (b) t= 20 s 500 1000 1500 2000 200 400 600 800 1000 1200 1400 1600 y´ [px] z´ [p x] 0 5 10 15 20 25 30 35 40 45 (c) t= 30 s 500 1000 1500 2000 200 400 600 800 1000 1200 1400 1600 y´ [px] z´ [p x] 0 5 10 15 20 25 30 35 40 45 (d) t= 40 s Figure 6.15: Development of shear strain of test CRF-03. The solid black lines represent the shape of the slope prior rainfall initiation. drop as the wetting front descends. In addition, the cohesive and frictional components of the soil strength increase with depth. In comparison to the slope failure observed upon increasing self-weight (Chapter 5), the volume of the initial slide of the Type A failure mode is smaller due to the reduced depth of the slip surface. The evolution of the shear zone is visualised in Figure 6.15. The highest shear strain in evident at the mid-height of the slope. Similar observations are made on self-weight induced slope failures. In the figures on cumulative soil displacements (Fig. 6.14) and shear strain (Fig. 6.15) the mobilised soil mass seems to follow a straight plane inclined at about 33 degree, arising in the figures of 20 s rainfall duration. Yet, under continuous rainfall the shear zone develops under a steeper angle of about 55 degree in the upper half of the slope towards the crest region. Based on the presented pre-failure deformation with maximum values in the top region, the rain-induced slope failure - 167 - Figure 7. Snap shots of shear strain for the rainfall duration of t = 10, 20, 30, 40 s [8] Soil mechanics handles the former, fluid dynamics the latter. A model unifying both regimes remains to be developed. Acknowledgements I am grateful for funding from the European Commission within its People Programme (Marie Curie Actions) of the Seventh Framework Programme FP7/2007-2013/ under Research Executive Agency Grant Agreement No. 289911 with the title: Multiscale Modelling of Landslides and Debris Flows (MUMOLADE). References [1] H. B. Mühlhaus, I. Vardoulakis. The thickness of shear bands in granular materials. Géotechnique 37(3):271–283, 1987. doi:10.1680/geot.1987.37.3.271. [2] S. B. Savage, K. Hutter. The motion of a finite mass of granular material down a rough incline. Journal of Fluid Mechanics 199:177–215, 1989. doi:10.1017/S0022112089000340. [3] Lazari, M. (2015) Strain localization analysis of multiphase geomaterials. PhD thesis, University of Padova, Italy. [4] Switała, B. M. (2016). Analysis of slope stabilisation with soil bioengineering methods. PhD thesis, University of Natural Resources and Life Sciences, Vienna, Austria. [5] Albaba, A. (2015) Discrete element modeling of the impact of granular debris flows on rigid and flexible structures. PhD thesis, Universite Grenoble Alpes, France. [6] Leonardi, A. (2015) Debris flow simulations with Discrete Elements and Lattice-Boltzmann, PhD thesis, ETH, Switzerland. [7] C. Peng, W. Wu, H.-s. Yu, C. Wang. A sph approach for large deformation analysis with hypoplastic constitutive model. Acta Geotechnica 10, 2015. doi:10.1007/s11440-015-0399-3. [8] Idinger, G. (2016) Experimental study on failure initiation in partially saturated slopes. PhD thesis, University of Natural Resources and Life Sciences, Vienna, Austria. 34 https://doi.org/10.1680/geot.1987.37.3.271 https://doi.org/10.1017/S0022112089000340 https://doi.org/10.1007/s11440-015-0399-3 Acta Polytechnica CTU Proceedings 29:30–34, 2020 1 Introduction 2 Strain localization and triggering of slope failure by FEM 3 Stability of vegetated slope by FEM 4 Discrete Element Simulation 5 A Hybrid Method combining LBM with DEM 6 Smooth Particle Hydrodynamics 6.1 Centrifuge tests on a scaled slope model 7 Conclusion Acknowledgements References