# Numerical Investigation of Impact Characteristics of Rigid Retaining Wall for Blocking Landslides: A Case Study of the Shum Wan Road Landslide

Authors: Yunyun Fan; Feng Zhang; Ying Liu; Dongyang Ma
Journal: Geofluids
 (2022)

---

## Abstract

Landslides pose a serious threat to the downstream infrastructure and personnel safety. For situations where the location of possible landslides is difficult to determine, or the potential landslides areas are widely distributed and scattered, passive protective structural measures need to be adopted to relieve the threat of landslides. Taking the Shum Wan Road landslide as an example, the discrete element method was used to reproduce and analyze the landslide movement. On this basis, this study conducts a further numerical investigation on the impact characteristics of the rigid retaining wall used to block landslides. The results show that the numerical back analysis can well represent the evolution of the landslide, and the numerical results are in good agreement with the actual situation, thus confirming its correctness and effectiveness. Under the condition of setting a rigid retaining wall, the landslide can be effectively blocked and form a final deposition behind the retaining wall. The impact on the rigid retaining wall varies with time, and the distribution of the impact in space is also uneven. Through the statistical analysis of the contact forces on the retaining wall, the distributions of static and peak stresses on the retaining wall can be obtained. Compared with the static stability, the spatial distribution of dynamic impact is more difficult to predict, especially the area of the peak impact has a significant influence on the design. Numerical investigation can efficiently and reliably provide abundant information of static stress and dynamic impact on the rigid retaining wall, which can provide important reference for the optimal design of the rigid retaining wall structure.


---

## Body

## 1. Introduction

Landslide is one of common natural geological disasters in mountainous areas, which often leads to damage to infrastructure and poses a serious threat to people’s lives and property safety [1, 2]. In practical engineering, when it is difficult to adopt accurate and comprehensive active measures, passive measures of blocking or evacuating are needed to deal with them [3, 4]. The engineering situations in which passive measures are applied generally include the following: the location of potentially dangerous slopes is difficult to determine, or the potentially dangerous areas with insufficient stability are widely distributed and scattered. Passive measures can protect infrastructure by blocking or changing the path of landslide movement [3]. Common passive protection measures include shed galleries [5], baffles [6], passive nets, and retaining walls [7–9]. Among them, retaining wall is a key structure that can provide reliable protection downstream.Due to the limited understanding of the process of landslide impacting protection structures [10], a general standard for the design of protective structures has not been established [11, 12], so the design of passive protection structures has to adopt empirical methods [6]. However, landslides often have high impact forces [13], and empirical methods often fail to meet practical safety requirements [4]. Unreasonable layout schemes often lead to structural damage or even failure accidents, which may cause more serious damage and losses [14]. In-depth research on the movement of landslides under control and their impact on protective structures is an important topic that must be faced in prevention and mitigation projects [12, 15].The landslide shows the characteristics of fluid during the movement, and its constitutive relation fits with the solid, so the landslide movement has the characteristics of both fluid and solid [16]. Landslides are controlled by different mechanisms in different movement stages. At low velocity, the landslide movement is dominated by friction [17], while as the velocity increases, the collision becomes progressively more pronounced [12, 18, 19]. Under the blocking action of the protective structure, the motion state of the landslide is changed [6], and its energy is also consumed [20]. In order to reveal the mechanism of impact and deposition in the process of landslide impacting protective structures, scholars have made further studies on the dynamic characteristics of the process [12], among which, the interaction between landslides and rigid retaining walls has received the most extensive attention [21]. The studies found that the impact force on the protective structure is related to the scale, impact angle, and velocity of the landslide [22]. These factors affecting the impact force are not only influenced by terrain conditions such as slope angle, slope friction, and sliding surface [22, 23] but also controlled by the composition and mechanical properties of geomaterials [23]. Some preliminary estimation methods have been proposed based on experiments and analysis [22]. The shape of the protective structure can impose different constraints on the landslide, which significantly affect the dissipation of energy, thereby achieving different protective effects [4]. Blocked by the protective structure, the deposition behind the retaining wall will form a dead zone [24], in which the mass no longer moves. Then, the mass in this zone changes the geomorphological conditions of the subsequent mass, which will make the subsequent mass run-up in the vertical direction [25]. Run-up is one of the important reasons for the peak impact force [18, 26]. Studies have shown that particle size and gradation will have a great influence on the run-up height and the peak impact force [11]. In addition to retaining structures, some open structures can slow down the landslide, promote energy dissipation, and change the trajectory of the landslide [6]. If the site allows, open multilayer structures are more effective for disaster mitigation [27]. The open mitigation structure can effectively reduce the landslide velocity and retain some geomaterials, reducing the impact on the downstream retaining structure and also bringing more optimization possibilities for the systematic protection configuration [4, 12, 28]. The interaction between landslides and protective structures is an important aspect of disaster protection [29]. However, this interaction remains to be further studied so far [30]. Obviously, the revealed key influencing mechanisms are very helpful for understanding the impact characteristics of retaining walls. These mechanisms can be used to comprehensively and quantitatively assess the forces acting on retaining walls, as well as their ability to block landslides and regulate deposition [12].Currently, available research methods include field test [31], laboratory test [32], theoretical analysis [33], and numerical simulation [27]. Field test and laboratory test are the basis for theoretical and numerical studies [34]. However, research is sometimes limited by experimental conditions and test equipment, and analytical solutions cannot face complex computational conditions, so numerical simulation becomes the most reliable means [18]. In the context of the rapid development of modern numerical theory and technology, numerical methods have almost become one of the main means to study landslide disasters due to their advantages in computational cost [9]. Numerical methods currently used to simulate landslide movement mainly include finite difference method (FDM) [30], finite volume method (FEM) [35], lattice Boltzmann method (LBM) [8], discrete element method (DEM) [27, 36, 37], smoothed particle hydrodynamics (SPH) [38], and material point method (MPM) [28, 39]. These numerical methods are mainly based on two types of theories, continuum and discrete media. The model of average integration along the depth derived from the continuum theory can simulate large-scale landslides [34] and has been widely used [30, 35]. This type of model can match different bottom resistances according to the geomaterial composition and properties of the landslide [35, 40]. At the same time, lateral pressure coefficients were introduced into the model to reflect the different motion states of geomaterials [24, 40]. Compared with the traditional fluid model, the bottom surface resistance and lateral pressure coefficient combined with the geomaterial properties are the distinguishing features of this model. Depth-averaged theory models have achieved great success in engineering applications, but some studies consider them inappropriate for landslide-structure interactions [4, 10]. This is because the model cannot take into account the vertical variation of physical quantities, thus limiting the ability to obtain some key information [41, 42]. The discrete element method proposed by Dr. Cundall is a representative method based on discrete medium theory. Using this method, it is convenient to investigate the impact characteristics of rigid walls and the energy dissipation process [6, 43], which can be used to evaluate the rationality of protective structures [32, 44] and discuss the possibility of optimal design of protective structures [12, 43]. The high configuration requirement of the computer is the main disadvantage of this method, but this has not affected the enthusiasm of scholars to use it to study the mechanism [15]. Compared with the methods based on continuum theory, the discrete element method has irreplaceable advantages in simulating particle morphology and gradation, as well as mesomechanisms.At present, some phenomena and mechanisms in the process of blocking landslides by protective structures have been obtained, which are meaningful for designing more reasonable protective structures. However, due to the complex environment in which protective structures are applied, the existing results still cannot fully meet the needs of the design and practice of protective structures [45]. Numerical investigation work will help improve the design of protective structures [12]. Therefore, it is necessary to conduct a more in-depth and extensive investigation and research on the impact acting on protective structures. In this study, the PFC3D (Particle Flow Code 3D) software based on the discrete element method is used to numerically investigate the impact on rigid retaining wall structures and to discuss how to use the numerical results to provide guidance for the practical optimization of protective structures.


