Numerical Simulation of the Erosion Effect Caused by the Impact of High-Velocity Landslide

Due to the complex composition consisting of solid particles and fuids with diferent physical properties, geophysical fows often show complex and diverse dynamic characteristics. For landslides with high water content, there are complex interactions between the solid and fuid phases. Terefore, it is difcult to grasp the dynamic characteristics and the disaster scale of this type of landslide, especially under complex terrain and ground conditions. Te drag efect is an important aspect of the interaction between the solid and liquid phases. We optimized the enhanced drag coefcient formula to further consider the efect of high-velocity movement. By considering the volume fraction relationships between diferent phases, a mechanical erosion rate model is utilized for multiphase fows. Based on the r.avafow numerical tool and the multiphase mass fow model, considering the interphase interaction characteristics of high-velocity liquefed landslides, we analyzed the infuence of the obstruction of buildings and their entrainment into the landslide on the dynamic characteristics and hazard range of the Shenzhen 2015 landslide. Tis provides a reference for the analysis of complex geophysical disasters based on the multiphase mass fow model. Importantly, we have demonstrated the reduced mobility of the considered erosive impact event, which is in line with the physical principle.


Introduction
Geophysical mass fows such as landslides and debris fows are often characterized by high fow velocities, long-runout distances, large inundation areas, and high impact forces [1,2], they represent one of the major natural hazards threatening life and property in mountainous areas.
Under the infuence of the factors such as global climate change and engineering activities, the occurrence and frequency of the geophysical mass fow remain high in some areas. Te complex material composition and obvious difference in mechanical properties make the geophysical mass fow often show complex and diverse dynamic characteristics [3,4]. Tese factors make it more difcult to predict the potential inundation areas of these disasters. For landslide events that occur near urban areas, there are generally interactions between landslide bodies and obstacles (e.g., buildings) typically in the runout and accumulation areas [5]. Tese interactions may signifcantly afect the dynamic characteristics and the fnal disaster scale of the landslide. In this case, reasonable consideration of the obstructive efect of buildings is very important to accurately predict the scale of the disaster.
Concerning interactions between the geophysical mass fow and obstacles, a lot of attention was attracted from researchers [6][7][8]. Liu et al. [9] obtained the loads acting on the buildings in the Shenzhen 2015 landslide based on the depth-integrated continuum model, and then analyzed the infuence of the landslide body on the structure. One-way coupling (fow to buildings) was considered in their study. By using the LS-DYNA software, the efect of building blockage on the mobility of the Shenzhen 2015 landslide and the associated energy dissipation mechanisms are evaluated by Luo et al. [10]. Te same landslide has also been studied by the coupled smoothed particle hydrodynamics (SPH) method and fnite element method (FEM) [11], an element erosion algorithm is adopted to simulate the destruction process of the buildings.
Te dynamic characteristics, such as fow depth and velocity of the landslide provide the basis for disaster assessment. During the interactions between the landslide body and the buildings, the damages caused by the landslide to the buildings are controlled by the characteristics of the buildings themselves, such as structural form, material strength, and foundation depth. Te accurate analysis of these factors has gone beyond the scope of geotechnical engineering issues. Te infuence of destructible buildings on the dynamic characters and ultimate accumulation extent of the geophysical mass fow is the main focus of this study. In particular, how to carry out relevant numerical simulations and analyses in an accessible way under a unifed, advanced physics-based modeling and computational framework.
For water-laden landslides, there exist complex interactions between the solid and fuid phases [12,13], which infuence the dynamic characteristics dramatically [14,15]. Terefore, the multiphase mass fow model is a more reasonable choice for analyzing landslides with high water content. Pudasaini [13] proposed a general twophase mass fow model, which includes various interactions between the solid and the fuid constituents. Besides buoyancy, the model also includes generalized drags between phases, enhanced non-Newtonian viscous stress, and virtual mass. Tese properties make the two-phase mass fow model of Pudasaini [13] and its extensions able to model a wide range of geophysical mass fows [4,[16][17][18][19].
Numerical modeling is a suitable method for complex geotechnical mass fow and is an important method for back-analyzing disaster events and exploring the mechanisms of disasters [2,14]. Te applicability and efectiveness of the numerical models should be verifed by experiments or site events. Te tool r.avafow [20] that will be utilized in this study is designed as a raster module of the GRASS GIS software. It employs the C programming language in the core computation module and the Python and R statistical languages for auxiliary work such as data transfer, analysis, and graphing. A multiphase mass fow model of Pudasaini and Mergili [4] with an enhanced generalized drag model [21] and a virtual mass force model [22] are utilized in the most recent version (V2.4 prerelease) of r.avafow. Tese models are suitable for the problems covered in this study. Furthermore, the mechanical erosion rate model of Pudasaini and Fischer [23] was implemented in the r.avafow tool and was utilized for modeling the 2015 Shenzhen landslide (China). Te destruction of the buildings caused by landslides is approximated by the erosion process. Trough this method, the infuences of destructible buildings on the landslide's dynamic characteristics and the fnal accumulation state are analyzed.

Mass and Momentum Conservation Equations.
Tere are some important aspects that should be considered in the analysis of the geophysical mass fow. Tese include the diverse composition of the material, the complex interactions between the phases, the possible liquefaction resulting from the high-velocity movement, and fnally, the variations of the internal components caused by erosion and entrainment along the path and the runout. Tese aspects are important characteristics and also the source of difculty in predicting the hazard scale of geophysical mass fows. Pudasaini and Mergili [4] have proposed a fexible multiphase mass fow model that includes enhanced drag force [21] and virtual mass force models [22]. By considering the diferences in mechanics and rheological characteristics between the coarse and fne solid particle and viscous fuid, this multiphase mass fow model can better represent the real geophysical mass fows, and it has wider and diverse application prospects as proven by their successful applications in accurately simulating complex cascading catastrophic landslide events [18,24].
Adopting the general application of the Pudasaini and Mergili [4] three-phase model, in this study, the buildings on the sliding path were regarded as the third phase, in addition to the solid and fuid phases of the landslide body. Te third phase represents the obstacles that can be eroded and entrained into the landslide body. Te third phase follows a Coulomb-viscoplastic [25] material behavior, which may have signifcant viscosity and yield strength. Te yield criterion helps distinguish between the fow and nonfow regions of the mass [4]. Terefore, it helps to obtain diverse depositional states of the geophysical mass fow. When the buildings are destabilized and erosion occurs, the obstacle phase (buildings) is entrained into the landslide body and the higher viscosity makes the obstacle phase hardly expand laterally [4]. Te damaged buildings can hardly slide with the landslide over a long distance. Terefore, a solid phase with Coulomb-viscoplastic rheology is suitable for approximating the efects of entrainable obstacles in the landslide simulation. During the erosion process, there are jumps in physical quantities such as velocity and density on the erosion interface. Te contribution of obstacle erosion to the mobility of the landslide is afected by the competition of physical quantities on both sides of the erosion interface. With the help of the erosive landslide mobility equation of Pudasaini and Krautblatter [26], the contribution of erosion to landslide mobility can be analyzed based on mechanical principles.
Based on the three-phase mass fow model of Pudasaini and Mergili [4] and including erosion and entrainment [23,26], the following model equations are adopted for the conservation of mass, for the solid (s) phase of a landslide, fuid (f ) phase of a landslide, and the solid phase of the obstacle (o): where following mechanical principles of Pudasaini and Fischer [23] and Pudasaini and Krautblatter [26], the velocities of the solid, fuid, and obstacle mass that have just been eroded from the bed (in the x and y directions) are denoted as , respectively, and are called the erosion velocities of the respective materials from the bed. Generally, the eroded mass does not directly gain the same velocity as the landslide body. It is also mechanically incorrect to simply set the velocity of the eroded mass to zero or ignore it [23]. Terefore, it is necessary to distinguish the velocity of the entrained mass from that of the landslide mass [26]. Pudasaini and Krautblatter [26] frst found and emphasized that the mobility of erosive landslides depends on the ratio of erosion velocity to landslide velocity. u vm s , u vm f , and u vm o are the virtual mass induced mass and momentum enhancements for the solid, fuid, and obstacle-solid phases, respectively.
are hydraulic pressure coefcients for solid, fuid, and obstaclesolid phase respectively, and K is the Earth pressure coeffcient, g z is the component of gravitational acceleration in the z-direction (similarly, g x and g y are the gravitational acceleration components in the x and y directions, respectively), and c f s � ρ f /ρ s and c f o � ρ f /ρ o are the density ratios of the corresponding phases. Due to the symmetry of model equations, the y-directional momentum equations have a similar form. Te defnitions of the source terms and associated parameters on the right-hand side of these momentum equations are presented below.

Te Source Terms and Relevant Parameters.
Te source terms in x-directional momentum equations can be expressed as [4,23] Shock and Vibration In these equations, C i,j DG is the drag coefcient between the two phases that are indicated in the superscript with i and j (where i, j � s, f, o, and i ≠ j), and τ nN is the enhanced non-Newtonian viscous stress [4]. Te value of the parameter J should be determined based on whether the fow has a laminar (J � 1) or turbulent (J � 2) character [13]. It can be seen from the above expressions of the source terms that the original model treats both fne particles and fuids as non-Newtonian fuids with a yield. High viscosity and non-Newtonian fuids can have lower mobility than solid particles. Buildings have a specifc strength, so they fail and collapse only when the impact stress exceeds their shear strength. In addition, buildings generally have higher density due to reinforcements, and the broken parts have lower mobility, which is similar to the low mobility of high-viscosity non-Newtonian fuids. Terefore, the high-viscosity, non-Newtonian fuids were utilized to represent buildings impacted by landslides in this study.
Te generalized and smooth enhanced drag coefcients that appear on the right-hand sides of (8)-(10) are analytically derived by Pudasaini [21]: In equation (11a), U s,f T represents the terminal velocity of a solid particle falling in the fuid, are parameters representing the proportion of the fuid-like and solid-like drag characteristics in the total drag force, respectively. P s,f ∈ (0, 1) combines the fuid-like and solidlike drag contributions between the solid and fuid phases. Te last item S s,f P � (P s,f /α s | + |(1| − |P s,f )/α f )K s,f in the denominator is the smoothing or damping function that removes the singularity when the drag coefcient evolves with the proportion of the phase. Tis new function [21] removes the singularity of the previous version [13] and ofers a smooth variation of the drag coefcient as the solid volume fraction evolves. Te variable K s,f depends on the total mass fux, which is inversely proportional to the drag coefcient. A higher value of K s,f means higher mass fux, which in turn leads to lower drag. Although the parameter P s,f can be regarded as a numerical parameter, it can also be related to physical quantities such as solid volume fraction [19,21]. Te other two drag coefcients (in equations (11b) and 11c) have a similar form to equations (11a) and refect the physics-based interaction between the two phases. We refer to Pudasaini and Mergili [4] for the expressions of the other drag coefcients.
Te parameter K refects the quantity of mass fux and may vary with the dynamic state of the mass fow [21]. Terefore, dynamic parameters such as fow rate, volume fraction, and shear rate have a potential infuence on the value of the parameter K. For example, in the accumulation zone, the velocity of the material fow is slowed down, and the material fux should decrease accordingly, which corresponds to an increase in drag. A feasible way to determine the value of K is to correlate it with the fow velocity, as suggested by Pudasaini [21]. An expression based on the hyperbolic tangent function is proposed in this study to construct a mass fux that depends on the dynamic state of the landslide. For simplicity, the value of mass fux will be related to the bulk velocity of solid-fuid mixture fows: where the variable v represents the bulk velocity of the landslide mixture, K max and K min are the upper and lower bounds of the mass fux, respectively. And tanh (v) is the hyperbolic tangent function which returns values close to -1 for v ⟶ − ∞ and values close to 1 for v ⟶ + ∞. Te parameter v o shifts the velocity in which the K v value changes. For example, when K min � 1.0 and K max � 2.0, for specifed values of v o and m, the curves for the variations of mass fux K are demonstrated in Figure 1. It can be found that the parameter m afects the width and smoothness of the smooth transition interval. As the value of m increases, a larger transition interval and a smoother transition will be obtained. As can be seen from Figure 1, velocity-based mass fux functions refect that the changes in K only occur within a certain velocity interval. When the velocity is below the lower bound of the interval, the mass fux remains unchanged, indicating that the increase in velocity has not yet caused further changes in the interaction state between the two phases. When the velocity exceeds the upper bound of the interval, the mass fux K no longer increases, which may correspond to a fully liquefed state. For the models of the virtual mass forces, efective viscosities, pressure-and rate-dependent Coulomb-viscoplastic rheological laws, and non-Newtonian fuid stresses, please see Pudasaini and Mergili [4].

Erosion-Entrainment Model. Experiments and feld observations show that erosion and subsequent entrainment
can signifcantly change the dynamics of geophysical mass fow [27][28][29]. Tis has been proven by Mergili et al. [18] and Shugar et al. [24] who applied the Pudasaini and Mergili [4] multiphase mass fow model to simulate those catastrophic natural landslide events. Generally, during the erosion process, the eroded material along the path is accelerated to a certain velocity, and part of the material is entrained into the landslide body. Tis fact is revealed for the frst time with the mechanical erosion and entrainment model by Pudasaini and Fischer [23] and Pudasaini and Krautblatter [26]. In this process, momentum exchange occurs between the eroded material and the landslide body. Similarly, when the landslide body encounters obstacles such as buildings, the obstacles are severely impacted and destroyed, and parts of the damaged objects are mobilized and entrained into the landslide body. Buildings fail only when the shear stress generated by the impact of the landslide exceeds their shear strength. Te damaged part is detached from the main part of the building. Te detached parts will move with the landslide only after satisfying certain mechanical conditions. Tis process is similar to the normal erosion process. Tis has been proven with mechanical erosion and entrainment model by Pudasaini and Krautblatter [26].
Depending on factors such as the water content, porosity, density, and friction of the eroded material, erosion and entrainment may lead to the acceleration or deceleration of landslides in the natural environment. Tis explicitly depends on how the erosion velocity is related to the fow velocity, i.e., u b � λ b u, as mechanically proven by Pudasaini and Krautblatter [26]. Treating collapsible buildings as a kind of erodible terrain provides a convenient method for analyzing the infuence of buildings on the landslide accumulation area.
Te erosion process of the geophysical mass fow may be driven by a variety of potential mechanisms [23,26]. One of the most common types of erosion is basal erosion, in which, the development of erosion is controlled mainly by the shear stress on the erosion interface. Te focus of this study is to consider the landslide impact on buildings, their failures, and on the characteristics of landslide motion in an appropriate way. Te buildings on the path of the landslide directly hinder the movement of the landslide, and the kinetic energy is consumed when the landslide body hits the building. If the buildings do not crumble and collapse throughout the whole process, they can be approximated as a part of the terrain in the simulation, e.g., by increasing  friction. However, the huge impact of the landslide often causes a part of the building to be broken and collapse, and even parts of the building will be detached from the foundation and move along with the landslide. After collapse or displacement, the height and shape of buildings change, and their impact on the subsequent landslide mass also changes. Terefore, to simulate the infuence of buildings on the landslide dynamics, it is necessary to consider the variation in the mobility of landslides caused by the impedance from the buildings. Te variations in the height, shape, and position of the collapsed buildings should also be considered.
Tis study utilizes the erosion rate formula proposed by Pudasaini and Krautblatter [26] to multiphase mass fow, that is, we treat buildings as erodible solids with Coulombviscoplastic characteristics in addition to the coarse-solid phase and fuid phase that make up the landslide mass. Te coarse-solid phase was named by Pudasaini and Mergili [4] within their multiphase mass fow model, which represents a class of solid matter that obeys the Coulomb friction law. In this study, this phase represents all solid particles in the landslide body, regardless of their particle sizes.
When landslide encounters obstacles (e.g., buildings), the landslide mass impacts the surface of the obstacles and fows over the surface with a high velocity ( Figure 2). Te damage and even collapse of obstacles caused by shearing can be approximated as a kind of erosion. Te momentum fux at the eroded interface is generated by the applied shear stress. During the erosion process, the presence of frictional net shear stress must be balanced by the net momentum fux [23]: where the subscript so represents the total solid phase mixture, which consists of the solid phase of a landslide and the solid phase from obstacles. Terefore, E so represents the erosion rate of the total solid phase mixture. All the solid phases at the shear interface obey Coulomb's friction law, regardless of whether the phase belongs to landslides or obstacles. Terefore, the net shear stress can be expressed as [13,23] where ζ is the slope angle ( Figure 2) and α m so and α b so are the solid mixture volume fractions on either side of the erosion interface, that is, so are the fuid-solid density ratios of the landslide body and bed, respectively, ρ m so and ρ b so are the solid mixture densities of the landslide body and the obstacle, respectively. ρ m so and ρ b so can be determined based on the volume fraction relationship between the two solid phases within the landslide and obstacle, respectively, that is, . μ m so and μ b so are the Coulomb friction coefcients of solid phases within the landslide and the obstacle, respectively, which can also be determined based on solid volume fractions, that is, . Te shear velocity can be expressed as the square root of the ratio between the net shear stress and the relative net density at the eroded interface, which is proportional to the mean fow velocity of the solid phase u s [23]. With this relationship, the erosion rate can fnally be expressed as where E 0 � E so has been realized, and ] is the shear velocity factor that represents the transformation factor between the efective erosional shear stress and the efective velocity jump at the erosion interface. For the fuid erosion rate, see Pudasaini and Fischer [23]. Tis erosion rate formula is suitable for multiphase mass fow and is mainly based on the mechanically correct and mathematically consistent erosion rate model of the twophase mass fow [23,26].  [23] and Pudasaini and Krautblatter [26]. Te fow depth of the landslide is h, the vertical coordinate of the top of the bed is z b , the average velocity of the landslide mass is u, while λ m so and λ b so are the erosion drift coefcients that relates the velocity at the bottom of the landslide and the velocity at the top of the bed to the average fow velocity. Shear stress generated by the mixture of the solid and obstacle phase on the bed is denoted as τ m so , and conversely, the shear stress produced by the bed on these two phases is τ b so .

Shock and Vibration
Te most outstanding contributions of Pudasaini and Mergili [4] are to extend the two-phase mass fow model to three phases and to take into account the complex interactions between the phases with distinct properties. In particular, the introduction of a fne-solid phase with Coulomb-viscoplastic mechanical behavior enables the model to refect the more complex and diverse dynamics of the real geophysical mass fow. Tis study directly and elegantly utilized the mechanically based, mathematically consistent erosion rate models developed by Pudasaini and Fischer [23] and Pudasaini and Krautblatter [26] to the multiphase dynamic model.
In the following study, the impact area of a landslide is mainly located in an industrial park, and most of the ground surface within the landslide movement range is hardened. Terefore, ground surface erosion caused by landslides is negligible, and only the erosion of the buildings on the runout path will be considered. Tis is modelled by considering E s = 0, E f = 0, and E o is given in equation (15). Tese expressions are then utilized in equations (1)-(7), modeling the erosion induced mass and momentum productions.

The 2015 Shenzhen Landslide and Its Modeling
A landslide occurred in the Guangming New District of Shenzhen (China) on December 20, 2015, which was caused by the instability of a large artifcial landfll [30]. Approximately 2.73 × 10 6 m 3 of construction spoil was mobilized with an impact area of 0.38 km 2 in this landslide (Figure 3(b)), and the maximum runout distance is about 1100 m [31]. Tis landslide destroyed 33 buildings, caused 73 fatalities, and four people missing. Te natural slope of the mountain in the inundation area is in the range of 25-35°. How the sliding body avalanches such a long distance on such relatively fat terrain and obtains such a high velocity is the main focus of the researchers [15,19,[30][31][32]. Te source area of the landslide was originally a quarry, within which a deep pit was formed due to long-term mining. Te bedrock of the pit is mainly Cretaceous granite rock with lower permeability [32]. During the abandonment period, a large amount of water accumulated in the bottom of the quarry pits, which was mainly rainfall and surface water. Years later, the quarry was used as a landfll for municipal construction waste. Te stagnant water at the bottom of the quarry pit was not pumped out before flling with construction waste. Tis results in higher water content in the landfll soil, especially near the bottom. Te flled soil is mainly silty soil and clay with a small amount of gravel grain [31].
Studies have shown that the high-velocity and longrunout distance characteristics of the Shenzhen 2015 landslide are mainly due to the high-water content and high initial pore pressure in the source area of the landslide [15,30,31,33]. Te low permeability of soil results in a signifcant increase in pore pressure during the sliding [15]. Terefore, when simulating this landslide, the single-phase fow model needs to use a very low basal friction coefcient to obtain results that are more consistent with reality [32,34]. Such an unrealistically low friction was needed in those simulations because erosion-induced net momentum production, which produced excess kinetic energy and resulted in higher mass fow mobility, could not be considered in those simulations. Tis has been proven in Pudasaini and Krautblatter [26] by presenting a novel mechanical model for landslide mobility with erosion. For the two-phase fow model, it is necessary to consider the high fuid phase ratio and the fuidization characteristics in terms of drag force [19].
At least 33 buildings have been damaged in this landslide since the impact area of the landslide is located in the industrial park. Terefore, the infuence of buildings is an important factor afecting the accumulation range of the landslide and should be considered reasonable. Due to the complexity of the problem, the infuences of buildings were not considered in most of the previous analyses [15,19,32,34]. However, if those buildings had not been mobilized, energy loss due to the impact would have been higher, resulting in reduced mobility. Tis important fact has been proven by Pudasaini and Krautblatter [26] with their erosion-induced landslide mobility model. In this study, the collapsible buildings are regarded as a low-mobility phase and the infuences of obstacles on landslides are approximated through the erosion-entrainment efect.

Numerical Simulation Tool.
In this study, the numerical tool named r.avafow (V2.4 Pre-release) was utilized. r.avafow [20] is a fexible and GIS-based open-source computational framework for the simulation of geomorphic mass fows and process chains [17,18]. Tis simulation tool is based on the Pudasaini and Mergili [4] multiphase mass fow model and is wrapped in Python and R language for data management, pre-and postprocessing tasks, while the core of the computational algorithm is implemented in the C programming language. In r.avafow, the custom erosion and entrainment model can be defned as a complementary function for the source items. Te maximum depth of entrainment can be defned by the user as a raster map. At the end of each time step, the erosion and entrained depths are updated to the basal and fow topography simultaneously. At the same time, fow momentum is updated accordingly.

Model Parameters of 2015 Shenzhen Landslide.
Te hydraulic pressure coefcient β x s (� |K x s |g z |(1| − |c f s )) in equation (5) shows that buoyancy reduces the solid normal and lateral stresses by a factor of (1| − |c f s ). Terefore, a greater buoyancy (pore fuid pressure) efect can be refected by a lower solid phase density value [13,35]. To account for the high-water content of the soil in the source area of the Shenzhen 2015 landslide, a relatively smaller solid phase density (ρ s ) should be adopted. Following the study of Qiao et al. [19], the densities ρ s � 1800 kg/m 3 and ρ f � 1100 kg/m 3 are used for the solid and fuid phases of the landslide body. Obstacles have a higher density. Te equivalent bulk density of the buildings cannot be equivalent to the density of concrete because there are large empty spaces inside the buildings. Terefore, a density value higher than the landslide density and lower than the concrete density is assigned to the obstacle, that is, ρ o � 2000 kg/m 3 . Based on previous feld surveys and analyses [30,31], the solid volume fraction of the soil in the landslide source area is taken as 0.58 in this study.
Drag force is an important factor afecting the interaction between multiple phases. Parameter P has a negative efect on the advection and dispersion of the fow [13]. It can be associated with specifc physical quantities, such as the solid-phase volume fraction. Following Pudasaini and Mergili [4] and Qiao et al. [19], P � α s is adopted in this study. By assuming a velocity drift parameter λ, Pudasaini [21] evaluates the upper bound of the mass fux contribution as κ⩽(α s | + |λ|α f )|u s |. For natural debris fows or fuidized landslides, the velocity of the solid phase is usually at the level of 10 ms − 1 [21], and the velocity drift parameter is a value close to unity, so the upper bound of K can be estimated as 10 ms − 1 . Te default value of the fux parameter K is 1.0 in r.avafow. For the landslides with high-water content in this study, to refect the infuence of the landslide motion state on the drag force, the fux parameter K in the form of equation (12) was adopted, and K min � 1.0 and K max � 2.0 were used.
In Figure 4, three drag coefcient curves with constant mass fux K are compared following the analytical, smooth, and enhanced drag model developed by Pudasaini [21]. Te diference between the drag coefcient curves of K � 1 and K � 2 gradually increases as the solid volume fraction falls within the usual range for the natural geophysical fow, i.e., 0.4 < α s < 0.6. Te drag coefcient curve with K � 0 corresponds to the old drag function of Pudasaini [13], which rapidly tends to an infnite value as the solid phase volume fraction gradually increases. Te new drag function [21]   Shock and Vibration provides a full-domain representation of the drag that conforms to the physics. When the bed with a pure solid phase is eroded, the entrained solid phase causes the solid phase volume fraction of the landslide body to increase signifcantly. Te drag feature of this dense regime can be described by the new drag function [21]. In this study, the drag curves with K � 1 and K � 2 can be regarded as the upper and lower bounds of the drag functions with velocitydependent mass fuxes. When the solid phase volume fraction is in the normal interval (0.4 < α s < 0.6), the variation in drag coefcient caused by the velocity factor is relatively small. Equation (12) can be regarded as a conservative approximation of the actual drag coefcient. By evaluating the existing simulation results of the Shenzhen 2015 landslide [15,19], the pore water pressure in the landslide body has a further increase in the velocity interval (10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20). Te relationship between the pore water pressure gradient and mass fux can be approximately regarded as a linear relationship [1]. Terefore, it can be approximately considered that in the interval (10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20), with the increase in velocity, the mass fux will increase to a certain extent. So, the middle value of the velocity interval v 0 is estimated to be 15 m/s, and m � 3 is adopted, which means the velocity transition curve is wider and fatter (Figure 1). To replace the original constant mass fux K with equation (12), the source code of r.avafow should be modifed accordingly.
Te soil in the source area has a high content of fne particles. Te sieving experimental data shows that the volume of solid particles with a particle size below 0.1 mm accounts for about 40% of the total solid particle volume. Tis particle composition leads to a relatively small internal friction angle and basal friction angle. Yin et al. [31] measured the internal friction angle of some soil samples in the source area to be 26°, and the measured internal friction angle used by Ouyang et al. [32] was 31.9°. If the bottom friction angles of coarse and fne particles are estimated to be 25°and 7.5°, respectively, the equivalent bottom friction angle of the total solid phase is about 18°according to the volume fraction of fne particles. Based on the above analysis, in this study, the internal friction angle and basal friction angle of the solid phase, which include both coarse and fne solid particles in the landslide body, are set as 31°a nd 18°, respectively. For fuid with a high content of fne sediment, the viscosity η f = 2 Pa.s was adopted in the simulation. Te terminal velocity UT is 0.1, and the particle Reynolds number Rep is 1. Te exponent J for drag in equation (11) takes the value 1 for laminar-type fow. Te viscosity η o = 10 2 Pa.s and yield strength τ yo = 40 Pa are adopted for the obstacle. Most of the other parameter values (e.g., those in the "special" parameter list) have adopted default values that have been verifed [8,13,17,18,26].
Based on the feld data and images from the 2015 Shenzhen landslide, the main damaged buildings on the digital topographic map were restored (Figure 5(a)). By comparing the digital elevation map of the source area before and after the landslide, the soil height of the source area can be determined, and the maximum height value is about 44.6 m. Te maximum height of the destroyed building was about 16.7 m. Buildings that were damaged or collapsed in the landslide are marked by solid blue lines in Figure 5(a), and these buildings will be treated as erodible bed in the simulation. Uncollapsed buildings are shown in Figure 5(b), and they will be treated as nonerodible bed in the simulation. For the erosion rate model, most of the parameters are the physical parameters of the landslide mass or basal bed, for example, density and volume fraction. Te fnal average accumulation thickness of the landslide is about 8 m, and the thickness of the landslide in fow is smaller than this value. Terefore, it can be considered that the distribution of landslide velocity in height conforms to the hypothesis of plug fow, that is, λ m s ≈ 1. Te erosion drift coefcient in the basal bed can be determined based on the physical quantities on both sides of the erosion interface with the mechanical erosion drift function [23], that is, Te shear velocity factor ] satisfes 1/ � ] √ � 0.05, it describes the rate at which the eroded mass is entrained into the fowing mass, or depositional mass is removed from the fowing mass. Te value of this parameter follows the recommendation of the software manual. Te basal friction angle of the bed is 8°, which is 10°s maller than that of the solid phase in the landslide. Due to the complexity of the buildings themselves, the value of this angle needs to be determined based on trial calculations.

