Due to the complex composition consisting of solid particles and fluids with different physical properties, geophysical flows often show complex and diverse dynamic characteristics. For landslides with high water content, there are complex interactions between the solid and fluid phases. Therefore, it is difficult to grasp the dynamic characteristics and the disaster scale of this type of landslide, especially under complex terrain and ground conditions. The drag effect is an important aspect of the interaction between the solid and liquid phases. We optimized the enhanced drag coefficient formula to further consider the effect of high-velocity movement. By considering the volume fraction relationships between different phases, a mechanical erosion rate model is utilized for multiphase flows. Based on the r.avaflow numerical tool and the multiphase mass flow model, considering the interphase interaction characteristics of high-velocity liquefied landslides, we analyzed the influence of the obstruction of buildings and their entrainment into the landslide on the dynamic characteristics and hazard range of the Shenzhen 2015 landslide. This provides a reference for the analysis of complex geophysical disasters based on the multiphase mass flow model. Importantly, we have demonstrated the reduced mobility of the considered erosive impact event, which is in line with the physical principle.

1. Introduction

Geophysical mass flows such as landslides and debris flows are often characterized by high flow 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 influence of the factors such as global climate change and engineering activities, the occurrence and frequency of the geophysical mass flow remain high in some areas. The complex material composition and obvious difference in mechanical properties make the geophysical mass flow often show complex and diverse dynamic characteristics [3, 4]. These factors make it more difficult 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]. These interactions may significantly affect the dynamic characteristics and the final disaster scale of the landslide. In this case, reasonable consideration of the obstructive effect of buildings is very important to accurately predict the scale of the disaster.

Concerning interactions between the geophysical mass flow and obstacles, a lot of attention was attracted from researchers [68]. 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 influence of the landslide body on the structure. One-way coupling (flow to buildings) was considered in their study. By using the LS-DYNA software, the effect of building blockage on the mobility of the Shenzhen 2015 landslide and the associated energy dissipation mechanisms are evaluated by Luo et al. [10]. The same landslide has also been studied by the coupled smoothed particle hydrodynamics (SPH) method and finite element method (FEM) [11], an element erosion algorithm is adopted to simulate the destruction process of the buildings.

The dynamic characteristics, such as flow 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. The accurate analysis of these factors has gone beyond the scope of geotechnical engineering issues. The influence of destructible buildings on the dynamic characters and ultimate accumulation extent of the geophysical mass flow is the main focus of this study. In particular, how to carry out relevant numerical simulations and analyses in an accessible way under a unified, advanced physics-based modeling and computational framework.

For water-laden landslides, there exist complex interactions between the solid and fluid phases [12, 13], which influence the dynamic characteristics dramatically [14, 15]. Therefore, the multiphase mass flow model is a more reasonable choice for analyzing landslides with high water content. Pudasaini [13] proposed a general two-phase mass flow model, which includes various interactions between the solid and the fluid constituents. Besides buoyancy, the model also includes generalized drags between phases, enhanced non-Newtonian viscous stress, and virtual mass. These properties make the two-phase mass flow model of Pudasaini [13] and its extensions able to model a wide range of geophysical mass flows [4, 1619].