## 2. Case Study

### 2.1. Numerical Method

In this study, the PFC3D software (version 3.10-258) based on the discrete element method is used to numerically simulate the movement and impact characteristics of the landslide. In the calculation of each time step of the discrete element method, it is first necessary to judge whether the discrete elements are in contact, and then, use the contact model equation to calculate the contact force between the discrete elements, and finally, calculate the velocity and position of each discrete element by using Newton’s second law (equation of motion) according to the force and moment acting on the discrete element. The contact and motion information of discrete elements are continuously updated, through cyclic calculation, and then, the dynamic solution of the studied object is realized.The contact equation describes the relationship between the relative displacement between elements and the contact force. The contact force between discrete elements can be decomposed into two components, the normal force and the shear force, which are perpendicular and parallel to the contact plane, respectively. They can be expressed by(1)Fin=knun,∆Fis=−ks∆uswhere Fin and ∆Fis are the normal contact force and shear force increment between discrete elements, respectively, kn and ks are the normal and shear contact stiffness, respectively, and un and ∆us are the normal displacement and relative shear displacement increment, respectively.Equation (1) is the most commonly used linear stiffness contact model in the discrete element method. The normal contact stiffness characterizes the relationship between the cumulative normal contact force and the cumulative normal displacement. The shear stiffness characterizes the relationship between the shear contact force and the shear displacement increment.In the linear stiffness model, the normal stiffness and shear stiffness between discrete elements can be expressed in Equation (2) as:
(2)kn=knAknBknA+knB,ks=ksAksBksA+ksB,where the superscripts A and B, respectively, represent two discrete elements in contact with each other.According to the shear force and the normal force, it can be determined whether the two contact discrete elements will slide, and the friction force at the contact is defined by the friction coefficientμ, which is equal to the minimum friction coefficient of the two contact discrete elements, that is, μ=minμA,μB. When the shear contact force Fis at the contact point is larger than the maximum static friction force Fsmax=μFin, there will be relative sliding between the two contact elements, and the sliding friction force at the contact point is Fis=minμA,μBFin.According to Newton’s second law, the translation and rotation of the discrete element need to be calculated at each time step. The motion equation for this purpose is shown in(3)Fi=mix¨i+βgẋi,Mi=Iiφ¨i+βgφ̇i,where i=1,2,3 represents the components along the three coordinate axes; mi is the mass of the discrete element, Fi is the unbalanced force component acting on the discrete element, Mi is the unbalanced moment component acting on the discrete element; xi is the translational displacement of the discrete element, φi is the angular displacement, Ii is the rotational inertia, and βg is the damping coefficient.Using the difference method, the relationship between translational and rotational variables can be obtained, as shown in(4)x¨it=1∆tẋit+∆t/2−ẋit−∆t/2,ẋit=12ẋit+∆t/2+ẋit−∆t/2,(5)φ¨it=1∆tφ̇it+∆t/2−φ̇it−∆t/2,φ̇it=12φ̇it+∆t/2+φ̇it−∆t/2.Using the motion Equation (3) and the relationship between the amounts of motion Equations (4) and (5), the velocity equations of translation and rotation can be obtained:
(6)ẋit+∆t/2=1−1/2∙βg/mi∙∆tẋit−∆t/2+Fit∙∆t/mi1+1/2∙βg/mi∙∆t,φ̇it+∆t/2=1−1/2∙βg/Ii∙∆tφ̇it−∆t/2+Mit∙∆t/Ii1+1/2∙βg/Ii∙∆t.At this time, the translational displacement and rotational displacement can be expressed as:(7)xit+∆t=xit+ẋit+∆t/2∙∆t,φit+∆t=φit+φit+∆t/2∙∆t.After performing the above iterative calculations for all discrete elements, the dynamic simulation results of the studied object are presented.


### 2.2. Shum Wan Road Landslide