Numerical Results and Analyses.
Due to special geological and topographic factors, the source mass of the Shenzhen 2015 landslide has a high-water content [31]. Terefore, once the landslide is triggered, it will attain a higher velocity soon. Due to the high-water content of the actual soil, a higher fuid volume fraction is given to the soil of the source area in the simulations. For a multiphase fow model, due to the high initial volume fraction of the fuid phase, a signifcant increase in buoyancy and a reduction in friction will be obtained. Terefore, once the simulation starts, the mass in the landslide source area will displace immediately. To compare the infuence of buildings on the dynamic characteristics of the landslide, the buildings within the infuence range of the landslide are regarded as erodible masses with Coulomb-viscoplastic characteristics in simulation (I). In simulation (II), the infuence of destroyed buildings within the impact range of a landslide will not be considered, that is, the destroyed buildings will be replaced by the ground surface. Buildings outside the impact range of a landslide will be regarded as nonerodible bed. In r.avafow, both erodible and nonerodible buildings are a part of the initial terrain. Trough the numerical simulation of r.avafow, the impact range and dynamic characteristics of the landslide at diferent times are obtained.

Accumulation Range of the Landslide.
In simulation (I), after the triggering, the landslide mass in the source area fowed out from the low-lying region on the north side of the landfll. Te front of the landslide was pushed forward quickly as it approached the buildings (t 25 s). At this time, only the left side of the leading front of the landslide impacted the erodible buildings ( Figure 6). After this, the Shock and Vibration sliding mass was further hindered by the buildings on the sliding front. Te southern edge of the buildings on the southern side of the industrial park is approximately on a straight line (green dashed line B-B′ in Figure 6), so the frst impedance from the buildings occurred on the west side of this dashed line, and the subsequent impedance occurred eastward along this dashed line. Te fow on the east side of the landslide was hindered last, so this part of the fow had the longest runout distance. Due to the obstruction of the buildings, the fow direction has been altered. At the same time, the impact of the landslide caused the mass of the buildings to be detached (the eroded mass is marked in green in Figure 6 at t � 60.0 s). With the subsequent accumulation of the landslide material at the sliding front, the accumulation height increases at the sliding front and the accumulation range continues to expand. Te resulting erosion extent is further developed. On the most western side of the industrial park, low-rise buildings were eroded and carried into the sliding mass. Te entrained obstacle phase has a smaller volume fraction in the local total fow height, which is similar to collapsed low-rise buildings being buried by a large amount of subsequent landslide mass ( Figure 6 at t � 100 s). At about t � 200 s, the morphology of the accumulation area was nearly formed. Te amount and velocity of the landslide mass that subsequently fowed into the accumulation area decreased signifcantly ( Figure 6). Due to the small bed slope in the source area (only a few degrees), a portion of the material remained in the source area at the end of the simulation. Te investigation showed that the average thickness of the remaining material in the source area was about 20.5 m [30].
In simulation (I), after the destruction, the accumulation is close to the original position of the destroyed buildings.
Te eroded obstacle phase slides away with the landslide only a little. Tis displacement feature is close to the actual situation.
For comparison, simulation (II) was carried out without considering the buildings (inside the region with a blue border in Figure 5(a)) that were damaged in the actual landslide, that is, simulation (II) used the digital elevation map shown in Figure 5(b). In addition, the obstacle phase need not be considered in simulation (II), only the solid and liquid phases in the landslide mass will be considered. Simulation (II) also does not need to account for erosion, and the erosion term in the model equations will be zeroed out. Except for these three aspects, the rest of the parameters are the same as in simulation (I).
In simulation (II) (Figure 7), the landslide slid about 222 m longer than in simulation (I) before encountering the erodible buildings. Te change in the fow direction of the landslide occurred later than in simulation (I). So, the landslide achieved signifcantly larger runout distances than in simulation (I) at the same time. Tis phenomenon is especially obvious on the right side of the landslide front (see Figures 6 and 7). Te obstacles in the red dashed circle in Figure 6 were broken by the impact of the landslide (the eroded obstacles are in green), and the landslide can only bypass the obstacles due to the obstruction. In Figure 7, the fow characteristics of the landslide difer from reality due to the absence of this obstruction. By comparing the actual fow coverage and accumulation ranges in the two simulations (represented by the red solid line in Figures 6  and 7), it is observed that simulation (I) obtains an accumulation range that is more consistent with the actual situation after considering the entrainment of obstructing buildings. Since the maximum velocity is gained at the foot of the mountain, and at this location in simulation (I), the landslide has not yet encountered the buildings. Terefore, simulation (II) and simulation (I) have the same maximum sliding velocity, which is 30.1 m/s. Tis value is very close to the maximum velocity value estimated by Yin et al. [31], which is 29.8 m/s.