Numerical modeling is a suitable method for complex geotechnical mass flow and is an important method for back-analyzing disaster events and exploring the mechanisms of disasters [2, 14]. The applicability and effectiveness of the numerical models should be verified by experiments or site events. The tool r.avaflow [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 flow 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.avaflow. These 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.avaflow tool and was utilized for modeling the 2015 Shenzhen landslide (China). The destruction of the buildings caused by landslides is approximated by the erosion process. Through this method, the influences of destructible buildings on the landslide’s dynamic characteristics and the final accumulation state are analyzed.

2. Multiphase Mass Flow Model with Erosion-Entrainment

2.1. Mass and Momentum Conservation Equations

There are some important aspects that should be considered in the analysis of the geophysical mass flow. These include the diverse composition of the material, the complex interactions between the phases, the possible liquefaction resulting from the high-velocity movement, and finally, the variations of the internal components caused by erosion and entrainment along the path and the runout. These aspects are important characteristics and also the source of difficulty in predicting the hazard scale of geophysical mass flows.

Pudasaini and Mergili [4] have proposed a flexible multiphase mass flow model that includes enhanced drag force [21] and virtual mass force models [22]. By considering the differences in mechanics and rheological characteristics between the coarse and fine solid particle and viscous fluid, this multiphase mass flow model can better represent the real geophysical mass flows, 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 fluid phases of the landslide body. The third phase represents the obstacles that can be eroded and entrained into the landslide body. The third phase follows a Coulomb-viscoplastic [25] material behavior, which may have significant viscosity and yield strength. The yield criterion helps distinguish between the flow and nonflow regions of the mass [4]. Therefore, it helps to obtain diverse depositional states of the geophysical mass flow. 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]. The damaged buildings can hardly slide with the landslide over a long distance. Therefore, a solid phase with Coulomb-viscoplastic rheology is suitable for approximating the effects 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. The contribution of obstacle erosion to the mobility of the landslide is affected 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 flow 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, fluid (f) phase of a landslide, and the solid phase of the obstacle (o):where h represents the total flow depth of the mixture, , , and represent solid, fluid, and entrained obstacle volume fraction, respectively. The vectors , , and represent depth-averaged velocities (in the x and y directions) for solid, fluid, and entrained obstacles, respectively. The basal topography is represented by b = b (x, y, t) that may evolve in space and also in time if erosion and entrainment are considered. Es, Ef, and Eo are the solid, fluid, and obstacle erosion rates, which represent the mass productions resulting from the entrainments of solid, fluid, and obstacles, respectively.

Following the pioneering works of Pudasaini and Mergili [4] on the basic momentum conservation equations without considering erosion and entrainment and the works of Pudasaini and Fischer [23] and Pudasaini and Krautblatter [26] on the mechanical models for the mobility of erosive landslides, the momentum conservation equations in the x-direction for the solid, fluid, and entrained obstacle phases are written as follows:where following mechanical principles of Pudasaini and Fischer [23] and Pudasaini and Krautblatter [26], the velocities of the solid, fluid, and obstacle mass that have just been eroded from the bed (in the x and y directions) are denoted as , , and , 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]. Therefore, it is necessary to distinguish the velocity of the entrained mass from that of the landslide mass [26]. Pudasaini and Krautblatter [26] first found and emphasized that the mobility of erosive landslides depends on the ratio of erosion velocity to landslide velocity. , , and are the virtual mass induced mass and momentum enhancements for the solid, fluid, and obstacle-solid phases, respectively. , , and are hydraulic pressure coefficients for solid, fluid, and obstacle-solid phase respectively, and K is the Earth pressure coefficient, is the component of gravitational acceleration in the z-direction (similarly, and are the gravitational acceleration components in the x and y directions, respectively), and and are the density ratios of the corresponding phases. Due to the symmetry of model equations, the y-directional momentum equations have a similar form. The definitions of the source terms and associated parameters on the right-hand side of these momentum equations are presented below.

2.2. The Source Terms and Relevant Parameters

The source terms in x-directional momentum equations can be expressed as [4, 23]

In these equations, is the drag coefficient 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]. The value of the parameter J should be determined based on whether the flow 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 fine particles and fluids as non-Newtonian fluids with a yield. High viscosity and non-Newtonian fluids can have lower mobility than solid particles. Buildings have a specific 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 fluids. Therefore, the high-viscosity, non-Newtonian fluids were utilized to represent buildings impacted by landslides in this study.

The generalized and smooth enhanced drag coefficients that appear on the right-hand sides of (8)–(10) are analytically derived by Pudasaini [21]:

In equation (11a), represents the terminal velocity of a solid particle falling in the fluid, and  =  are parameters representing the proportion of the fluid-like and solid-like drag characteristics in the total drag force, respectively. (0, 1) combines the fluid-like and solid-like drag contributions between the solid and fluid phases. The last item  =  in the denominator is the smoothing or damping function that removes the singularity when the drag coefficient evolves with the proportion of the phase. This new function [21] removes the singularity of the previous version [13] and offers a smooth variation of the drag coefficient as the solid volume fraction evolves. The variable depends on the total mass flux, which is inversely proportional to the drag coefficient. A higher value of means higher mass flux, which in turn leads to lower drag. Although the parameter can be regarded as a numerical parameter, it can also be related to physical quantities such as solid volume fraction [19, 21]. The other two drag coefficients (in equations (11b) and 11c) have a similar form to equations (11a) and reflect the physics-based interaction between the two phases. We refer to Pudasaini and Mergili [4] for the expressions of the other drag coefficients.