On August 13, 1995, a landslide took place at the hillside above Shum Wan Road, Hong Kong. It caused the collapse of a 30 m long section of Nam Long Shan Road, as shown in Figure1. The literatures reported that the landslide debris crossed Shum Wan Road and damaged three shipyards and a factory near the seafront [46–48]. Prior to the landslide, the hillside was densely vegetated and had an overall gradient of about 27°. The geology at the landslide area comprised a thin mantle of colluvium overlying partially weathered fine-ash to coarse-ash crystal tuff. The landslide resulted in a 70 m high scar, with a width varying from about 50 m just below Nam Long Shan Road to about 90 m above Shum Wan Road. Figure 2 shows the map of the landslide based on topographic survey, geological mapping, and field observations. The upper part of the landslide surface was concave in shape and was up to about 12 m in depth below the prefailure ground surface. The landslide released about 2.6∗104m3 of soil and rock. The debris was deposited on Shum Wan Road and the reclaimed land, spreading over an area of about 0.5 ha. The failure was caused principally by the presence of weak layers in the ground, ingress of water during prolonged heavy rainfall, and partial blockage of the drainage system.Figure 1
The Shum Wan Road landslide on August 13 1995 [46].Figure 2
The map of the Sum Wan Road landslide (modified from Knill 1996).


### 2.3. Back Analysis of the Landslide

The terrain, sliding surface, and disaster scope of the Shum Wan Road landslide are given in the literature [46–48]. A geometric terrain model was established according to the contour of the terrain where the Shum Wan Road landslide is located, and the terrain was discretized by triangular wall elements, as shown in Figure 3.Figure 3
Numerical model of the Shum Wan Road landslide.To obtain a discrete element model for simulating the landslide, it is necessary to use ball element to fill the landslide body. Considering the simulation effect and calculation cost, a ball element with a radius of 0.75 m was selected for the numerical simulation according to the scale of the landslide. The outer surface of the landslide is determined by the ground surface before failure and sliding surface, which are discretized by wall elements, thus forming an enclosed space. Ball elements are filled in the enclosed space formed by wall elements (which is equivalent to the initial landslide) and reach equilibrium under gravity. Finally, a numerical model with discrete elements was established as shown in Figure3.According to the information given in the literature [46–48], the friction coefficient between the ball and the wall can be determined, and the density of the ball can be determined through mass equivalent [43, 49]. However, the stiffness and damping coefficient of discrete elements need to be obtained through numerical back analysis. As a result, the mesocalculation parameters of discrete elements that conform to the landslide movement are shown in Table 1.Table 1
Mesoparameters used in the numerical simulation of Shum Wan Road landslide.
Ball elementWall elementBall diameter (m)Ball density (kg/m3)Normal stiffness (N/m)Shear stiffness (N/m)Friction coefficientNormal stiffness (N/m)Shear stiffness (N/m)Friction coefficient1.5026501e101e100.2141e101e100.214When the landslide starts to move, delete the wall elements that act as a constraint on the upper surface of the landslide body (these wall elements are only used to establish a numerical model of the landslide), and then, the ball elements used to simulate the unstable geomaterial start to move down the slope under gravity. Through numerical simulation, the movement process of the Shum Wan Road landslide was obtained, as shown in Figure4.Figure 4
Flow depth contour map of the Shum Wan Road landslide (m).
(a)
Initial moment of the landslide movement(b)
Numerical result at the 5th second(c)
Numerical result at the 10th second(d)
Numerical result of the final depositionIt can be seen from Figure4 that when the slope enters failure state, the landslide begins to move downward along the slope (Figure 4(a)). The movement is mainly restricted and affected by the terrain. The landslide is mainly accelerated in the early stage of movement, the front edge moves faster, and the rear edge moves relatively slowly, and the largest flow depth area is still at the back of the landslide (Figure 4(b)). The movement process of the landslide developed rapidly. As the landslide entered a gentle section, the sliding force it received was reduced, and the bottom resistance increased, so the front edge of the landslide started to decelerate. The landslide almost reached the actual disaster scope within 10 s (Figure 4(c)), and the final deposition was basically formed within 50 s (Figure 4(d)). During this process, the main mass of the landslide went downstream, and most of the mass deposited at the foot of the hillside, thus causing great harm to the road and the infrastructures downstream.In the numerical results shown in Figure4, the landslide movement is generally within the actual disaster scope, and the edge of the landslide is consistent with the actual disaster scope. However, the disaster scope obtained by numerical simulation in the downstream area full of buildings and facilities is different from the actual situation to a certain extent. The reason for this difference is the existence of these buildings and facilities. The buildings and facilities downstream of Shum Wan Road were destroyed by the landslide, but in fact, they also played a role in blocking the landslide.But it is a pity that due to the lack of information on the distribution of buildings and facilities before the occurrence of the landslide as well as their structure and materials, their blocking effects on the landslide cannot be considered in the calculation (even if the information is mastered, it is still extremely difficult to realize a numerical simulation that carefully considers the coupling effect of landslide impact and structural collapse). It can be seen from the numerical results that although the blocking effect of buildings and facilities is not considered, the numerical simulation has generally given a true representation of the movement process of the Shum Wan Road landslide. The numerical results obtained are correct and reliable, and the numerical calculations are also effective and credible. The huge damage to the downstream infrastructure caused by the Shum Wan Road landslide reminds us that it is necessary to make further specific investigations into the impact of the landslide.


### 2.4. Impact Characteristics of Rigid Retaining Wall