Evolution of Erosion Characteristics.
For simulation (I), the evolution characteristics of the erosion extent during the erosion development stage are shown in Figure 8. Te buildings within the impact range can be divided into three groups by the yellow dashed lines (Figure 8). Te buildings in Group A are located in front of the landslide movement, and the erosion development direction (cyan arrow in Figure 8) is nearly the same as the movement direction. Buildings in Group C are located on the east side of the landslide movement, and the erosion development direction is almost perpendicular to the direction of the landslide movement. Buildings in Group B are located sporadically and both frontal and lateral erosion have occurred. From t � 35 s to t � 40 s, the erosion extent continued to expand gradually. At t � 50 s, most of the buildings in Group A were eroded. In Figure 8, near the south of the erodible buildings (locally enlarged area), due to impedance, the velocity vector arrows of the fowing mass show a signifcant change in fow directions. Tis change is consistent with the shape of the fnal accumulation (represented by the solid red line in Figure 8). Te evolution characteristics of erosion can be further compared in the profle section (yellow line in Figure 6 or Figure 3(b)). In Figure 9, the fow height is denoted by h o for the mobile obstacle phase eroded from the buildings. Once they have lost mobility and accumulated, they will become a part of the bed again. Obstacles that have not eroded and those that have lost mobility are considered parts of the bed, whose top elevation is denoted as z b . In the early stage (Figure 9), erosion develops progressively along the direction of the landslide movement. Te climb of the landslide along the building surface results in a greater amount of erosion at the top of the buildings. In the later stage (Figure 9), at t � 50 s, most of the front buildings have eroded. At 100 s, the building has been completely eroded. Since the eroded obstacle phase follows a viscoplastic material criterion with viscosity and yield strength, it is able to maintain a steeper slope and less deformation after erosion. By comparing the h o at t � 100 s and t � 150 s, it can be found that, eroded obstacles displace very slowly in the later stage.
3.6. Characteristics of the Velocity. Te velocity characteristics of the solid-phase along the profle line (yellow line in Figure 6 or Figure 3(b)), in the two simulations are  compared in Figure 10. It can be found that, at the same moment, the solid mass of the landslide obtains the maximum velocity at the foot of the slope (the red point in Figure 10). Before the landslide encounters the erodible obstacle (t � 25 s), the velocity distribution in both the simulations is the same. After the landslide encounters an erodible obstacle, two curves (corresponding to t � 30 s and t � 35 s, respectively) have a sharp drop in the velocity at the leading front of the landslide. At t � 25 s, the impedance of the buildings resulted in a slightly smaller maximum runout distance in simulation (I). At t � 35 s, as the erosion continued, the decline in the velocity of the leading front of the landslide increased further. Te gap in the landslide leading front between simulation (II) and simulation (I) increased further.
In simulation (I), the sliding of the landslide was hindered when the landslide impacted the buildings. Te erosion process of the obstacle is determined by the relationship between the net stress and the net momentum fux on both sides of the erosion interface. Te evolution curves of the total kinetic energy of the landslide are shown in Figure 11 (right vertical axis). In simulation (I), the landslide was hindered by the building since t � 30 s, at which time the landslide gained the maximum kinetic energy. In simulation (II), the total kinetic energy of the landslide maintained a relatively rapid increase until t � 35 s. After that, due to the infuence of terrain and frictional resistance, the growth rate of total kinetic energy decreased slightly. At t � 40 s, the landslide hit the buildings and the total kinetic energy began to drop. Te comparison found that the total kinetic energy in simulation (I) ended the upward trend early due to the impedance of the buildings, and the total kinetic energy in simulation (I) was always smaller than that in simulation (II) due to the subsequent continuous impedance.
Te red curve in Figure 11 represents the evolution of the erosion volume (left vertical axis). Te erosion volume (representing the amount of destroyed buildings) started to increase since t � 25 s, and the landslide eroded the part of the buildings on the left side of the sliding direction  Figure 9: Te evolution characteristics of the erosion height and bed elevation on the profle section (along the yellow dashed line in Figure 6). ( Figure 6). As the landslide progresses, the landslide mass eroded a large number of buildings in front of it, and the eroded volume increases signifcantly. Te erosion volume reaches its maximum value at t � 90 s and then decreases slowly. Te decrease in eroded volume is caused by the accumulation of the obstacle phase. Te maximum total erosion volume is about 50 × 10 4 m 3 , which is about 18.3% of the source material volume, which is substantial for the situation considered here with mobilized buildings. Terefore, the accumulation range and morphology in simulation (II), which does not consider this volume, are quite diferent from the actual ones.

