Abstract
Majiagou landslide, a major ancient landslide in Three Gorges Reservoir region, is located in the high earthquake area of southwest China. The 2013 Badong earthquake caused an obvious deformation of landslide monitored by the sliding inclinometer. A strong earthquake may induce the reactivation of ancient landslide. So, it is necessary to research the seismic dynamic response of Majiagou landslide. For this purpose, discontinuous deformation analysis (DDA), improved by introducing the artificial joint and viscous boundary, is applied in this study. The displacements at monitoring points caused by Badong earthquake are calculated and compared with the field data, verifying the numerical method and model. Further, a strong earthquake with the peak acceleration of 1 g is assumed to act on the landside, the initiation and evolution process of landslide is simulated, and the movement features of landslide are discussed. The dynamic failure of landslide and the local amplification of seismic wave can be embodied, indicating that the improved DDA provides an alternative approach for analyzing the seismic dynamic response of jointed rock.
1. Introduction
Earthquakes mostly occur at the junction of plates, where the terrain is quite variable and large number of slopes are distributed widely. Further, earthquakes may trigger the sliding failure of slopes, resulting in serious threats to human lives and natural environment. In the past few decades, many such disasters have been reported [1–6]. For example, the 1999 Chi-Chi earthquake () induced 9272 landslides, which caused 2400 deaths and over 8000 injured. More than 20,000 lives were taken by 60,104 landslides triggered by the 2008 Wenchuan earthquake (). Recently, the 2013 Lushan earthquake in Sichuan () caused 3883 landslides, which killed 196 people and injured 11,500 people. Therefore, for disaster reduction and prevention, it is crucial to study the dynamic response of slope and evaluate its stability or postfailure impact under the action of earthquake.
In the field of geotechnical engineering, numerical simulation is a powerful tool for investigating the mechanical behavior of rock mass. The involved methods can be roughly divided into two categories: continuous methods and discontinuous methods. These continuous methods are usually subjected to small deformation assumption and difficult to handle numerous discontinuities in rock mass. Conversely, the discontinuous methods are capable of modeling large deformation behavior of rock mass containing discontinuities. Among them, smoothed particle hydrodynamics (SPH), discrete element method (DEM), and discontinuous deformation analysis (DDA) are widely used for the simulation of earthquake-induced landslides. For the run-out analysis of flow-like mass, SPH simulation is suitable (Dai et al. [7], Mao et al. [8], and Zhang et al. [9]), while for the landslide controlled by joints, DEM or DDA method is applicable. Tang et al. [10], Lu et al. [11], Zhang et al. [12], and Luo et al. [13] applied DEM to simulate the kinematic behaviors of earthquake landslides and analyze the landslide formation mechanism and impact areas. Zhou et al. [14] employed the FDM-DEM-coupled method to study the dynamic process of Donghekou landslide triggered by Wenchuan earthquake. DEM is force-based approach, and exhibits two limitations: artificial damping and very small time step, but DDA as energy-based approach can overcome these limitations. Hatzor et al. [15] firstly verified the ability of DDA to simulate the seismic dynamic behavior of blocky rock and predicted the failure mode of jointed rock slope. Thereafter, Wu et al. [16, 17] proposed a novel algorithm of seismic input for DDA simulation and applied the seismic DDA to model the kinematic behavior of landslides induced by Chi-Chi earthquake. Zhang et al. [18, 19], Song et al. [20], and Peng et al. [21] enriched the DDA method and studied the mobility features of Daguangbao and Donghekou landslides triggered by Wenchuan earthquake.
The aforementioned studies usually assumed the coseismic landslides as block systems and mainly concentrated on their kinematic characteristics, e.g., the run-out distance, impact area, or deposit topography. Notably, few of them researched the initiation mechanism of landslide and the vital role of discontinuities in seismic wave propagation. Aiming at this, the improved DDA method introducing the artificial joint and viscous boundary is employed in this paper to study the reactivation process and dynamic response of Majiagou landslide under earthquake and analyzes its failure mechanism and movement patterns.
2. Geological Setting of Majiagou Landslide
The Majiagou landslide is located on the left bank of Zhaxi River, a tributary of Yangtze River. Figure 1 shows (a) the location, (b) an air photo, and (c) the geological section of Majiagou landslide along the main sliding direction, 291°. The landslide is developed on the front edge of a giant ancient landslide deposit and bounded by seasonal gullies. The sliding mass has a tongue shape, a length of 537 m, an average width of 180 m, and a maximal thickness of 29.5 m. The altitude difference of landslide is about 160 m, decreasing from 290 m to 130 m above sea level. The slope terrain can be divided into three parts. The gradient of the upper part varies from 10° to 26°, the middle part’s is about 12° to 20°, and the lower part’s ranges from 30° to 45°.
(a) Location
(b) Air photo
(c) Geological section along the sliding direction
Explored by the drill holes and shallow wells, the sliding zone is in the contact belt between the surface deposit and the underlying bedrock. The sliding surface has an average gradient of 18°. The upper sliding mass is mainly composed of sand and mudstone debris mixed with gravel. The bedrock of Majiagou area is the interbedded argillaceous siltstone and silty mudstone of Jurassic Suining Formation. Among them, the mudstone can easily form weak interlayer after weathering and softening with water. The bedrock with the occurrence of 270°~290°∠15°~18° (dip direction/dip angle), dips gentler than, and in the same direction as, the slope, which obviously facilitates a landslide phenomenon in terms of geological structure.
In order to grasp the deformation features and evolution trend of Majiagou landslide, two deep displacement monitoring holes are, respectively, set at altitudes 177 m (ZK1) and 190 m (ZK2), and long-term regular observations on landslide deformation are carried out.
3. The Improved DDA Method
In 1988, a novel discrete numerical method, namely, discontinuous deformation analysis (DDA), is developed by Dr. Shi [22] to study the statics and dynamics of rock block system. The rock mass is assumed as a block assembly delimited by discontinuities. The Newton’s second law of motion is employed to describe the motion of each block, and the contact penalty springs are applied to avoid the penetration between blocks. Different from DEM, DDA derives the governing equations by minimizing the total potential energy of block system and implicitly solves them.
3.1. DDA Governing Equations
The governing equations of DDA have the following form: where
and denote the generalized displacement and force vectors; , , and are the mass, damping and stiffness matrices, respectively; in the current DDA program, no damping is considered, while the dynamic coefficient (valued from 0 to 1) for transferring velocity between time steps is used to dissipate energy [23–25]; is the interval of time step; is the initial velocity of the current time step, that is, the terminal velocity of the previous time step.
Replacing and of Equation (1) and rearranging it, the DDA governing equations can be written as: where is a 6th order submatrix, depends on the material properties of block , and is related to the contact state between block and ; and are 6th order vectors, is the generalized displacement vector of block , and is the generalized force vector acting on block .
Equation (3) is repeatedly solved to achieve no-penetration state of block system by adjusting the contact springs between adjacent blocks.
3.2. Fracturing Modeling
The original DDA method deals with the blocky rock. However, many rock masses are not completely discrete, but intermittent or even intact. In order to expand the research object of DDA, the authors [26] incorporated the artificial joint companied with its fracturing algorithm into DDA to model the continuous-discontinuous process of rock mass. The rock mass is cut into triangular block system by real and artificial joints, and the artificial joints bond adjacent blocks by the cohesive force. The cohesive force constantly evolve as the rock mass moves and deforms, and in the normal direction, it can be depicted by Figure 2.
When the normal cohesive force at an artificial joint decreases to zero, the artificial joint will occur tensile failure and becomes a real joint. The shear failure of artificial joint can be judged by Mohr-Coulomb criterion: in which is the tangential contact force; is the length of artificial joint; and are cohesion and inner friction angle of rock, respectively.
3.3. Viscous Boundary
Generally, the numerical calculation model is a finite domain taken from an engineering site, and the artificial boundary may cause the reflection of stress wave, which is inconsistent with the reality. The viscous boundary proposed by Lysmer et al. [27] can remedy this defect by absorbing the wave energy at the model’s boundary. The authors introduced it into DDA for solving the dynamic problems [28].
As displayed in Figure 3, two pairs of normal and shear dashpots are applied at each boundary block. The viscous tractions produced by the dashpots are opposite and proportional to the velocity of block, and they can be denoted as follows: where is the density of rock; and are the normal and shear velocity components of block, respectively; and are the p-wave and s-wave velocities, respectively.
Based on the principle of minimum potential energy, the contributions of viscous tractions to the governing equations can be obtained. The derivation process is detailed in literature [28].
4. Dynamic Response of Majiagou Landslide
4.1. Numerical Model and Parameters
Based on the field survey, the calculation model of Majiagou landslide is established (as shown in Figure 4). In the model, the bedrock is regarded as a complete rock block, several weak bedding planes of mudstone are treated as joints, and the sliding body is divided into triangular blocks bonded by the artificial joints. The mechanical parameters satisfying the Weibull’s distribution are assigned to the sliding blocks to consider the material heterogeneity. Some monitoring points are set up on the surface and inside of landslide. It should be noted that the negative displacements, velocities, and accelerations in the -direction obtained by the simulation represent the downhill movements.
As for the seismic input, the original DDA directly applies the time-dependent acceleration on each block as body force, which may induce incorrect absolute movements of blocks. Some efforts have been made on the seismic input method for DDA [29–31]. Commonly, an alternative method is applying the seismic acceleration only to the base block with a large mass. In this simulation, we also adopt this method, that is, a seismic wave is applied to the large base block. The viscous boundaries are implemented at the model’s lateral boundaries.
The material parameters listed in Table 1 are used in the calculation. The shear strength parameters of sliding belt in natural state are obtained by consolidation quick shear tests, and those in saturated state are from direct shear tests after several wet-dry cycles [32]. The sliding body is a mixture of soil and gravel, and its shear strength parameters are obtained through large-scale direct shear tests [32]. Due to the lack of test data, the tensile strength parameters are estimated from the cohesions.
In addition, DDA dynamic analysis also involves some very important calculation parameters, such as the time interval, maximum displacement ratio, and dynamic coefficient. Some researchers validated DDA by analytical solutions and shaking table experiments and found that when the time interval is small, the maximum displacement ratio is close to the time interval, and the dynamic coefficient is greater than 0.95, more reasonable results can be obtained [33–35]. Based on this, the calculations in this paper set the time interval to 0.001 s, the maximum displacement ratio to 0.001, and the dynamic coefficient to 0.98.
4.2. The Badong Earthquake
On December 16, 2013, an earthquake with Ms 5.1 occurred in Badong county, Hubei province, China. An obvious deformation was monitored at the toe of Majiagou landslide. Figure 5 plots the displacement curves measured by the sliding inclinometers at ZK1 and ZK2 from December 14, 2013, to April 22, 2014, and the displacement variation before and after the earthquake. It can be seen that the depths of sliding surface at ZK1 and ZK2 were 30 m and 26 m, respectively. After the earthquake, the landslide produced a downhill deformation, which has the same trend as before the earthquake, and gradually recovered until April 22, 2014. The earthquake-induced displacement at ZK1 decreased linearly with depth, and the maximum of 2.4 mm occurred on the surface. However, since ZK2 was in the water level fluctuation zone, the displacement was more complicated, and the maximum of 4.6 mm occurs 8 m below the surface.
(a) The displacement at ZK1
(b) The displacement at ZK2
(c) The displacement variation before and after the earthquake
According to the magnitude and epicenter distance of Badong earthquake, the peak seismic acceleration in the Majiagou area can be calculated as 0.24 m/s2 [36]. Adjust the amplitude of the recorded seismic wave and convert it to the direction of sliding surface to obtain the inputted seismic wave (as shown in Figure 6). In December, the water level of Three Gorges Reservoir is 175 m, and correspondingly, the sliding body below 175 m is set as saturated. Figure 7 displays the earthquake-caused displacement at ZK1 and ZK2 simulated by DDA. We can find that the simulated displacement distribution along the depth is in good agreement with the monitoring results. The horizontal displacement distribution in the landslide, given in Figure 8, indicates that Badong earthquake caused Majiagou landslide to slide downhill as a whole, and the middle surface was deformed greatly. However, it needs to be confirmed by more monitoring data.
4.3. A Strong Earthquake
The Majiagou landslide lies in an earthquake-prone zone. A strong earthquake may trigger landslide activity and cause serious disaster. Therefore, predicting the reactivation process and understanding the failure features and mechanism of landslide are critical. Besides, the effect of the fluctuating reservoir water level cannot be ignored either. Herein, a strong earthquake is assumed to act on Majiagou landslide, its failure processes under 145 m and 175 m water level are simulated, and the dynamic responses at the different parts of landslide are investigated.
In this simulation, considering that the horizontal seismic wave has a greater impact, a synthetic seismic wave with the peak acceleration of 1 g, plotted in Figure 9, acts horizontally on Majiagou landslide. As shown in Figures 10 and 11, the landslide is reactivated under the strong shaking. Firstly, the landslide slightly slides downward under the initial seismic force, and several tensile cracks are initiated at the top part. Then, as the shear strength is degraded by the reservoir water, the toe of landslide is partially fragmented and disconnected from the middle part under the intense shaking. After that, some intersecting cracks are generated at the middle of landslide due to the loss of resistance from the toe. Finally, fragments are continuously produced at the toe front and fall into the Zhaxi river, and the gully between the toe and the middle of landslide is further widen. After the earthquake, the landslide is roughly disintegrated into four sections. Under the condition of 175 m water level, the landslide failure is more serious, and even block the Zhaxi river. This is closely related to the large range of sliding zone deteriorated by the reservoir water.
Figure 12 displays the horizontal velocity histories at the monitoring points of landslide under the earthquake and 175 m water lever. It can be observed that the sliding mass movements can be classified into four patterns: (1)Monitoring point S1 is at the front of landslide toe. After it is fragmented, it is free to fall into the river and climb to the opposite bank. It moves very violently and has the maximum velocity of 33.4 m/s at 16 s(2)Monitoring points S2~S4 are located at the toe and middle of landslide. At about 10 s, the toe is broken, resultantly possesses high momentum and starts to slide as a unit at the velocity of >10 m/s. They have four significant velocity cycles. In the first cycle, S3 reaches the peak velocity of 14.1 m/s at 17 s. In the second cycle, S2 reaches the peak velocity of 13.5 m/s at 26.5 s. The movement of S4 at the landslide middle is in the same pattern, but much slower(3)Monitoring points S5~S8 are located at the middle and top of landslide. Between 0 and 22 s, they have great movement with relatively low velocities of <4 m/s. After 22 s, they have the velocities around 0, thus basically stay static(4)Monitoring points D1~D5 are inside the landslide. They move almost synchronously with the same waveform. Among them, D5 at the landslide top has a higher peak velocity of 2.6 m/s
(a) Monitoring point S1
(b) Monitoring points S2~S4
(c) Monitoring points S5~S8
(d) Monitoring points D1~D5
Figure 13 compares the peak velocities at the monitoring points under 145 m and 175 m water level. It demonstrates that the sliding resistance at the landslide toe decreases as the water level rises, resulting in a faster slide. On the landslide surface, the velocities at the toe are faster, while the velocities at the middle and top have no much difference. Inside the landslide, the velocities increase roughly linearly with the altitudes. Figure 14 displays the peak accelerations at the monitoring points under 145 m and 175 m water level. We can observe the local amplification of seismic wave, indicating that the proposed method enables the simulation of wave propagation in rock mass. The maximum acceleration amplification factor is at the toe front, reaching 20 and 7.2 at the 175 m and 145 m water level, respectively. While at the top shallow, the amplification factor is minimum, 1.7 and 1.9 at the water level of 175 m and 145 m, respectively. The amplification factors inside the landslide increase with the altitudes, reflecting the amplification effect of bedding planes on the seismic wave.
5. Conclusions
In this paper, the improved DDA method introducing the artificial joint and viscous boundary, is employed to investigate the dynamic response of Majiagou landslide. The deformation of landslide induced by the 2013 Badong earthquake is analyzed, and the numerical results agree well with the monitoring data, indicating that the numerical method is applicable. Based on this, considering a hypothetical powerful earthquake and the fluctuation of reservoir water level, the failure and movement process of Majiagou landslide is simulated. Some numerical observations improve the knowledge of the reactivation mechanism and kinetic characteristics of landslide. (1)The failure initiated from the tension cracks at the landslide top. Then, the landslide toe deteriorated by the reservoir water is broken under the continuous strong shaking. Subsequently, the landslide middle loses the resistance of toe and fractures. Finally, the landslide is broken up into several sections, and some fragments fall from the toe and block the Zhaxi river. It can be concluded that the reactivation of Majiagou landslide is closely related to the strength degradation of landslide toe by the reservoir water. Therefore, the water level is a vital factor for the stability of reservoir bank landslides(2)The sliding mass has various movement patterns. The landslide toe, especially its front, moves at high velocities, which have several obvious acceleration-deceleration cycles. The middle of landslide moves greatly in the early and midterm of earthquake and then gradually attains the stable state on its own sliding resistance. The movement of bedrock is basically synchronized, while the bedding planes can speed up it(3)The seismic wave is significantly amplified on the landslide surface. The larger amplification factors are concentrated at the landslide toe. At the middle and top of landslide, the amplification factors are not sensitive to the altitude, but at the shallow, the amplification factors are smaller. The amplification effect of bedding planes inside the landslide on the seismic wave is noticeable
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge and greatly appreciate the support for this study by the National Key Research and Development Plan (2017YFC1501304) and National Natural Science Foundation of China (41772328 and 11672360).