When the position of the unstable slope cannot be accurately determined, or the unstable slope is widely distributed, reinforcing the slope costs too much but cannot ensure an effective protection of objects. In this case, it can also be considered to protect the objects by setting a rigid retaining wall upstream.The Shum Wan Road Landslide is a typical roadside landslide. There are many infrastructures in the lower reaches of Shum Wan Road, and the employees working there are also threatened and affected by the landslide. In the environment of Shum Wan Road Landslide, rigid retaining wall is a representative protective structure. Moreover, according to the characteristics and scale of the landslide, it is appropriate to use a rigid retaining wall as a protective structure. Therefore, the Shum Wan Road Landslide was selected as a case to study the impact characteristics of rigid retaining walls.The completed numerical back analysis has led to a deeper understanding of the landslide movement process. Numerical back analysis also makes it possible to further study the impact of landslides on roads and to in-depth study of the mechanical characteristics of blocking structures along roads. Numerical investigations are important basis for the design of protective structures.To understand the impact of landslides on the road, vertical rigid retaining walls are set along the road to block the landslide. The main function of these rigid retaining walls is to block the moving mass to ensure the safety of downstream infrastructures and personnel. The setting of the rigid retaining wall is shown in Figure5.Figure 5
Rigid retaining wall along the road.The numerical investigation of the landslide process under the condition of setting a rigid retaining wall is carried out again, in which the rigid retaining wall adopts the parameters used for the wall elements in Table1. The impact on the rigid retaining wall was recorded, and the landslide movement process and the time histories of the impact were obtained as shown in Figures 6 and 7.Figure 6
Flow depth contour map of landslide process blocked by rigid retaining wall.
(a)
Initial moment of the landslide movement(b)
Numerical result at the 7th second(c)
Numerical result at the 15th second(d)
Numerical result of the final depositionFigure 7
The time histories of the impact.
(a)
The time history of the force (N)(b)
The time history of the moment (N∗m)From the results shown in Figure6, it can be seen that the landslide moves along the slope surface under the gravity after failure (Figure 6(a)). The landslide approached the road in about 7 s (Figure 6(b)). After that, blocked by the rigid retaining wall, the process of the landslide showed an obvious difference from the process without the rigid retaining wall. The moving mass forms a preliminary deposition under the blocking of the rigid retaining wall (Figure 6(c)), and then, the subsequent mass can be affected by the deposition formed before, expand along the formed deposition in the direction of the road, and gradually form the final deposition (Figure 6(d)). The time histories of the impact acting on the retaining wall shown in Figure 7 correspond to the process of blocking the landslide. The numerical method used in this investigation is the discrete element method, and the impact curve of the rigid retaining wall has certain volatility. However, the changing trends of the force and moment on the rigid retaining wall along the road are clear, and the trends of the two time-history curves are in good agreement, and they have a good correspondence with the dynamic and static loads of the landslide. Due to the arrival of the landslide, the rigid retaining wall began to be impacted at 7.4 s. At 27.1 s, the deposition of geotechnical materials was basically completed, and the force also tended to be stable. During this period, the force on the rigid retaining wall fluctuated due to the high-frequency collision between particles, and this fluctuation also exists objectively in the actual impact process. The time history curves of the impact on the rigid retaining wall give an important hint, that is, for the blocking structure, we should not only pay attention to the static force of the landslide stable deposition on it, but also cannot ignore the impact caused by the landslide. In order to change the motion state of the landslide, the impact peak value that the rigid retaining wall should bear is an important factor that must be considered in the engineering design. Since an obvious run-up phenomenon has not been found in the numerical process of this case, the peaks of the two impact time-history curves of the rigid retaining wall are generated by the dynamic impact of the landslide. The time history curves of the impact on the structure are an important reference data for structural design, and numerical simulation is a convenient and reliable way to obtain these data.After the landslide stopped moving, a stable deposition was formed, which produced a stable pressure on the rigid retaining wall. Compared with the dynamic impact during the landslide movement, the final pressure of the deposition on the rigid retaining wall is long-term and is the primary factor that needs to be paid attention to. In Figure6(d), it can be seen that the deposition with high middle and low sides is finally formed behind the rigid retaining wall, so the distribution of the force on the rigid retaining wall along the road direction is uneven. After the deposition became stable, the discrete elements formed a contact force distribution on the rigid retaining wall, as shown in Figure 8.Figure 8
Contact force distribution on rigid retaining wall in stable state (N).From the contact force distribution shown in Figure8, the outline of the contact distribution is basically similar to the shape of the deposition. The contact force on the retaining wall is discrete, and these discrete contact forces are caused by the contact between discrete elements and the wall. Although the distribution of the position and magnitude of the contact force on the retaining wall looks a little messy under the transmission of the force chain, these discrete forces still form certain regularity as a whole. For example, the overall law of contact force is that the bottom of the middle area is more stressed, and the edge of the deposition is less stressed. In the example of this study, although the contact force on the rigid retaining wall is mainly in the horizontal direction, it is formed by force transmission under gravity, so the magnitude of the contact force is closely related to the local height of the deposition.Since the numerical method divides the landslide into many discrete elements, the simulation is inconsistent with the actual continuous state of the landslide. Although numerical simulations yield discrete contact force distributions on rigid retaining walls, they cannot be directly applied to the design of protective structures. These contact forces need to be counted in an equivalent way, so as to obtain the results required by the design, which are more helpful to engineering applications. The stable force and moment on each section along the road are counted in sections, and the results are shown in Figure9.Figure 9
Stable load distributions on rigid retaining wall along the road.
(a)
Distribution of the force (N)(b)
Distribution of the moment (N∗m)The deposition generally has the characteristics of high in the middle and low on both sides. The distributions of the force and moment on the rigid retaining wall along the road direction are uneven, but their distribution laws are similar. These distributions are mainly controlled by the deposition height, and the distribution characteristics of stable force on the retaining wall along the road are fully demonstrated in Figure9. Compared with the force on each section, the moment is more significantly affected by the height. This is because the soil pressure gradually increases in the vertical direction, so that the moment on each rigid retaining wall is generally in proportion to the square of the deposition height.Further, it is more urgent to obtain a continuous stress distribution that is closer to the actual situation, so the contact force on the rigid retaining wall in the stable state is converted into the stress distribution by means of statistical averaging, as shown in Figure10.Figure 10
The normal stress distribution on rigid retaining wall in stable state (Pa).From the stress distribution on the rigid retaining wall in the stable state shown in Figure10, it can be seen that the law of stress distribution and contact force is similar. The stress is larger at the bottom of the middle area and smaller at the edges of both sides. Compared with the contact force distribution results that cannot be directly applied to engineering, the stress distribution provides more abundant information and can provide a basis for optimized protective structure design. In addition, in order to investigate the nonnegligible dynamic impact, the historical maximum stress on the retaining wall was recorded and counted, and the maximum normal stress distribution in history was obtained, as shown in Figure 11.Figure 11
The maximum normal stress distribution on the rigid retaining wall in history (Pa).Comparing Figures10 and 11, it can be seen that the impact stress on the rigid retaining wall is much larger than that in the stable state. The price of the rigid retaining wall forcing the landslide to stop in a short time is that it has to withstand a larger impact. Interestingly, there is a high stress area in the center of the retaining wall, which is the impact area of geomaterials with high velocity. It is difficult to determine these impact areas and the impact stress on them without numerical methods. Combined with the distribution of the impact area, it can be considered to paste some soft materials for energy dissipation on the upstream surface of the retaining wall, which can reduce the peak of the impact force by prolonging the impact load time. Of course, there may be some differences between the numerical calculation and the actual situation, which allows the scope of mitigation materials to be appropriately enlarged. Under complex terrain conditions, numerical simulations can be used to obtain the impact characteristics of the protective structure that are consistent with the actual situation. The results of these impact characteristics can provide great help for the design of protective structures.