The parameter reflects the quantity of mass flux and may vary with the dynamic state of the mass flow [21]. Therefore, dynamic parameters such as flow rate, volume fraction, and shear rate have a potential influence on the value of the parameter . For example, in the accumulation zone, the velocity of the material flow is slowed down, and the material flux should decrease accordingly, which corresponds to an increase in drag. A feasible way to determine the value of is to correlate it with the flow velocity, as suggested by Pudasaini [21]. An expression based on the hyperbolic tangent function is proposed in this study to construct a mass flux that depends on the dynamic state of the landslide. For simplicity, the value of mass flux will be related to the bulk velocity of solid-fluid mixture flows:where the variable represents the bulk velocity of the landslide mixture, and are the upper and lower bounds of the mass flux, respectively. And tanh is the hyperbolic tangent function which returns values close to -1 for and values close to 1 for . The parameter shifts the velocity in which the value changes. For example, when  = 1.0 and  = 2.0, for specified values of and , the curves for the variations of mass flux are demonstrated in Figure 1. It can be found that the parameter m affects 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 flux functions reflect that the changes in only occur within a certain velocity interval. When the velocity is below the lower bound of the interval, the mass flux 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 flux no longer increases, which may correspond to a fully liquefied state. For the models of the virtual mass forces, effective viscosities, pressure-and rate-dependent Coulomb-viscoplastic rheological laws, and non-Newtonian fluid stresses, please see Pudasaini and Mergili [4].

2.3. Erosion-Entrainment Model

Experiments and field observations show that erosion and subsequent entrainment can significantly change the dynamics of geophysical mass flow [2729]. This has been proven by Mergili et al. [18] and Shugar et al. [24] who applied the Pudasaini and Mergili [4] multiphase mass flow 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. This fact is revealed for the first 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. The damaged part is detached from the main part of the building. The detached parts will move with the landslide only after satisfying certain mechanical conditions. This process is similar to the normal erosion process. This 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. This explicitly depends on how the erosion velocity is related to the flow velocity, i.e., ub = λbu, as mechanically proven by Pudasaini and Krautblatter [26]. Treating collapsible buildings as a kind of erodible terrain provides a convenient method for analyzing the influence of buildings on the landslide accumulation area.

The erosion process of the geophysical mass flow 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. The 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. The 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. Therefore, to simulate the influence 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. The variations in the height, shape, and position of the collapsed buildings should also be considered.