Discussion
Geophysical mass fows in nature often encounter various natural or artifcial obstacles on their movement path, such as houses, retaining dams, and diversion dikes. When analyzing the dynamic characteristics of the geophysical mass fow, artifcial obstacles that have not been destroyed can be regarded as rigid terrain. But when a part of the obstacle breaks down and moves, the interaction between the obstacle and the geophysical mass fow becomes complicated. In this case, a numerical tool based on a unifed computing architecture becomes necessary. Tis complex analysis is made possible by the multiphase geophysical mass fow model [4], the mechanical erosion rate model [23] and the erosive landslide mobility theory [26], and the associated numerical implementation in r.avafow.
In the present study, the solid phase representing the obstacle follows Coulomb-viscoplastic mechanical behavior, and the development of its erosion is governed by the competition between the net stress and the net momentum fux on both sides of the erosion interface. Te erosion drift coefcient λ b so is determined based on the physical quantities on both sides of the erosion interface, i.e., λ b so � λ m so /[1| + |ρ b so α b so /(ρ m so |α m so )]. Since ρ b so > ρ m so , α b so > α m so , and λ m so ≈ 1, then λ b so < 1/2. As has been demonstrated by Pudasaini and Krautblatter [26], when λ b so < 1/2, the landslide loses kinetic energy during the erosion process, that is, the obstacle hinders the sliding of the landslide. Tis demonstrates the functionality of the state of the erosioninduced mobility presented by Pudasaini and Krautblatter [26]. Tis clearly manifests the practical and engineering importance of this research work that practitioners can take direct beneft from it.
Most of the parameter values in the erosion model are determined by the physical properties of the material. Due to the complex characteristics of the eroded object, the bottom friction angle of the obstacle is determined based on the physical characteristics and trial simulations. Figure 12 shows the depositional characteristics of obstacles on the profle section with diferent yield strengths. It is observed that the depositional characteristics of eroded obstacles are afected by their yield strength. Te mass with lower yield strength will be pushed further forward in the sliding direction of the landslide. When the yield strength exceeds a certain value, the depositional position of the obstacles is almost unchanged.