## 2.1. Numerical Method

In this study, the PFC3D software (version 3.10-258) based on the discrete element method is used to numerically simulate the movement and impact characteristics of the landslide. In the calculation of each time step of the discrete element method, it is first necessary to judge whether the discrete elements are in contact, and then, use the contact model equation to calculate the contact force between the discrete elements, and finally, calculate the velocity and position of each discrete element by using Newton’s second law (equation of motion) according to the force and moment acting on the discrete element. The contact and motion information of discrete elements are continuously updated, through cyclic calculation, and then, the dynamic solution of the studied object is realized.The contact equation describes the relationship between the relative displacement between elements and the contact force. The contact force between discrete elements can be decomposed into two components, the normal force and the shear force, which are perpendicular and parallel to the contact plane, respectively. They can be expressed by(1)Fin=knun,∆Fis=−ks∆uswhere Fin and ∆Fis are the normal contact force and shear force increment between discrete elements, respectively, kn and ks are the normal and shear contact stiffness, respectively, and un and ∆us are the normal displacement and relative shear displacement increment, respectively.Equation (1) is the most commonly used linear stiffness contact model in the discrete element method. The normal contact stiffness characterizes the relationship between the cumulative normal contact force and the cumulative normal displacement. The shear stiffness characterizes the relationship between the shear contact force and the shear displacement increment.In the linear stiffness model, the normal stiffness and shear stiffness between discrete elements can be expressed in Equation (2) as:
(2)kn=knAknBknA+knB,ks=ksAksBksA+ksB,where the superscripts A and B, respectively, represent two discrete elements in contact with each other.According to the shear force and the normal force, it can be determined whether the two contact discrete elements will slide, and the friction force at the contact is defined by the friction coefficientμ, which is equal to the minimum friction coefficient of the two contact discrete elements, that is, μ=minμA,μB. When the shear contact force Fis at the contact point is larger than the maximum static friction force Fsmax=μFin, there will be relative sliding between the two contact elements, and the sliding friction force at the contact point is Fis=minμA,μBFin.According to Newton’s second law, the translation and rotation of the discrete element need to be calculated at each time step. The motion equation for this purpose is shown in(3)Fi=mix¨i+βgẋi,Mi=Iiφ¨i+βgφ̇i,where i=1,2,3 represents the components along the three coordinate axes; mi is the mass of the discrete element, Fi is the unbalanced force component acting on the discrete element, Mi is the unbalanced moment component acting on the discrete element; xi is the translational displacement of the discrete element, φi is the angular displacement, Ii is the rotational inertia, and βg is the damping coefficient.Using the difference method, the relationship between translational and rotational variables can be obtained, as shown in(4)x¨it=1∆tẋit+∆t/2−ẋit−∆t/2,ẋit=12ẋit+∆t/2+ẋit−∆t/2,(5)φ¨it=1∆tφ̇it+∆t/2−φ̇it−∆t/2,φ̇it=12φ̇it+∆t/2+φ̇it−∆t/2.Using the motion Equation (3) and the relationship between the amounts of motion Equations (4) and (5), the velocity equations of translation and rotation can be obtained:
(6)ẋit+∆t/2=1−1/2∙βg/mi∙∆tẋit−∆t/2+Fit∙∆t/mi1+1/2∙βg/mi∙∆t,φ̇it+∆t/2=1−1/2∙βg/Ii∙∆tφ̇it−∆t/2+Mit∙∆t/Ii1+1/2∙βg/Ii∙∆t.At this time, the translational displacement and rotational displacement can be expressed as:(7)xit+∆t=xit+ẋit+∆t/2∙∆t,φit+∆t=φit+φit+∆t/2∙∆t.After performing the above iterative calculations for all discrete elements, the dynamic simulation results of the studied object are presented.


## 2.2. Shum Wan Road Landslide