This study utilizes the erosion rate formula proposed by Pudasaini and Krautblatter [26] to multiphase mass flow, that is, we treat buildings as erodible solids with Coulomb-viscoplastic characteristics in addition to the coarse-solid phase and fluid phase that make up the landslide mass. The coarse-solid phase was named by Pudasaini and Mergili [4] within their multiphase mass flow 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 flows over the surface with a high velocity (Figure 2). The damage and even collapse of obstacles caused by shearing can be approximated as a kind of erosion. The momentum flux 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 flux [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. Therefore, Eso 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. Therefore, the net shear stress can be expressed as [13, 23]where is the slope angle (Figure 2) and and are the solid mixture volume fractions on either side of the erosion interface, that is, (i = m, b). and are the fluid-solid density ratios of the landslide body and bed, respectively, and are the solid mixture densities of the landslide body and the obstacle, respectively. and can be determined based on the volume fraction relationship between the two solid phases within the landslide and obstacle, respectively, that is, (i = m, b). and are the Coulomb friction coefficients of solid phases within the landslide and the obstacle, respectively, which can also be determined based on solid volume fractions, that is, (i = m, b).

The 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 flow velocity of the solid phase us [23]. With this relationship, the erosion rate can finally be expressed aswhere E0 = Eso has been realized, and is the shear velocity factor that represents the transformation factor between the effective erosional shear stress and the effective velocity jump at the erosion interface. For the fluid erosion rate, see Pudasaini and Fischer [23].

This erosion rate formula is suitable for multiphase mass flow and is mainly based on the mechanically correct and mathematically consistent erosion rate model of the two-phase mass flow [23, 26].

The most outstanding contributions of Pudasaini and Mergili [4] are to extend the two-phase mass flow model to three phases and to take into account the complex interactions between the phases with distinct properties. In particular, the introduction of a fine-solid phase with Coulomb-viscoplastic mechanical behavior enables the model to reflect the more complex and diverse dynamics of the real geophysical mass flow. This 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. Therefore, ground surface erosion caused by landslides is negligible, and only the erosion of the buildings on the runout path will be considered. This is modelled by considering Es = 0, Ef = 0, and Eo is given in equation (15). These expressions are then utilized in equations (1)–(7), modeling the erosion induced mass and momentum productions.

3. 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 artificial landfill [30]. Approximately 2.73 × 106 m3 of construction spoil was mobilized with an impact area of 0.38 km2 in this landslide (Figure 3(b)), and the maximum runout distance is about 1100 m [31]. This landslide destroyed 33 buildings, caused 73 fatalities, and four people missing. The 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 flat terrain and obtains such a high velocity is the main focus of the researchers [15, 19, 3032].

The source area of the landslide was originally a quarry, within which a deep pit was formed due to long-term mining. The 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 landfill for municipal construction waste. The stagnant water at the bottom of the quarry pit was not pumped out before filling with construction waste. This results in higher water content in the landfill soil, especially near the bottom. The filled soil is mainly silty soil and clay with a small amount of gravel grain [31].

Studies have shown that the high-velocity and long-runout 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]. The low permeability of soil results in a significant increase in pore pressure during the sliding [15]. Therefore, when simulating this landslide, the single-phase flow model needs to use a very low basal friction coefficient 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 flow mobility, could not be considered in those simulations. This has been proven in Pudasaini and Krautblatter [26] by presenting a novel mechanical model for landslide mobility with erosion. For the two-phase flow model, it is necessary to consider the high fluid phase ratio and the fluidization 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. Therefore, the influence of buildings is an important factor affecting the accumulation range of the landslide and should be considered reasonable. Due to the complexity of the problem, the influences 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. This 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 influences of obstacles on landslides are approximated through the erosion-entrainment effect.

3.1. Numerical Simulation Tool

In this study, the numerical tool named r.avaflow (V2.4 Pre-release) was utilized. r.avaflow [20] is a flexible and GIS-based open-source computational framework for the simulation of geomorphic mass flows and process chains [17, 18]. This simulation tool is based on the Pudasaini and Mergili [4] multiphase mass flow 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.avaflow, the custom erosion and entrainment model can be defined as a complementary function for the source items. The maximum depth of entrainment can be defined 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 flow topography simultaneously. At the same time, flow momentum is updated accordingly.

3.2. Model Parameters of 2015 Shenzhen Landslide

The hydraulic pressure coefficient in equation (5) shows that buoyancy reduces the solid normal and lateral stresses by a factor of . Therefore, a greater buoyancy (pore fluid pressure) effect can be reflected 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/m3 and ρf = 1100 kg/m3 are used for the solid and fluid phases of the landslide body. Obstacles have a higher density. The equivalent bulk density of the buildings cannot be equivalent to the density of concrete because there are large empty spaces inside the buildings. Therefore, a density value higher than the landslide density and lower than the concrete density is assigned to the obstacle, that is, ρo = 2000 kg/m3. Based on previous field 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 affecting the interaction between multiple phases. Parameter P has a negative effect on the advection and dispersion of the flow [13]. It can be associated with specific 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 flux contribution as . For natural debris flows or fluidized 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 can be estimated as 10 ms−1. The default value of the flux parameter is 1.0 in r.avaflow. For the landslides with high-water content in this study, to reflect the influence of the landslide motion state on the drag force, the flux parameter in the form of equation (12) was adopted, and  = 1.0 and  = 2.0 were used.