Conclusions
Te drag efect between the solid and liquid phases is an important aspect that afects the dynamic characteristics of the multiphase mass fow. To consider the evolution of the drag coefcient [21] within a landslide, a new formula for the mass fux parameter for Pudasaini's drag coefcient was proposed. Tis new formula refects the infuence of the landslide velocity on mass fux parameters. Terefore, it is suitable for the simulation of high-velocity landslides with higher water content.
In this study, the mechanical erosion rate model of Pudasaini and Fischer [23] was utilized for the multiphase mass fow. Te parameters in the erosion rate model are easily determined as general physical properties. For the special problem of eroding buildings considered here, buildings are described with Coulomb-viscoplastic material. Te goal is to approximate the movable impedance from collapsible buildings with an erodible basal bed. To obtain an accumulation range that is consistent with the actual landslide, the basal friction coefcient of the entrained buildings should be smaller than that of the landslide body. Te entrainment phenomenon has been explained by the theory of erosive landslide mobility of Pudasaini and Krautblatter [26].
Te multiphase mass fow model [4] which is applied to a GIS-based computational tool r.avafow has been extended to simulate the complex dynamics of the 2015 Shenzhen landslide event with the process-based computing framework. For disaster prevention and mitigation, it is very important to obtain reliable fow mechanical characteristics and an accumulation range of geophysical mass fow through efective numerical tools. Tis study provides a method to approximately consider the infuence of erodible obstacles such as buildings and constitutes a useful reference for the prediction and assessment of geophysical disasters associated with landslides and debris fows. 14 Shock and Vibration

Data Availability
Te data supporting this study are available upon request from the corresponding author.

Conflicts of Interest
Te authors declare that there are no conficts of interest regarding the publication of this paper.