On August 13, 1995, a landslide took place at the hillside above Shum Wan Road, Hong Kong. It caused the collapse of a 30 m long section of Nam Long Shan Road, as shown in Figure1. The literatures reported that the landslide debris crossed Shum Wan Road and damaged three shipyards and a factory near the seafront [46–48]. Prior to the landslide, the hillside was densely vegetated and had an overall gradient of about 27°. The geology at the landslide area comprised a thin mantle of colluvium overlying partially weathered fine-ash to coarse-ash crystal tuff. The landslide resulted in a 70 m high scar, with a width varying from about 50 m just below Nam Long Shan Road to about 90 m above Shum Wan Road. Figure 2 shows the map of the landslide based on topographic survey, geological mapping, and field observations. The upper part of the landslide surface was concave in shape and was up to about 12 m in depth below the prefailure ground surface. The landslide released about 2.6∗104m3 of soil and rock. The debris was deposited on Shum Wan Road and the reclaimed land, spreading over an area of about 0.5 ha. The failure was caused principally by the presence of weak layers in the ground, ingress of water during prolonged heavy rainfall, and partial blockage of the drainage system.Figure 1
The Shum Wan Road landslide on August 13 1995 [46].Figure 2
The map of the Sum Wan Road landslide (modified from Knill 1996).


## 2.3. Back Analysis of the Landslide

The terrain, sliding surface, and disaster scope of the Shum Wan Road landslide are given in the literature [46–48]. A geometric terrain model was established according to the contour of the terrain where the Shum Wan Road landslide is located, and the terrain was discretized by triangular wall elements, as shown in Figure 3.Figure 3
Numerical model of the Shum Wan Road landslide.To obtain a discrete element model for simulating the landslide, it is necessary to use ball element to fill the landslide body. Considering the simulation effect and calculation cost, a ball element with a radius of 0.75 m was selected for the numerical simulation according to the scale of the landslide. The outer surface of the landslide is determined by the ground surface before failure and sliding surface, which are discretized by wall elements, thus forming an enclosed space. Ball elements are filled in the enclosed space formed by wall elements (which is equivalent to the initial landslide) and reach equilibrium under gravity. Finally, a numerical model with discrete elements was established as shown in Figure3.According to the information given in the literature [46–48], the friction coefficient between the ball and the wall can be determined, and the density of the ball can be determined through mass equivalent [43, 49]. However, the stiffness and damping coefficient of discrete elements need to be obtained through numerical back analysis. As a result, the mesocalculation parameters of discrete elements that conform to the landslide movement are shown in Table 1.Table 1
Mesoparameters used in the numerical simulation of Shum Wan Road landslide.
Ball elementWall elementBall diameter (m)Ball density (kg/m3)Normal stiffness (N/m)Shear stiffness (N/m)Friction coefficientNormal stiffness (N/m)Shear stiffness (N/m)Friction coefficient1.5026501e101e100.2141e101e100.214When the landslide starts to move, delete the wall elements that act as a constraint on the upper surface of the landslide body (these wall elements are only used to establish a numerical model of the landslide), and then, the ball elements used to simulate the unstable geomaterial start to move down the slope under gravity. Through numerical simulation, the movement process of the Shum Wan Road landslide was obtained, as shown in Figure4.Figure 4
Flow depth contour map of the Shum Wan Road landslide (m).
(a)
Initial moment of the landslide movement(b)
Numerical result at the 5th second(c)
Numerical result at the 10th second(d)
Numerical result of the final depositionIt can be seen from Figure4 that when the slope enters failure state, the landslide begins to move downward along the slope (Figure 4(a)). The movement is mainly restricted and affected by the terrain. The landslide is mainly accelerated in the early stage of movement, the front edge moves faster, and the rear edge moves relatively slowly, and the largest flow depth area is still at the back of the landslide (Figure 4(b)). The movement process of the landslide developed rapidly. As the landslide entered a gentle section, the sliding force it received was reduced, and the bottom resistance increased, so the front edge of the landslide started to decelerate. The landslide almost reached the actual disaster scope within 10 s (Figure 4(c)), and the final deposition was basically formed within 50 s (Figure 4(d)). During this process, the main mass of the landslide went downstream, and most of the mass deposited at the foot of the hillside, thus causing great harm to the road and the infrastructures downstream.In the numerical results shown in Figure4, the landslide movement is generally within the actual disaster scope, and the edge of the landslide is consistent with the actual disaster scope. However, the disaster scope obtained by numerical simulation in the downstream area full of buildings and facilities is different from the actual situation to a certain extent. The reason for this difference is the existence of these buildings and facilities. The buildings and facilities downstream of Shum Wan Road were destroyed by the landslide, but in fact, they also played a role in blocking the landslide.But it is a pity that due to the lack of information on the distribution of buildings and facilities before the occurrence of the landslide as well as their structure and materials, their blocking effects on the landslide cannot be considered in the calculation (even if the information is mastered, it is still extremely difficult to realize a numerical simulation that carefully considers the coupling effect of landslide impact and structural collapse). It can be seen from the numerical results that although the blocking effect of buildings and facilities is not considered, the numerical simulation has generally given a true representation of the movement process of the Shum Wan Road landslide. The numerical results obtained are correct and reliable, and the numerical calculations are also effective and credible. The huge damage to the downstream infrastructure caused by the Shum Wan Road landslide reminds us that it is necessary to make further specific investigations into the impact of the landslide.


## 2.4. Impact Characteristics of Rigid Retaining Wall