In Figure 4, three drag coefficient curves with constant mass flux are compared following the analytical, smooth, and enhanced drag model developed by Pudasaini [21]. The difference between the drag coefficient curves of  = 1 and  = 2 gradually increases as the solid volume fraction falls within the usual range for the natural geophysical flow, i.e., 0.4 < αs < 0.6. The drag coefficient curve with  = 0 corresponds to the old drag function of Pudasaini [13], which rapidly tends to an infinite value as the solid phase volume fraction gradually increases. The new drag function [21] 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 significantly. The drag feature of this dense regime can be described by the new drag function [21]. In this study, the drag curves with  = 1 and  = 2 can be regarded as the upper and lower bounds of the drag functions with velocity-dependent mass fluxes. When the solid phase volume fraction is in the normal interval (0.4 < αs < 0.6), the variation in drag coefficient caused by the velocity factor is relatively small. Equation (12) can be regarded as a conservative approximation of the actual drag coefficient.

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–20). The relationship between the pore water pressure gradient and mass flux can be approximately regarded as a linear relationship [1]. Therefore, it can be approximately considered that in the interval (10–20), with the increase in velocity, the mass flux will increase to a certain extent. So, the middle value of the velocity interval is estimated to be 15 m/s, and m = 3 is adopted, which means the velocity transition curve is wider and flatter (Figure 1). To replace the original constant mass flux with equation (12), the source code of r.avaflow should be modified accordingly.

The soil in the source area has a high content of fine particles. The 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. This 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 fine 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 fine 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 fine solid particles in the landslide body, are set as 31° and 18°, respectively. For fluid with a high content of fine sediment, the viscosity ηf = 2 Pa.s was adopted in the simulation. The terminal velocity UT is 0.1, and the particle Reynolds number Rep is 1. The exponent J for drag in equation (11) takes the value 1 for laminar-type flow. The viscosity ηo = 102 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 verified [8, 13, 17, 18, 26].

Based on the field 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. The 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. The final average accumulation thickness of the landslide is about 8 m, and the thickness of the landslide in flow is smaller than this value. Therefore, it can be considered that the distribution of landslide velocity in height conforms to the hypothesis of plug flow, that is, . The erosion drift coefficient 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, . The shear velocity factor satisfies  = 0.05, it describes the rate at which the eroded mass is entrained into the flowing mass, or depositional mass is removed from the flowing mass. The value of this parameter follows the recommendation of the software manual. The basal friction angle of the bed is 8°, which is 10° smaller 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.

3.3. 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]. Therefore, once the landslide is triggered, it will attain a higher velocity soon. Due to the high-water content of the actual soil, a higher fluid volume fraction is given to the soil of the source area in the simulations. For a multiphase flow model, due to the high initial volume fraction of the fluid phase, a significant increase in buoyancy and a reduction in friction will be obtained. Therefore, once the simulation starts, the mass in the landslide source area will displace immediately. To compare the influence of buildings on the dynamic characteristics of the landslide, the buildings within the influence range of the landslide are regarded as erodible masses with Coulomb-viscoplastic characteristics in simulation (I). In simulation (II), the influence 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.avaflow, both erodible and nonerodible buildings are a part of the initial terrain. Through the numerical simulation of r.avaflow, the impact range and dynamic characteristics of the landslide at different times are obtained.

3.4. Accumulation Range of the Landslide