When the position of the unstable slope cannot be accurately determined, or the unstable slope is widely distributed, reinforcing the slope costs too much but cannot ensure an effective protection of objects. In this case, it can also be considered to protect the objects by setting a rigid retaining wall upstream.The Shum Wan Road Landslide is a typical roadside landslide. There are many infrastructures in the lower reaches of Shum Wan Road, and the employees working there are also threatened and affected by the landslide. In the environment of Shum Wan Road Landslide, rigid retaining wall is a representative protective structure. Moreover, according to the characteristics and scale of the landslide, it is appropriate to use a rigid retaining wall as a protective structure. Therefore, the Shum Wan Road Landslide was selected as a case to study the impact characteristics of rigid retaining walls.The completed numerical back analysis has led to a deeper understanding of the landslide movement process. Numerical back analysis also makes it possible to further study the impact of landslides on roads and to in-depth study of the mechanical characteristics of blocking structures along roads. Numerical investigations are important basis for the design of protective structures.To understand the impact of landslides on the road, vertical rigid retaining walls are set along the road to block the landslide. The main function of these rigid retaining walls is to block the moving mass to ensure the safety of downstream infrastructures and personnel. The setting of the rigid retaining wall is shown in Figure5.Figure 5
Rigid retaining wall along the road.The numerical investigation of the landslide process under the condition of setting a rigid retaining wall is carried out again, in which the rigid retaining wall adopts the parameters used for the wall elements in Table1. The impact on the rigid retaining wall was recorded, and the landslide movement process and the time histories of the impact were obtained as shown in Figures 6 and 7.Figure 6
Flow depth contour map of landslide process blocked by rigid retaining wall.
(a)
Initial moment of the landslide movement(b)
Numerical result at the 7th second(c)
Numerical result at the 15th second(d)
Numerical result of the final depositionFigure 7
The time histories of the impact.
(a)
The time history of the force (N)(b)
The time history of the moment (N∗m)From the results shown in Figure6, it can be seen that the landslide moves along the slope surface under the gravity after failure (Figure 6(a)). The landslide approached the road in about 7 s (Figure 6(b)). After that, blocked by the rigid retaining wall, the process of the landslide showed an obvious difference from the process without the rigid retaining wall. The moving mass forms a preliminary deposition under the blocking of the rigid retaining wall (Figure 6(c)), and then, the subsequent mass can be affected by the deposition formed before, expand along the formed deposition in the direction of the road, and gradually form the final deposition (Figure 6(d)). The time histories of the impact acting on the retaining wall shown in Figure 7 correspond to the process of blocking the landslide. The numerical method used in this investigation is the discrete element method, and the impact curve of the rigid retaining wall has certain volatility. However, the changing trends of the force and moment on the rigid retaining wall along the road are clear, and the trends of the two time-history curves are in good agreement, and they have a good correspondence with the dynamic and static loads of the landslide. Due to the arrival of the landslide, the rigid retaining wall began to be impacted at 7.4 s. At 27.1 s, the deposition of geotechnical materials was basically completed, and the force also tended to be stable. During this period, the force on the rigid retaining wall fluctuated due to the high-frequency collision between particles, and this fluctuation also exists objectively in the actual impact process. The time history curves of the impact on the rigid retaining wall give an important hint, that is, for the blocking structure, we should not only pay attention to the static force of the landslide stable deposition on it, but also cannot ignore the impact caused by the landslide. In order to change the motion state of the landslide, the impact peak value that the rigid retaining wall should bear is an important factor that must be considered in the engineering design. Since an obvious run-up phenomenon has not been found in the numerical process of this case, the peaks of the two impact time-history curves of the rigid retaining wall are generated by the dynamic impact of the landslide. The time history curves of the impact on the structure are an important reference data for structural design, and numerical simulation is a convenient and reliable way to obtain these data.After the landslide stopped moving, a stable deposition was formed, which produced a stable pressure on the rigid retaining wall. Compared with the dynamic impact during the landslide movement, the final pressure of the deposition on the rigid retaining wall is long-term and is the primary factor that needs to be paid attention to. In Figure6(d), it can be seen that the deposition with high middle and low sides is finally formed behind the rigid retaining wall, so the distribution of the force on the rigid retaining wall along the road direction is uneven. After the deposition became stable, the discrete elements formed a contact force distribution on the rigid retaining wall, as shown in Figure 8.Figure 8
Contact force distribution on rigid retaining wall in stable state (N).From the contact force distribution shown in Figure8, the outline of the contact distribution is basically similar to the shape of the deposition. The contact force on the retaining wall is discrete, and these discrete contact forces are caused by the contact between discrete elements and the wall. Although the distribution of the position and magnitude of the contact force on the retaining wall looks a little messy under the transmission of the force chain, these discrete forces still form certain regularity as a whole. For example, the overall law of contact force is that the bottom of the middle area is more stressed, and the edge of the deposition is less stressed. In the example of this study, although the contact force on the rigid retaining wall is mainly in the horizontal direction, it is formed by force transmission under gravity, so the magnitude of the contact force is closely related to the local height of the deposition.Since the numerical method divides the landslide into many discrete elements, the simulation is inconsistent with the actual continuous state of the landslide. Although numerical simulations yield discrete contact force distributions on rigid retaining walls, they cannot be directly applied to the design of protective structures. These contact forces need to be counted in an equivalent way, so as to obtain the results required by the design, which are more helpful to engineering applications. The stable force and moment on each section along the road are counted in sections, and the results are shown in Figure9.Figure 9
Stable load distributions on rigid retaining wall along the road.
(a)
Distribution of the force (N)(b)
Distribution of the moment (N∗m)The deposition generally has the characteristics of high in the middle and low on both sides. The distributions of the force and moment on the rigid retaining wall along the road direction are uneven, but their distribution laws are similar. These distributions are mainly controlled by the deposition height, and the distribution characteristics of stable force on the retaining wall along the road are fully demonstrated in Figure9. Compared with the force on each section, the moment is more significantly affected by the height. This is because the soil pressure gradually increases in the vertical direction, so that the moment on each rigid retaining wall is generally in proportion to the square of the deposition height.Further, it is more urgent to obtain a continuous stress distribution that is closer to the actual situation, so the contact force on the rigid retaining wall in the stable state is converted into the stress distribution by means of statistical averaging, as shown in Figure10.Figure 10
The normal stress distribution on rigid retaining wall in stable state (Pa).From the stress distribution on the rigid retaining wall in the stable state shown in Figure10, it can be seen that the law of stress distribution and contact force is similar. The stress is larger at the bottom of the middle area and smaller at the edges of both sides. Compared with the contact force distribution results that cannot be directly applied to engineering, the stress distribution provides more abundant information and can provide a basis for optimized protective structure design. In addition, in order to investigate the nonnegligible dynamic impact, the historical maximum stress on the retaining wall was recorded and counted, and the maximum normal stress distribution in history was obtained, as shown in Figure 11.Figure 11
The maximum normal stress distribution on the rigid retaining wall in history (Pa).Comparing Figures10 and 11, it can be seen that the impact stress on the rigid retaining wall is much larger than that in the stable state. The price of the rigid retaining wall forcing the landslide to stop in a short time is that it has to withstand a larger impact. Interestingly, there is a high stress area in the center of the retaining wall, which is the impact area of geomaterials with high velocity. It is difficult to determine these impact areas and the impact stress on them without numerical methods. Combined with the distribution of the impact area, it can be considered to paste some soft materials for energy dissipation on the upstream surface of the retaining wall, which can reduce the peak of the impact force by prolonging the impact load time. Of course, there may be some differences between the numerical calculation and the actual situation, which allows the scope of mitigation materials to be appropriately enlarged. Under complex terrain conditions, numerical simulations can be used to obtain the impact characteristics of the protective structure that are consistent with the actual situation. The results of these impact characteristics can provide great help for the design of protective structures.


## 3. Discussion

This study shows that the discrete element method can make a good back analysis of the movement process of the landslide, and the impact characteristics can be obtained by this method. Due to the theoretical assumptions and simplifications of the discrete element method, it is impossible for numerical simulations to fully represent all the mechanisms in the dynamic process of landslides. We need to discover them and improve the numerical models to make the simulations closer to the real situation. For example, a high-speed landslide may shovel weak material along the path, making itself bulkier, thereby increasing its mobility. Although some studies have attempted to simulate the entrainment phenomenon along the path, the targeted discrete element model still needs to be further improved [49]. At the same time, due to the various forms of protective structures, it is necessary to combine the characteristics of different structures to study the coupling effect of landslides and structures [12], which can not only better evaluate the working condition of the structure, but also the possible consequences after the failure process.Due to the lack of the available data in the literature, the topographic map of the area where the landslide is located is also limited, resulting in the deposition being restricted by the lateral boundary under the blocking of the retaining walls set along the roadside (for the purpose of avoiding the escape of discrete elements). The boundary constraints on discrete elements can be seen on the left side of Figure6, but the influence of this constraint is limited because the height of the deposition at the boundary is relatively small. In addition, there is a protruding “ridge” on the slope, and some protection works can be built along the “ridge” in order to effectively utilize the terrain and maximize the protection effect, so that the landslide could not affect the greater range.In this study, the numerical simulation method is used to obtain the impact process of the Shum Wan Road landslide on the rigid retaining wall along the road. The main reason for assuming that the protective body is rigid is that in the case where the composition of the geomaterials is mainly soil with small particle size rather than rocks with large particle size, it is not suitable to use a flexible mesh structure for protection. Furthermore, the rigid retaining wall has been widely used and is convenient in construction. This study is the basis for the design of the protective structure. However, further work still needs to be done to obtain a protective structure that can actually be used in engineering practice, including the selection of appropriate building materials (generally reinforced concrete), as well as specific material ratios and reinforcement designs. Besides, stability studies on the foundations of these protective structures are also indispensable. If these works are neglected, the protective structures may overturn without structure failure.Protection engineering is systematic. In addition to passive protection, active protection work should also be done, including some engineering treatments through field investigation and calculation analysis, such as local drainage and reinforcement of the slope. When conditions permit, an open structure to slow down the landslide velocity and conduct a layer-by-layer obstruction can be supplemented on the slope, so as to reduce the impact and risk of damage to the rigid blocking structure at the end. Although the scale of the Shum Wan Road landslide is not very large, the density of important buildings and people downstream is high. In addition to setting up effective protective structures, relevant contingency plans should also be formulated in advance to deal with worse situations beyond design expectations. The research of the rigid retaining wall under the impact of landslide will still be limited in the future application, and the limitation is derived from the numerical method used. Due to the limitations of current computing power, it is impossible for discrete elements to exactly match the actual particle size. At the same time, the determination method of discrete element size and microscopic parameters that match the actual situation is still not mature. These will restrict the application of the method in future prediction research. Nevertheless, with the help of case calibration and verification, the discrete element method is still a reliable and effective method for in-depth study of the complex process of landslide impact on protection structures.


## 4. Conclusion

Rigid retaining wall is a common passive engineering structure used to block landslides and has been widely used because of its economical and reliable properties. However, due to the limited knowledge of rigid retaining walls under the impact of landslides, more conservative solutions have to be adopted for retaining walls under complex conditions, which limits the space for structural optimization.In this study, the discrete element method is used to numerically simulate the Shum Wan Road landslide in Hong Kong, China, and the movement process and development characteristics of the landslide are reproduced. The numerical results obtained are in good agreement with the actual situation, so that the correctness and effectiveness of the calculation are verified.Based on the back analysis of the landslide, rigid retaining walls were set up along Shum Wan Road to investigate the characteristics of the impact on it. A deposition was formed under the blocking of the retaining wall, and the downstream infrastructure and personnel safety can be guaranteed, thereby eliminating the threat. The time history curves of the impact on the rigid retaining wall can clearly show the impact peak caused by the dynamic force of the landslide, which is an important factor for protective structures that cannot be ignored. At the same time, the space force on the retaining wall is also uneven. Although the contact force acting on the retaining wall obtained by the discrete element method is discontinuous, the force distribution of the segmented retaining wall can be obtained by using statistical approach, as well as the distributions of static and peak stresses, which are very important for the design of the retaining wall. The investigation and processing of data in this study enables the obtained information of static and dynamic forces acting on the retaining wall to be used for further optimization of the protective structure. For a complex terrain, the temporal and spatial distributions of dynamic impact from landslides are difficult to predict. Numerical investigations based on discrete element method can obtain this information economically and efficiently and thus can be regarded as a reliable method to study the mechanical characteristics of protective structures.



---

*Source: 7835780-2022-05-19.xml*