In simulation (I), after the triggering, the landslide mass in the source area flowed out from the low-lying region on the north side of the landfill. The 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 sliding mass was further hindered by the buildings on the sliding front. The 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 first impedance from the buildings occurred on the west side of this dashed line, and the subsequent impedance occurred eastward along this dashed line. The flow on the east side of the landslide was hindered last, so this part of the flow had the longest runout distance. Due to the obstruction of the buildings, the flow 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. The 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. The entrained obstacle phase has a smaller volume fraction in the local total flow 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. The amount and velocity of the landslide mass that subsequently flowed into the accumulation area decreased significantly (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. The 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. The eroded obstacle phase slides away with the landslide only a little. This 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. The change in the flow direction of the landslide occurred later than in simulation (I). So, the landslide achieved significantly larger runout distances than in simulation (I) at the same time. This phenomenon is especially obvious on the right side of the landslide front (see Figures 6 and 7). The 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 flow characteristics of the landslide differ from reality due to the absence of this obstruction. By comparing the actual flow 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. Therefore, simulation (II) and simulation (I) have the same maximum sliding velocity, which is 30.1 m/s. This value is very close to the maximum velocity value estimated by Yin et al. [31], which is 29.8 m/s.

3.5. Evolution of Erosion Characteristics

For simulation (I), the evolution characteristics of the erosion extent during the erosion development stage are shown in Figure 8. The buildings within the impact range can be divided into three groups by the yellow dashed lines (Figure 8). The 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 flowing mass show a significant change in flow directions. This change is consistent with the shape of the final accumulation (represented by the solid red line in Figure 8).

The evolution characteristics of erosion can be further compared in the profile section (yellow line in Figure 6 or Figure 3(b)). In Figure 9, the flow height is denoted by ho 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 zb. In the early stage (Figure 9), erosion develops progressively along the direction of the landslide movement. The 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 ho 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

The velocity characteristics of the solid-phase along the profile 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. The 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. The erosion process of the obstacle is determined by the relationship between the net stress and the net momentum flux on both sides of the erosion interface. The 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 influence 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. The 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.

The red curve in Figure 11 represents the evolution of the erosion volume (left vertical axis). The 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 6). As the landslide progresses, the landslide mass eroded a large number of buildings in front of it, and the eroded volume increases significantly. The erosion volume reaches its maximum value at t = 90 s and then decreases slowly. The decrease in eroded volume is caused by the accumulation of the obstacle phase. The maximum total erosion volume is about 50 × 104 m3, which is about 18.3% of the source material volume, which is substantial for the situation considered here with mobilized buildings. Therefore, the accumulation range and morphology in simulation (II), which does not consider this volume, are quite different from the actual ones.

4. Discussion

Geophysical mass flows in nature often encounter various natural or artificial obstacles on their movement path, such as houses, retaining dams, and diversion dikes. When analyzing the dynamic characteristics of the geophysical mass flow, artificial 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 flow becomes complicated. In this case, a numerical tool based on a unified computing architecture becomes necessary. This complex analysis is made possible by the multiphase geophysical mass flow model [4], the mechanical erosion rate model [23] and the erosive landslide mobility theory [26], and the associated numerical implementation in r.avaflow.

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 flux on both sides of the erosion interface. The erosion drift coefficient is determined based on the physical quantities on both sides of the erosion interface, i.e., . Since , , and , then . As has been demonstrated by Pudasaini and Krautblatter [26], when , the landslide loses kinetic energy during the erosion process, that is, the obstacle hinders the sliding of the landslide. This demonstrates the functionality of the state of the erosion-induced mobility presented by Pudasaini and Krautblatter [26]. This clearly manifests the practical and engineering importance of this research work that practitioners can take direct benefit 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 profile section with different yield strengths. It is observed that the depositional characteristics of eroded obstacles are affected by their yield strength. The 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.

5. Conclusions

The drag effect between the solid and liquid phases is an important aspect that affects the dynamic characteristics of the multiphase mass flow. To consider the evolution of the drag coefficient [21] within a landslide, a new formula for the mass flux parameter for Pudasaini’s drag coefficient was proposed. This new formula reflects the influence of the landslide velocity on mass flux parameters. Therefore, 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 flow. The 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. The 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 coefficient of the entrained buildings should be smaller than that of the landslide body. The entrainment phenomenon has been explained by the theory of erosive landslide mobility of Pudasaini and Krautblatter [26].

The multiphase mass flow model [4] which is applied to a GIS-based computational tool r.avaflow 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 flow mechanical characteristics and an accumulation range of geophysical mass flow through effective numerical tools. This study provides a method to approximately consider the influence of erodible obstacles such as buildings and constitutes a useful reference for the prediction and assessment of geophysical disasters associated with landslides and debris flows.

Data Availability

The data supporting this study are available upon request from the corresponding author.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


This study was supported by Anhui Provincial Natural Science Foundation (grant no. 2108085ME191), the Natural Science Foundation of the Higher Education Institutions of Anhui Province (grant no. KJ2019A0128), the Talent Introduction Foundation of Anhui University of Science and Technology (grant no. 12062), the Youth Foundation Project of Anhui University of Science and Technology (grant no. 11887), and the Graduate Innovation Foundation of Anhui University of Science and Technology (grant no. 2021CX2045).