# The Impact of Microearthquakes Induced by Reservoir Water Level Rise on Stability of Rock Slope

Authors: Dongliang Li; Xinrong Liu; Xingwang Li; Yongquan Liu
Journal: Shock and Vibration
 (2016)

---

## Abstract

In order to study the impact of frequent microearthquakes induced by water level rise on the stability of rock bedded slopes in the Three Gorges Reservoir (TGR) area, Zhaoshuling Landslide (a representative slope) is selected to study. Safety factors based on probability statistics andF
L
A
C
3
D are used for numerical simulation (under the operating condition that five earthquakes of Intensity IV are applied to slope in succession after water level rises from 145 m to 175 m). Then the slope’s dynamic stability characteristics and failure mechanism are analyzed. The study shows that slope deformation is evidently the result of thrust load. The deformation starts from the steeply dipping segment in the middle part of slip mass and is controlled by the soft interlayer. Shear failure tends to occur along the soft interlayer and the horizontal slip displacement increases from the rear to the front of the slope. The steeply dipping segment shows a general downslide trend. Although the gentle slope platform on the rear edge is relatively stable, it is vulnerable to tensile fractures which are precursors of landslide. Under the same failure probability, as the number of microearthquake occurrences increases, the safety factor of slope under microearthquake action decreases gradually.


---

## Body

## 1. Introduction

The Three Gorges Project (TGP) is a hydraulic engineering project on the Yangtze River that has attracted worldwide attention. It involves the construction of a long and narrow channel reservoir over 600 km long and 1-2 km wide [1–3]. However, this huge artificial water project leads to a number of pressing issues in addition to great benefits to the economic development of China. According to relevant data [4–6], the length of the unstable segment in the TGR area is 385 km, and 1190 slumped masses (with a total volume of around 3.4 billion cubic meters) have been identified along the reservoir bank of the TGR. In the event of river blockage and water level rise due to failure and instability of the reservoir bank slope, catastrophic consequences would be caused. For this reason, scholars have made relevant research on the stability of reservoir bank slopes in the TGR area. Zuo et al. [7] analyzed the relationships between the variations in the reservoir water level and width of the bank slope’s structural plane and the slope stability with the Strength Reduction Method. Luo et al. [8] found that the horizontal and vertical displacements of a slope abruptly increased in the Three Gorges area after reservoir impoundment and rainfall. Based on the comparison of the position, mineral and structure between incompetent beds, and sliding zones of two giant deep-seated translational rock landslides, Chai et al. [9] discussed the correlations between incompetent beds and giant landslide. With Anlesi Landslide as a typical case, Jian et al. [10] studied the mechanism of gentle dip translational rock slides in the TGR area. Based on spatial data mining and knowledge discovery, Wang et al. [11] carried out landslide mechanism analysis in the Three Gorges.In addition to the direct influence on the stability of reservoir bank slopes, the variation in the reservoir water level can also cause changes of the stress field and seepage field in the rock mass, thus inducing earthquakes [12–15]. According to relevant literature and monitoring data [16–22], the frequency of earthquake activities in the reservoir area soars after impoundment, surging up to 15 occurrences per day. The earthquake motion with such a high frequency and long term is inevitably a great challenge for the stability of the slope in the TGR area. For the stability of rock slopes under the earthquake action, Tang et al. [23] analyzed the relations between the rock landslide distribution and five main factors (fault proximity, epicenter proximity, peak ground acceleration, slope angle, and lithology). Gischig et al. [24] presented a series of 2D distinct-element numerical models which were aimed at clarifying interactions between earthquakes and large rock slope instabilities. Guo et al. [25] collected 54 landslides with travel distances of 347–4170 m triggered by Wenchuan Earthquake in 2008 to discuss the effectiveness of various influential factors on landslide travel distance, and his research results revealed that rock type, sliding source volume, and slope transition angle were the predominant factors on landslide travel distance. Besides, Xu [26] studied the stability of steep rock slopes under microearthquakes and Jiang [27] analyzed the deformation mechanism and dynamic response characteristics of the first typical loess landslide segment in the TGR area under microearthquakes. However, little research has been conducted on the impact of reservoir water level fluctuation and frequent microearthquake occurrences on the dynamic stability of the bedded rock slope.To sum up, with Zhaoshuling Landslide (Badong County, China) in the TGR area as a representative case, numerical simulation is conducted for the rock bedded slope with a soft interlayer under multiple microearthquakes induced by reservoir water level rise and the dynamic stability and failure mechanism of such slopes are analyzed based on the slope displacement deformation and dynamic response to provide theoretical bases for landslide control, slope stabilization, and disaster prevention.


## 2. Calculation of Dynamic Safety Factor

### 2.1. Time-History Curve of Dynamic Safety Factor

The safety factor of the slope sliding failure mode under the earthquake action can be determined with the limit equilibrium theory by calculating the earthquake inertia force on the rock mass at random timet.Assume thatF
→
t is the resultant force of gravity G
→
t and earthquake inertia force W
→
t on the rock mass at time t. The resultant force is projected onto the potential sliding direction to obtain the sliding force S
t as follows: (1)
S
t
=
F
→
t
·
i
→
,where i
→ is the unit vector of the potential sliding direction. F
→
t is projected onto the normal direction of the sliding surface to obtain the normal reacting force.According to the Coulomb Friction Laws, the slide-resistant forceR
t can be determined as follows: (2)
R
t
=
F
→
t
·
n
→
·
tan
⁡
φ
t
+
A
·
c
t
,where n
→ is the normal unit vector of the sliding surface, A is the contact area between the sliding mass and the rock mass below it, and φ
t and c
t are the internal friction angle and cohesion of the sliding surface at time t, respectively.Therefore, the dynamic safety factor can be expressed as follows:(3)
k
t
=
F
→
t
·
n
→
·
tan
⁡
φ
t
+
A
·
c
t
F
→
t
·
i
→
.Based on the formula above, a computation program can be developed using Fish language inF
l
a
c
3
D to obtain the time-history curve of dynamic safety factor for the slope under seismic action, thus providing a basis for stability analysis of slopes under frequent microearthquakes.


### 2.2. Safety Factor Based on Probability Statistics

Earthquake is normally assumed to be a random variable of time. Accordingly, the dynamic safety factor of slope under seismic action can also be seen as a random function of time, which enables the dynamic stability evaluation index of slope to be determined based on the reliability theory [28].Assume that, under seismic action, when the minimum value of the time-history curve for the dynamic safety factor is used as the evaluation index, the highest level of safety is expected (i.e., the probability of failure is 0); when the maximum value of the time-history curve is used as the evaluation index, the slope stability may be overestimated and the highest level of risk is expected (i.e., the probability of failure is 1). Based on the assumption above as well as the fundamentals of probability and statistics, if the safety factor is taken as the critical index, the evaluation index for dynamic stability of the slope can be transformed into the probability of the safety factor belowF
R due to the earthquake.As shown in Figure1, the area of the shaded part is the probability of failure P
f. Based on the principle of probability analysis, the reliability β can be determined as follows: (4)
β
=
μ
F
-
F
R
σ
F
,where μ
F is the mean of safety factors and σ
F is the standard deviation of safety factors. According to the formula above, (5)
F
R
=
μ
F
-
β
σ
F
.Figure 1
Distribution of dynamic safety factors.The time-history curve for the dynamic safety factor of the slope is discretized and the following formulas can be obtained with the statistical method:(6)
μ
F
=
1
N
∑
i
=
1
N
F
t
i
,
σ
F
=
1
N
-
1
∑
i
=
1
N
F
t
i
-
μ
F
2
,where N is the total number of discrete segments of the time axis and F
t
i is the average safety factor in the ith segment. When F
t
i is in normal distribution, a one-to-one correspondence can be established between β and probability of failure P
f. See Table 1.Table 1
Correspondences betweenP
f and β.
P
f
β
0.5
0
0.25
0.67
0.1
1.28
0.05
1.65
0.01
2.33
0.001
3.10
0.0001
3.72
0.00001
4.25


## 2.1. Time-History Curve of Dynamic Safety Factor

The safety factor of the slope sliding failure mode under the earthquake action can be determined with the limit equilibrium theory by calculating the earthquake inertia force on the rock mass at random timet.Assume thatF
→
t is the resultant force of gravity G
→
t and earthquake inertia force W
→
t on the rock mass at time t. The resultant force is projected onto the potential sliding direction to obtain the sliding force S
t as follows: (1)
S
t
=
F
→
t
·
i
→
,where i
→ is the unit vector of the potential sliding direction. F
→
t is projected onto the normal direction of the sliding surface to obtain the normal reacting force.According to the Coulomb Friction Laws, the slide-resistant forceR
t can be determined as follows: (2)
R
t
=
F
→
t
·
n
→
·
tan
⁡
φ
t
+
A
·
c
t
,where n
→ is the normal unit vector of the sliding surface, A is the contact area between the sliding mass and the rock mass below it, and φ
t and c
t are the internal friction angle and cohesion of the sliding surface at time t, respectively.Therefore, the dynamic safety factor can be expressed as follows:(3)
k
t
=
F
→
t
·
n
→
·
tan
⁡
φ
t
+
A
·
c
t
F
→
t
·
i
→
.Based on the formula above, a computation program can be developed using Fish language inF
l
a
c
3
D to obtain the time-history curve of dynamic safety factor for the slope under seismic action, thus providing a basis for stability analysis of slopes under frequent microearthquakes.


## 2.2. Safety Factor Based on Probability Statistics

Earthquake is normally assumed to be a random variable of time. Accordingly, the dynamic safety factor of slope under seismic action can also be seen as a random function of time, which enables the dynamic stability evaluation index of slope to be determined based on the reliability theory [28].Assume that, under seismic action, when the minimum value of the time-history curve for the dynamic safety factor is used as the evaluation index, the highest level of safety is expected (i.e., the probability of failure is 0); when the maximum value of the time-history curve is used as the evaluation index, the slope stability may be overestimated and the highest level of risk is expected (i.e., the probability of failure is 1). Based on the assumption above as well as the fundamentals of probability and statistics, if the safety factor is taken as the critical index, the evaluation index for dynamic stability of the slope can be transformed into the probability of the safety factor belowF
R due to the earthquake.As shown in Figure1, the area of the shaded part is the probability of failure P
f. Based on the principle of probability analysis, the reliability β can be determined as follows: (4)
β
=
μ
F
-
F
R
σ
F
,where μ
F is the mean of safety factors and σ
F is the standard deviation of safety factors. According to the formula above, (5)
F
R
=
μ
F
-
β
σ
F
.Figure 1
Distribution of dynamic safety factors.The time-history curve for the dynamic safety factor of the slope is discretized and the following formulas can be obtained with the statistical method:(6)
μ
F
=
1
N
∑
i
=
1
N
F
t
i
,
σ
F
=
1
N
-
1
∑
i
=
1
N
F
t
i
-
μ
F
2
,where N is the total number of discrete segments of the time axis and F
t
i is the average safety factor in the ith segment. When F
t
i is in normal distribution, a one-to-one correspondence can be established between β and probability of failure P
f. See Table 1.Table 1
Correspondences betweenP
f and β.
P
f
β
0.5
0
0.25
0.67
0.1
1.28
0.05
1.65
0.01
2.33
0.001
3.10
0.0001
3.72
0.00001
4.25


## 3. Project Overview

Zhaoshuling Landslide is an area with gentle terrain located in eastern New Urban District of Badong County, China (Figure2). Landslides there with an elevation below 400 m have two-stage platforms at their foreparts. This gentle slope zone is of great importance for Badong’s urban area with a shortage of planned construction land (see Figure 3 for Zhaoshuling Landslide and New Badong County). Therefore, the stability of Zhaoshuling Landslide is directly related to the safety of Riverside Road in the new county. Zhaoshuling is located in a transitional zone between Wuxia Gorge and Ailing Gorge in the middle section of the Three Gorges and is a relatively wide section in the valley area. The topographic profile shows gentle slope platforms alternating with steep slopes. In the aspect of terrain, the rear part and the riverside at the front edge are steep, and the middle part includes three platforms of gentle slopes with an average gradient of 10
°
~
15
°. Between the gentle slopes are steep slopes with a gradient exceeding 35° which is even larger in some parts. Also, gully walls on both sides of the gully in the area are very steep. Some parts of the walls have a gradient exceeding 50° and are potential unstable parts where slope deformation occurs. The axis of rock folds in Zhaoshuling and its surrounding areas is nearly oriented east and west, which is the same as the orientation of steep slopes. Steep slopes are mainly dip slopes. For lithology, the bedrock in the area is mainly argillaceous limestone and marl of Triassic Badong Formation (T
2
b) and silty mudstone and siltstone. It is a typical sliding-prone stratum in the Three Gorges area [29]. Figure 4 is the simplified topographic map of the Zhaoshuling Landslide area.Figure 2
Location of Zhaoshuling Landslide.Figure 3
Zhaoshuling Landslide and the New Badong County.Figure 4
Simplified topographic map of the Zhaoshuling Landslide area.


### 3.1. Characteristics of Landslide Form

Ancient Zhaoshuling Landslide is a bedded rock landslide with its front edge being close to the Yangtze River. Areas with the elevation of 150~200 m belong to the gentle slope platform of stage I; areas with the elevation of 350~425 m belong to the landslide platform of stage II; areas with the elevation of 475~500 m belong to the gentle slope platform at stage III. The remaining areas are slopes with a gradient of 25~50°. The landslide is terraced in the east-west orientation like its surrounding terrains.


### 3.2. Characteristics of Material Structure

Zhaoshuling Landslide is an extra large landslide inT
2
b
3 and T
2
b
2 strata and materials mainly come from the T
2
b
3 stratum. The slumped mass layer in the landslide surface mainly consists of block fractures, cataclastic rock, broken stone containing mud, and broken stone. The bedrock mainly consists of reddish-purple silty mudstone of T
2
b
2 and argillaceous siltstone. The rock mass is in a layer structure with a basically normal sequence of strata. Multiple layers of weak fractured zones can be found in the exposed T
2
b
3 stratum in the landslide area. The lowest sliding zone is around the T
2
b
2
/
T
2
b
3 interface. Characterized by a thin rear edge and thick front edge, this sliding zone is basically the same as the topographic relief in terms of shape.


### 3.3. Characteristics of Geological Structure

The Zhaoshuling area has a special geological structure that the lithologic interface ofT
2
b
3
/
T
2
b
2 is higher than that of other surrounding areas. This is the result of local tectonic uplift caused by two conjugated tensional-shear faults which are nearly oriented south and north. In addition to the characteristics of undulation in the near south-north orientation, the interface is also characterized by undulation in the east-south orientation which is related to tectonism in the near south-north orientation.


### 3.4. Characteristics of Soil in Sliding Zones

The form of sliding zones in Zhaoshuling Landslide is undulating under the control of stratum lithology and attitudes. The front sliding zone lifts up with an inclination of 8~11° and a length of 140~160 m. Materials in the sliding zone are mainly grayish-yellow, grayish-green, and red breccia soil. Some parts of the sliding zone are mylonitized and relatively compact. Its natural volume-weight is 22.9~24.6 kN/m3, larger than that of common gravelly soil. There are many subsliding zones and bedded crushed zones in the soft stratum. Figure 5 is the simplified geological section map of Zhaoshuling Landslide.Figure 5
Simplified geological section map of Zhaoshuling Landslide (plotted based on data of Comprehensive Investigation Institute, China, Changjiang Water Resources Commission (1995)). (1) Limestone and limestone containing mud; (2) mudstone, siltstone, and argillaceous siltstone; (3) soft stratum; (4) landslide deposits; (5) red mud broken stone with soil; (6) gravelly soil.


## 3.1. Characteristics of Landslide Form

Ancient Zhaoshuling Landslide is a bedded rock landslide with its front edge being close to the Yangtze River. Areas with the elevation of 150~200 m belong to the gentle slope platform of stage I; areas with the elevation of 350~425 m belong to the landslide platform of stage II; areas with the elevation of 475~500 m belong to the gentle slope platform at stage III. The remaining areas are slopes with a gradient of 25~50°. The landslide is terraced in the east-west orientation like its surrounding terrains.


## 3.2. Characteristics of Material Structure

Zhaoshuling Landslide is an extra large landslide inT
2
b
3 and T
2
b
2 strata and materials mainly come from the T
2
b
3 stratum. The slumped mass layer in the landslide surface mainly consists of block fractures, cataclastic rock, broken stone containing mud, and broken stone. The bedrock mainly consists of reddish-purple silty mudstone of T
2
b
2 and argillaceous siltstone. The rock mass is in a layer structure with a basically normal sequence of strata. Multiple layers of weak fractured zones can be found in the exposed T
2
b
3 stratum in the landslide area. The lowest sliding zone is around the T
2
b
2
/
T
2
b
3 interface. Characterized by a thin rear edge and thick front edge, this sliding zone is basically the same as the topographic relief in terms of shape.


## 3.3. Characteristics of Geological Structure

The Zhaoshuling area has a special geological structure that the lithologic interface ofT
2
b
3
/
T
2
b
2 is higher than that of other surrounding areas. This is the result of local tectonic uplift caused by two conjugated tensional-shear faults which are nearly oriented south and north. In addition to the characteristics of undulation in the near south-north orientation, the interface is also characterized by undulation in the east-south orientation which is related to tectonism in the near south-north orientation.


## 3.4. Characteristics of Soil in Sliding Zones

The form of sliding zones in Zhaoshuling Landslide is undulating under the control of stratum lithology and attitudes. The front sliding zone lifts up with an inclination of 8~11° and a length of 140~160 m. Materials in the sliding zone are mainly grayish-yellow, grayish-green, and red breccia soil. Some parts of the sliding zone are mylonitized and relatively compact. Its natural volume-weight is 22.9~24.6 kN/m3, larger than that of common gravelly soil. There are many subsliding zones and bedded crushed zones in the soft stratum. Figure 5 is the simplified geological section map of Zhaoshuling Landslide.Figure 5
Simplified geological section map of Zhaoshuling Landslide (plotted based on data of Comprehensive Investigation Institute, China, Changjiang Water Resources Commission (1995)). (1) Limestone and limestone containing mud; (2) mudstone, siltstone, and argillaceous siltstone; (3) soft stratum; (4) landslide deposits; (5) red mud broken stone with soil; (6) gravelly soil.


## 4. Numerical Calculation

### 4.1. Establishment of Numerical Model

A numerical model is established based on the typical geological section in the main sliding direction of the landslide. The model is 1200 m in length and its rear edge has an elevation of 475 m. The controlling structural plane of the slope is simulated using a thin-layer element, while other joint fissures are taken into account by reducing the strength and deformation modulus of the rock mass. In order to ensure the accuracy of dynamic analysis, the size defined by the mesh should be controlled. The mesh size is limited by the shortest wavelength of the input seismic wave. To ensure the authenticity of wave propagation in the media, the maximum sizeΔ
l of the element must be less than 1/10~1/8 of the minimum wavelength. According to the elastic wave propagation theory, the propagation velocities of the longitudinal wave and transverse wave are calculated with the following formula: (7)
C
p
=
E
1
-
μ
ρ
1
+
μ
1
-
2
μ
C
s
=
E
2
ρ
1
+
μ
,where E is the elastic modulus, G is the shear modulus, ρ is the density, and μ is Poisson’s ratio of the material.From formula (7), the propagation velocities of the longitudinal wave and transverse wave in different rock strata are obtained. With the computation speed and the proposed frequency distribution of seismic waves taken into consideration, f
=
10 Hz is used to control mesh accuracy and the maximum mesh size can be determined by the following formula: (8)
Δ
l
=
1
10
·
C
s
f
.The final established numerical model is shown in Figure6. The distribution of rock is shown in Figure 7.Figure 6
Finite element model.Figure 7
Distribution of rock materials.The ideal elastic-plastic model is used. Assuming that the rock materials follow the Mohr-Coulomb strength criterion, the physical and mechanical parameters selected are shown in Table2.Table 2
Physical and mechanical parameters of numerical model.
Rock
ν
K (m·d−1)
E
d (GPa)
c (kPa)
c
s
a
t (kPa)
φ (°)
φ
s
a
t (°)
ρ (kg·m−3)
ρ
s
a
t (kg·m−3)
Sliding body
0.3
10
13
190
130
24
22
2000
2200
Weak layer
0.4
3
0.8
20
12
20
18
2100
2300
Silty mudstone
0.35
0.5
10
200
120
23
21
2300
2590
Argillaceous limestone
0.2
5
16
260
200
32
30
2550
2690
Limestone
0.18
—
32
400
350
35
32
2580
2680
ν is Poisson’s ratio, K is permeability coefficient, E
d is dynamic elastic modulus, c is cohesion, c
s
a
t is saturated cohesion, φ is internal friction angle, φ
s
a
t is saturated internal friction angle, ρ is natural density, and ρ
s
a
t is saturation density.


### 4.2. Seismic Wave Input

Reservoir-induced earthquakes are usually of low magnitude with concentrated epicenters and short duration (usually 5–10 s). Therefore, the time-history curve of a previously recorded real earthquake acceleration is used in this paper and the record of foreshocks is cut. The waveband in 7 seconds during 26 s~33 s is selected as the seismic waveform in this paper (see Figure8). Then, the seismic wave is reduced according to the peak acceleration of the corresponding magnitude. The acceleration time-history curve thus obtained is used for numerical calculation in the study.Figure 8
Original seismic wave time-history curve.
(a)
Acceleration time-history curve
(b)
Velocity time-history curve
(c)
Displacement time-history curveThe seismic wave shown in Figure8 cannot be used directly for numerical calculation. As can be seen from Figure 8(c), the final displacement of the seismic wave is nonzero after an earthquake cycle. It is not the result of residual deformation due to failure of objects, but the result of the incomplete waveform selected artificially or recording errors. Therefore, baseline calibration should be performed on the acceleration time-history curve to eliminate the influence of acceleration error on velocity and displacement, so that the final values of velocity and displacement time-history curves are zero. See Figure 9 for the time-history curve of seismic wave after filtering, baseline calibration, and reduction to Intensity IV (peak acceleration is 0.22 m/s2).Figure 9
Original seismic wave time-history curve (Intensity IV).
(a)
Acceleration time-history curve (Intensity IV)
(b)
Velocity time-history curve (Intensity IV)
(c)
Displacement time-history curve (Intensity IV)AsF
l
a
c
3
D software is used for calculation, two dynamic boundary conditions are provided, that is, the free-field boundary and viscous boundary. In this study, the viscous boundary and local damping (damping coefficient is 0.15) [30] are used for analysis. For this boundary condition, tangential and normal dampers are applied on the boundary to realize the absorption of incident wave energy. In other words, dampers produce tangential and normal force to offset the stress caused by reflected waves. Accordingly, the seismic wave should be input in the form of stress time-history. The expression for conversion of velocity time-history to stress time-history is as follows: (9)
σ
n
=
-
2
ρ
C
p
v
n
,
σ
s
=
-
2
ρ
C
s
v
s
,where σ
n is normal stress, σ
s is shear stress, ρ is medium density, C
p and C
s are the propagation velocity of P-wave and S-wave in the medium, respectively, and v
n and v
s are the velocity in vertical and horizontal directions, respectively.Reservoir-induced earthquakes of the same magnitude may occur multiple times successively within a short period of time. Given this characteristic, five earthquakes of Intensity IV (acceleration amplitude is 0.22 m/s2) are applied successively to simulate the frequent occurrence of microearthquakes.


### 4.3. Numerical Calculation

The following operating condition is mainly considered during numerical simulation: the reservoir water level rises to 175 m from 145 m. The calculation mainly involves dead weight, seepage, and dynamic response analysis.(1) Dead Weight Calculation. During dead weight calculation, gravity acceleration is set first, followed by full constraint on the bottom surface of the model and constraint on the normal displacement around the model. The stress field under gravity is obtained first. The nephogram of the vertical and horizontal stress fields of the slope under natural state is shown in Figure 10.Figure 10
Initial stress field nephogram of the slope under natural state.
(a)
Vertical and horizontal stress field
(b)
Horizontal stress field(2) Seepage Calculation. During seepage calculation, the deep rock mass is thick-layer and compact argillaceous limestone, and as its permeability is weaker than that of the upper rock mass, it can be regarded as a relative water-resisting layer of the slope. The front part of the slope is in direct contact with the reservoir water and thus its pore water pressure is set according to the actual water level. The rear part of the slope is less affected by water level fluctuation and thus the pore water pressure of the nodes can be set according to the groundwater depth revealed during drilling exploration. With bilateral seepage of the model, the steady seepage fields of the slope at water levels of 145 m and 175 m are obtained, which provides a theoretical basis for the subsequent dynamic analysis. See Figure 11 for calculation results of the seepage fields.Figure 11
Calculation results of the seepage fields.
(a)
Water level is 145 m
(b)
Water level is 175 m(3) Dynamic Response Calculation. Considering that the induced earthquake occurs later than the variation in water level, all plastic zones, displacement, and velocity vectors are reset to zero before loading of the earthquake. For the lateral boundary whose normal direction is in x-axis, due to the buffer effect of dampers, the coupling effect between the free-field mesh and the model can be codetermined by (a) the stress generated from the velocity difference between the free-field boundary node and the original model boundary node and (b) the nodal force of the free-field boundary node itself. This boundary condition can be expressed as follows: (10)
F
x
=
-
ρ
C
p
v
x
m
-
v
x
f
f
A
+
σ
x
x
f
f
Δ
S
y
,
F
y
=
-
ρ
C
s
v
y
m
-
v
y
f
f
A
+
σ
x
y
f
f
Δ
S
y
,
F
z
=
-
ρ
C
s
v
z
m
-
v
z
f
f
A
+
σ
x
z
f
f
Δ
S
y
,where ρ is the material density; C
p and C
s are the longitudinal and transverse wave velocity of the lateral boundary, respectively; A is the element area corresponding to F
x, F
y, and F
z being solved; v
x
m, v
y
m, and v
z
m  are the node velocity of the model boundary in x, y, and z directions, respectively; v
x
f
f, v
y
f
f, and v
z
f
f  are the node velocity of the free-field boundary in x, y, and z directions, respectively; σ
x
x
f
f, σ
x
y
f
f, and σ
x
z
f
f are normal stress in x
x direction, shear stress in x
y direction, and shear stress in x
z direction of the free-field mesh node, respectively; and Δ
S
y is the influence area of the free-field mesh node.


## 4.1. Establishment of Numerical Model

A numerical model is established based on the typical geological section in the main sliding direction of the landslide. The model is 1200 m in length and its rear edge has an elevation of 475 m. The controlling structural plane of the slope is simulated using a thin-layer element, while other joint fissures are taken into account by reducing the strength and deformation modulus of the rock mass. In order to ensure the accuracy of dynamic analysis, the size defined by the mesh should be controlled. The mesh size is limited by the shortest wavelength of the input seismic wave. To ensure the authenticity of wave propagation in the media, the maximum sizeΔ
l of the element must be less than 1/10~1/8 of the minimum wavelength. According to the elastic wave propagation theory, the propagation velocities of the longitudinal wave and transverse wave are calculated with the following formula: (7)
C
p
=
E
1
-
μ
ρ
1
+
μ
1
-
2
μ
C
s
=
E
2
ρ
1
+
μ
,where E is the elastic modulus, G is the shear modulus, ρ is the density, and μ is Poisson’s ratio of the material.From formula (7), the propagation velocities of the longitudinal wave and transverse wave in different rock strata are obtained. With the computation speed and the proposed frequency distribution of seismic waves taken into consideration, f
=
10 Hz is used to control mesh accuracy and the maximum mesh size can be determined by the following formula: (8)
Δ
l
=
1
10
·
C
s
f
.The final established numerical model is shown in Figure6. The distribution of rock is shown in Figure 7.Figure 6
Finite element model.Figure 7
Distribution of rock materials.The ideal elastic-plastic model is used. Assuming that the rock materials follow the Mohr-Coulomb strength criterion, the physical and mechanical parameters selected are shown in Table2.Table 2
Physical and mechanical parameters of numerical model.
Rock
ν
K (m·d−1)
E
d (GPa)
c (kPa)
c
s
a
t (kPa)
φ (°)
φ
s
a
t (°)
ρ (kg·m−3)
ρ
s
a
t (kg·m−3)
Sliding body
0.3
10
13
190
130
24
22
2000
2200
Weak layer
0.4
3
0.8
20
12
20
18
2100
2300
Silty mudstone
0.35
0.5
10
200
120
23
21
2300
2590
Argillaceous limestone
0.2
5
16
260
200
32
30
2550
2690
Limestone
0.18
—
32
400
350
35
32
2580
2680
ν is Poisson’s ratio, K is permeability coefficient, E
d is dynamic elastic modulus, c is cohesion, c
s
a
t is saturated cohesion, φ is internal friction angle, φ
s
a
t is saturated internal friction angle, ρ is natural density, and ρ
s
a
t is saturation density.


## 4.2. Seismic Wave Input

Reservoir-induced earthquakes are usually of low magnitude with concentrated epicenters and short duration (usually 5–10 s). Therefore, the time-history curve of a previously recorded real earthquake acceleration is used in this paper and the record of foreshocks is cut. The waveband in 7 seconds during 26 s~33 s is selected as the seismic waveform in this paper (see Figure8). Then, the seismic wave is reduced according to the peak acceleration of the corresponding magnitude. The acceleration time-history curve thus obtained is used for numerical calculation in the study.Figure 8
Original seismic wave time-history curve.
(a)
Acceleration time-history curve
(b)
Velocity time-history curve
(c)
Displacement time-history curveThe seismic wave shown in Figure8 cannot be used directly for numerical calculation. As can be seen from Figure 8(c), the final displacement of the seismic wave is nonzero after an earthquake cycle. It is not the result of residual deformation due to failure of objects, but the result of the incomplete waveform selected artificially or recording errors. Therefore, baseline calibration should be performed on the acceleration time-history curve to eliminate the influence of acceleration error on velocity and displacement, so that the final values of velocity and displacement time-history curves are zero. See Figure 9 for the time-history curve of seismic wave after filtering, baseline calibration, and reduction to Intensity IV (peak acceleration is 0.22 m/s2).Figure 9
Original seismic wave time-history curve (Intensity IV).
(a)
Acceleration time-history curve (Intensity IV)
(b)
Velocity time-history curve (Intensity IV)
(c)
Displacement time-history curve (Intensity IV)AsF
l
a
c
3
D software is used for calculation, two dynamic boundary conditions are provided, that is, the free-field boundary and viscous boundary. In this study, the viscous boundary and local damping (damping coefficient is 0.15) [30] are used for analysis. For this boundary condition, tangential and normal dampers are applied on the boundary to realize the absorption of incident wave energy. In other words, dampers produce tangential and normal force to offset the stress caused by reflected waves. Accordingly, the seismic wave should be input in the form of stress time-history. The expression for conversion of velocity time-history to stress time-history is as follows: (9)
σ
n
=
-
2
ρ
C
p
v
n
,
σ
s
=
-
2
ρ
C
s
v
s
,where σ
n is normal stress, σ
s is shear stress, ρ is medium density, C
p and C
s are the propagation velocity of P-wave and S-wave in the medium, respectively, and v
n and v
s are the velocity in vertical and horizontal directions, respectively.Reservoir-induced earthquakes of the same magnitude may occur multiple times successively within a short period of time. Given this characteristic, five earthquakes of Intensity IV (acceleration amplitude is 0.22 m/s2) are applied successively to simulate the frequent occurrence of microearthquakes.


## 4.3. Numerical Calculation

The following operating condition is mainly considered during numerical simulation: the reservoir water level rises to 175 m from 145 m. The calculation mainly involves dead weight, seepage, and dynamic response analysis.(1) Dead Weight Calculation. During dead weight calculation, gravity acceleration is set first, followed by full constraint on the bottom surface of the model and constraint on the normal displacement around the model. The stress field under gravity is obtained first. The nephogram of the vertical and horizontal stress fields of the slope under natural state is shown in Figure 10.Figure 10
Initial stress field nephogram of the slope under natural state.
(a)
Vertical and horizontal stress field
(b)
Horizontal stress field(2) Seepage Calculation. During seepage calculation, the deep rock mass is thick-layer and compact argillaceous limestone, and as its permeability is weaker than that of the upper rock mass, it can be regarded as a relative water-resisting layer of the slope. The front part of the slope is in direct contact with the reservoir water and thus its pore water pressure is set according to the actual water level. The rear part of the slope is less affected by water level fluctuation and thus the pore water pressure of the nodes can be set according to the groundwater depth revealed during drilling exploration. With bilateral seepage of the model, the steady seepage fields of the slope at water levels of 145 m and 175 m are obtained, which provides a theoretical basis for the subsequent dynamic analysis. See Figure 11 for calculation results of the seepage fields.Figure 11
Calculation results of the seepage fields.
(a)
Water level is 145 m
(b)
Water level is 175 m(3) Dynamic Response Calculation. Considering that the induced earthquake occurs later than the variation in water level, all plastic zones, displacement, and velocity vectors are reset to zero before loading of the earthquake. For the lateral boundary whose normal direction is in x-axis, due to the buffer effect of dampers, the coupling effect between the free-field mesh and the model can be codetermined by (a) the stress generated from the velocity difference between the free-field boundary node and the original model boundary node and (b) the nodal force of the free-field boundary node itself. This boundary condition can be expressed as follows: (10)
F
x
=
-
ρ
C
p
v
x
m
-
v
x
f
f
A
+
σ
x
x
f
f
Δ
S
y
,
F
y
=
-
ρ
C
s
v
y
m
-
v
y
f
f
A
+
σ
x
y
f
f
Δ
S
y
,
F
z
=
-
ρ
C
s
v
z
m
-
v
z
f
f
A
+
σ
x
z
f
f
Δ
S
y
,where ρ is the material density; C
p and C
s are the longitudinal and transverse wave velocity of the lateral boundary, respectively; A is the element area corresponding to F
x, F
y, and F
z being solved; v
x
m, v
y
m, and v
z
m  are the node velocity of the model boundary in x, y, and z directions, respectively; v
x
f
f, v
y
f
f, and v
z
f
f  are the node velocity of the free-field boundary in x, y, and z directions, respectively; σ
x
x
f
f, σ
x
y
f
f, and σ
x
z
f
f are normal stress in x
x direction, shear stress in x
y direction, and shear stress in x
z direction of the free-field mesh node, respectively; and Δ
S
y is the influence area of the free-field mesh node.


## 5. Dynamic Response Analysis

### 5.1. Slope Displacement Analysis

Figure12 shows the nephogram of slope displacement in X direction upon water level rise after each earthquake. Figure 13 shows the plastic state of the sliding mass corresponding to Figure 12. As can be seen from Figure 12, slope deformation mainly occurs in the potential sliding area, characterized by slide along the soft interlayer, while in other parts only minor deformation is found. Comparison among Figures 12(a) through 12(d) indicates that the displacement contour lines follow a similar distribution pattern; relatively large deformation is observed in the gentle slope platform in the front part of the potential sliding mass, and the horizontal displacement maximums appear at the position where the sliding surface is exposed to the free face. The displacement of the steeply dipping segment in the middle of the sliding mass, however, is reduced evidently compared with the front gentle slope platform. This indicates that the slope tends to be separated by shear failure along the soft interlayer, and the uplifted segment in the front supports the steeply dipping segment in the middle, prevents it from sliding, and limits its sliding space. Displacement of the rear gentle slope is relatively small. It decreases from the rear part to the front part, indicating the existence of traction from the steeply dipping segment in the middle and the possibility of sliding along the soft interlayer. Based on the accumulated displacement after five earthquake inputs, the distribution of deformation follows a certain pattern. To be specific, the accumulated horizontal displacement of the front-most part of the sliding mass reaches 24.65 mm; the accumulated horizontal displacement of the front gentle slope is around 18 mm; the accumulated horizontal displacement of the middle steep segment is 12–18 mm; and the accumulated horizontal displacement of the rock mass in the rear edge is around 8 mm. According to the horizontal displacement values, the displacement increases by the largest increment during the first earthquake and then by smaller increments in subsequent earthquakes, indicating a slower rate of increase.Figure 12
Nephogram of slope displacement inX direction.Figure 13
Plastic state of the sliding mass.From the plastic state of the potential sliding mass, tensile zones are mainly located at the front shear part and the rear gentle slope platform and shear zones at the front uplifted segment and the upper part of the rear gentle slope platform. No plastic yield is found in the rock stratum of the middle steeply dipping segment. On one hand, this indicates that the steeply dipping segment shows a general downslide trend under the action of microearthquakes and the rear gentle slope platform is relatively stable, leading to tensile stress in the rear edge of the steeply dipping segment and the yielding of many elements there by tension. As a result, tensile fractures may occur. Meanwhile, as the front gentile slope platform prevents sliding, thrust from the upper rock mass leads to damage of the upper and lower surfaces of the uplifted segment by shear. On the other hand, plastic deformation is found mainly in the middle and upper part of the rear gentle slope platform, while no shear yield is found in the bottom surface. This also indicates that sliding of the middle steeply dipping segment produces traction to the rear gentle slope platform which thus presents a downslide trend, resulting in continuous tensile zones in the rear edge of the slope. That is to say, the rear edge is under both tensile force and shear force.After the first earthquake, some elements of the potential sliding mass are still under shear (shear-n); after the second earthquake, all elements regain their elasticity and no connected plastic zones have ever been developed. It indicates that, under this operating condition, the slope remains in a stable state, regardless of local failure and a certain amount of permanent displacement.Under this operating condition, the water level rises so rapidly that the water pressure variation is not completely transferred into the slope, and the variation in pore water pressure is detected mainly on the surface part in contact with water. The huge additional water pressure functions like a presser foot to some extent. The seepage unbalance vector thus caused points to the interior of the slope and has little impact on slope stability.


### 5.2. Analysis of Main Sliding Zone

Figure14 shows the plastic state distribution of elements in the main sliding zone after five earthquakes. Shear failure has been found in most of the elements, but the sliding zone at the bottom surface of the rear gentle slope platform is never destroyed, which indicates that this part is of great slide-resistant capacity. Thus, there is relatively small possibility of sliding failure in positions above it. In case of overall failure of the slope, the middle steeply dipping segment first slides and deforms and then gradually overcomes the resistance of the front uplifted segment, resulting in overall instability. As the steeply dipping segment slides down, tensile stress occurs in its rear edge and large tensile fractures appear. The rock mass at the rear part is then subject to traction and may as a result slide down. After five earthquakes, elements in the main sliding zone regain their elasticity in the end, indicating that the slope is still stable under this operating condition.Figure 14
Plastic state of the main sliding zone after five earthquakes.The comparison above shows that deformation of the slope is controlled by the soft interlayer and is mainly manifested as bedding slide of the slope along the soft interlayer. Throughout the calculation process, no obvious sign of overall failure is found in the slope due to the following factors: (1) the main sliding zone has been adjusted after multiple earthquakes and thus has certain stability; (2) the main sliding zone is buried deep between the soft rock and block-layered rock, with no surface directly connected to the free face, and the uplifted segment in the front part plays a role of sliding prevention; and (3) microearthquakes affect slope stability mainly through long-term coupling with other factors.


### 5.3. Analysis of Slope Dynamic Response Amplitudes

From the analysis above, it is evident that the slope deformation mainly occurs during the first seismic wave input and is greatly affected by the priming effect. After entering the stable response phase, the displacement time-history curve repeats itself to a certain extent.F
l
a
c
3
D’s built-in programming language, FISH language, is used to record the acceleration and velocity amplitude of each node of the slope. The nephogram of amplitudes gives a general overview of dynamic response characteristics of the slope.Figure15 presents the acceleration amplitude nephogram of the slope after five earthquake inputs. It can be found that the acceleration amplitudes are generally larger in the main sliding zone. This is because the bedding of this zone is a substance differentiation plane where the media have relatively low elastic moduli. When the seismic wave propagates to this bedding, strong reflection occurs, followed by strong dynamic response that easily causes its further degradation. After the second earthquake, the acceleration amplitude of the slope becomes flat and the slope enters the stable response phase. From the distribution of acceleration amplitudes in the main sliding zone, the steeply dipping segment and the area below it have significantly larger acceleration amplitudes than the upper part, indicating that the upper rock mass is relatively stable and the lower rock mass has a downslide trend.Figure 15
Acceleration amplitude nephogram of the slope.Figure16 shows the nephogram of velocity amplitudes after five earthquake inputs. It can be found that velocity amplitudes are also larger in the main sliding zone, especially at the shear opening in the front edge where the maximum velocity amplitude appears, indicating that the slope tends to have shear failure along the main sliding zone. After the second earthquake, the distribution of velocity amplitudes becomes stable, with the larger values found near the main sliding zone and the free face. The distribution of velocity amplitudes in the main sliding zone is consistent with that of the acceleration amplitudes. All these indicate that the middle and lower soft interlayer have stronger dynamic response than the upper part and thus are more likely to experience degradation under long-term microearthquake action.Figure 16
Velocity amplitude nephogram of the slope.Under frequent earthquakes of Intensity IV induced by reservoir water level rise, the dynamic response of the slope has the following characteristics: (1) the soft interlayer absorbs and reflects much of the seismic wave. The reflected seismic wave is then superimposed with the incident wave, causing the development of an obvious amplification effect in the soft interlayer; as a result, the response amplitudes of the nodes in the soft interlayer are generally larger, easily leading to further degradation of the soft interlayer; (2) failure mode of the slope is closely related to its dynamic response characteristics; the dynamic response value changes greatly in possible failure positions; (3) under the operating condition of five earthquakes, the slope enters the stable response phase after a short adaptation period, which demonstrates that microearthquakes have limited influences on the slope.Figure17 shows the nephogram of shear strain increments and shear strain rates. It can be found that shear strain increments are concentrated in the middle of the soft interlayer and are small in values; the maximum appears at the front shear opening and no connected plastic zones are developed. Shear strain rates are higher in the steeply dipping segment than in other segments, but the overall level remains low and the slope is in a stable state, proving that the impact of frequent microearthquakes on the slope is finitely convergent.Figure 17
Nephogram of shear strain increments and shear strain rates of the slope after five earthquakes.
(a)
Final nephogram of shear strain increments
(b)
Final nephogram of shear strain ratesIn conclusion, after five earthquakes of Intensity IV upon reservoir water level rise, the slope deformation shows a convergent trend without obvious connection in the plastic zone. Its dynamic response becomes stable after a short adaptation period, indicating that microearthquakes have certain influences on slope stability.


### 5.4. Dynamic Safety Factor

The dynamic safety factor time-history curve of the slope after each earthquake can be obtained according to Section2.1. The section between the minimum safety factor and the maximum safety factor is discretized into a set of smaller equal intervals (Δ
d
=
0.01) and instantaneous safety factors of the slope at each moment are counted to obtain the number of safety factors in each discrete interval and thus the distribution probability in each discrete interval.As can be found from Figure18, safety factors of the slope are in an approximately normal distribution. Based on the basic method of mathematical statistics, the statistical parameter average μ
F and standard deviation σ
F of the probability distribution of safety factors after each microearthquake input can be determined. Based on formula (5) and the acceptable failure probability, the safety factor can be determined. See Table 3.Table 3
Safety factors under different numbers of earthquakes and failure probabilities.
P
f
n
n
=
1
n
=
2
n
=
3
n
=
4
n
=
5
0.05
1.085
1.073
1.062
1.053
1.045
0.01
1.074
1.061
1.051
1.042
1.034
0.001
1.060
1.049
1.038
1.029
1.021
0.0001
1.050
1.038
1.028
1.019
1.011
0.00001
1.041
1.0293
1.019
1.011
1.003
n is the number of intensity IV microearthquake inputs.Figure 18
Distribution of safety factors.
(a)
n
=
1
(b)
n
=
2
(c)
n
=
3
(d)
n
=
4
(e)
n
=
5Figure19 shows that, under the same probability of failure, as the number of earthquakes increases, the safety factor of the slope under microearthquakes of Intensity IV gradually decreases.Figure 19
Curves of the number of earthquake inputs and safety factors under different failure probabilities.


## 5.1. Slope Displacement Analysis

Figure12 shows the nephogram of slope displacement in X direction upon water level rise after each earthquake. Figure 13 shows the plastic state of the sliding mass corresponding to Figure 12. As can be seen from Figure 12, slope deformation mainly occurs in the potential sliding area, characterized by slide along the soft interlayer, while in other parts only minor deformation is found. Comparison among Figures 12(a) through 12(d) indicates that the displacement contour lines follow a similar distribution pattern; relatively large deformation is observed in the gentle slope platform in the front part of the potential sliding mass, and the horizontal displacement maximums appear at the position where the sliding surface is exposed to the free face. The displacement of the steeply dipping segment in the middle of the sliding mass, however, is reduced evidently compared with the front gentle slope platform. This indicates that the slope tends to be separated by shear failure along the soft interlayer, and the uplifted segment in the front supports the steeply dipping segment in the middle, prevents it from sliding, and limits its sliding space. Displacement of the rear gentle slope is relatively small. It decreases from the rear part to the front part, indicating the existence of traction from the steeply dipping segment in the middle and the possibility of sliding along the soft interlayer. Based on the accumulated displacement after five earthquake inputs, the distribution of deformation follows a certain pattern. To be specific, the accumulated horizontal displacement of the front-most part of the sliding mass reaches 24.65 mm; the accumulated horizontal displacement of the front gentle slope is around 18 mm; the accumulated horizontal displacement of the middle steep segment is 12–18 mm; and the accumulated horizontal displacement of the rock mass in the rear edge is around 8 mm. According to the horizontal displacement values, the displacement increases by the largest increment during the first earthquake and then by smaller increments in subsequent earthquakes, indicating a slower rate of increase.Figure 12
Nephogram of slope displacement inX direction.Figure 13
Plastic state of the sliding mass.From the plastic state of the potential sliding mass, tensile zones are mainly located at the front shear part and the rear gentle slope platform and shear zones at the front uplifted segment and the upper part of the rear gentle slope platform. No plastic yield is found in the rock stratum of the middle steeply dipping segment. On one hand, this indicates that the steeply dipping segment shows a general downslide trend under the action of microearthquakes and the rear gentle slope platform is relatively stable, leading to tensile stress in the rear edge of the steeply dipping segment and the yielding of many elements there by tension. As a result, tensile fractures may occur. Meanwhile, as the front gentile slope platform prevents sliding, thrust from the upper rock mass leads to damage of the upper and lower surfaces of the uplifted segment by shear. On the other hand, plastic deformation is found mainly in the middle and upper part of the rear gentle slope platform, while no shear yield is found in the bottom surface. This also indicates that sliding of the middle steeply dipping segment produces traction to the rear gentle slope platform which thus presents a downslide trend, resulting in continuous tensile zones in the rear edge of the slope. That is to say, the rear edge is under both tensile force and shear force.After the first earthquake, some elements of the potential sliding mass are still under shear (shear-n); after the second earthquake, all elements regain their elasticity and no connected plastic zones have ever been developed. It indicates that, under this operating condition, the slope remains in a stable state, regardless of local failure and a certain amount of permanent displacement.Under this operating condition, the water level rises so rapidly that the water pressure variation is not completely transferred into the slope, and the variation in pore water pressure is detected mainly on the surface part in contact with water. The huge additional water pressure functions like a presser foot to some extent. The seepage unbalance vector thus caused points to the interior of the slope and has little impact on slope stability.


## 5.2. Analysis of Main Sliding Zone

Figure14 shows the plastic state distribution of elements in the main sliding zone after five earthquakes. Shear failure has been found in most of the elements, but the sliding zone at the bottom surface of the rear gentle slope platform is never destroyed, which indicates that this part is of great slide-resistant capacity. Thus, there is relatively small possibility of sliding failure in positions above it. In case of overall failure of the slope, the middle steeply dipping segment first slides and deforms and then gradually overcomes the resistance of the front uplifted segment, resulting in overall instability. As the steeply dipping segment slides down, tensile stress occurs in its rear edge and large tensile fractures appear. The rock mass at the rear part is then subject to traction and may as a result slide down. After five earthquakes, elements in the main sliding zone regain their elasticity in the end, indicating that the slope is still stable under this operating condition.Figure 14
Plastic state of the main sliding zone after five earthquakes.The comparison above shows that deformation of the slope is controlled by the soft interlayer and is mainly manifested as bedding slide of the slope along the soft interlayer. Throughout the calculation process, no obvious sign of overall failure is found in the slope due to the following factors: (1) the main sliding zone has been adjusted after multiple earthquakes and thus has certain stability; (2) the main sliding zone is buried deep between the soft rock and block-layered rock, with no surface directly connected to the free face, and the uplifted segment in the front part plays a role of sliding prevention; and (3) microearthquakes affect slope stability mainly through long-term coupling with other factors.


## 5.3. Analysis of Slope Dynamic Response Amplitudes

From the analysis above, it is evident that the slope deformation mainly occurs during the first seismic wave input and is greatly affected by the priming effect. After entering the stable response phase, the displacement time-history curve repeats itself to a certain extent.F
l
a
c
3
D’s built-in programming language, FISH language, is used to record the acceleration and velocity amplitude of each node of the slope. The nephogram of amplitudes gives a general overview of dynamic response characteristics of the slope.Figure15 presents the acceleration amplitude nephogram of the slope after five earthquake inputs. It can be found that the acceleration amplitudes are generally larger in the main sliding zone. This is because the bedding of this zone is a substance differentiation plane where the media have relatively low elastic moduli. When the seismic wave propagates to this bedding, strong reflection occurs, followed by strong dynamic response that easily causes its further degradation. After the second earthquake, the acceleration amplitude of the slope becomes flat and the slope enters the stable response phase. From the distribution of acceleration amplitudes in the main sliding zone, the steeply dipping segment and the area below it have significantly larger acceleration amplitudes than the upper part, indicating that the upper rock mass is relatively stable and the lower rock mass has a downslide trend.Figure 15
Acceleration amplitude nephogram of the slope.Figure16 shows the nephogram of velocity amplitudes after five earthquake inputs. It can be found that velocity amplitudes are also larger in the main sliding zone, especially at the shear opening in the front edge where the maximum velocity amplitude appears, indicating that the slope tends to have shear failure along the main sliding zone. After the second earthquake, the distribution of velocity amplitudes becomes stable, with the larger values found near the main sliding zone and the free face. The distribution of velocity amplitudes in the main sliding zone is consistent with that of the acceleration amplitudes. All these indicate that the middle and lower soft interlayer have stronger dynamic response than the upper part and thus are more likely to experience degradation under long-term microearthquake action.Figure 16
Velocity amplitude nephogram of the slope.Under frequent earthquakes of Intensity IV induced by reservoir water level rise, the dynamic response of the slope has the following characteristics: (1) the soft interlayer absorbs and reflects much of the seismic wave. The reflected seismic wave is then superimposed with the incident wave, causing the development of an obvious amplification effect in the soft interlayer; as a result, the response amplitudes of the nodes in the soft interlayer are generally larger, easily leading to further degradation of the soft interlayer; (2) failure mode of the slope is closely related to its dynamic response characteristics; the dynamic response value changes greatly in possible failure positions; (3) under the operating condition of five earthquakes, the slope enters the stable response phase after a short adaptation period, which demonstrates that microearthquakes have limited influences on the slope.Figure17 shows the nephogram of shear strain increments and shear strain rates. It can be found that shear strain increments are concentrated in the middle of the soft interlayer and are small in values; the maximum appears at the front shear opening and no connected plastic zones are developed. Shear strain rates are higher in the steeply dipping segment than in other segments, but the overall level remains low and the slope is in a stable state, proving that the impact of frequent microearthquakes on the slope is finitely convergent.Figure 17
Nephogram of shear strain increments and shear strain rates of the slope after five earthquakes.
(a)
Final nephogram of shear strain increments
(b)
Final nephogram of shear strain ratesIn conclusion, after five earthquakes of Intensity IV upon reservoir water level rise, the slope deformation shows a convergent trend without obvious connection in the plastic zone. Its dynamic response becomes stable after a short adaptation period, indicating that microearthquakes have certain influences on slope stability.


## 5.4. Dynamic Safety Factor

The dynamic safety factor time-history curve of the slope after each earthquake can be obtained according to Section2.1. The section between the minimum safety factor and the maximum safety factor is discretized into a set of smaller equal intervals (Δ
d
=
0.01) and instantaneous safety factors of the slope at each moment are counted to obtain the number of safety factors in each discrete interval and thus the distribution probability in each discrete interval.As can be found from Figure18, safety factors of the slope are in an approximately normal distribution. Based on the basic method of mathematical statistics, the statistical parameter average μ
F and standard deviation σ
F of the probability distribution of safety factors after each microearthquake input can be determined. Based on formula (5) and the acceptable failure probability, the safety factor can be determined. See Table 3.Table 3
Safety factors under different numbers of earthquakes and failure probabilities.
P
f
n
n
=
1
n
=
2
n
=
3
n
=
4
n
=
5
0.05
1.085
1.073
1.062
1.053
1.045
0.01
1.074
1.061
1.051
1.042
1.034
0.001
1.060
1.049
1.038
1.029
1.021
0.0001
1.050
1.038
1.028
1.019
1.011
0.00001
1.041
1.0293
1.019
1.011
1.003
n is the number of intensity IV microearthquake inputs.Figure 18
Distribution of safety factors.
(a)
n
=
1
(b)
n
=
2
(c)
n
=
3
(d)
n
=
4
(e)
n
=
5Figure19 shows that, under the same probability of failure, as the number of earthquakes increases, the safety factor of the slope under microearthquakes of Intensity IV gradually decreases.Figure 19
Curves of the number of earthquake inputs and safety factors under different failure probabilities.


## 6. Conclusions

In this study, numerical simulation is conducted for frequent microearthquakes (Intensity IV) induced by reservoir water level rise at the rock bedded slope with a soft interlayer in the TGR area, and calculation results are analyzed based on the displacement deformation and dynamic response of the slope. The following conclusions are obtained:(1)
Although continuous microearthquake action has led to local failure and a certain amount of permanent displacement of the slope, the slope is eventually in a stable state. Under microearthquake action, the steeply dipping segment shows a general downslide trend, while the gentle slope platform at the rear edge is relatively stable. As a result, tensile stress develops at the rear edge of steeply dipping segment, leading to yielding of the rock mass in this position, which may further lead to the development of tensile fractures, an early sign of landslide.
(2)
The horizontal slip displacement of the slope increases from the rear to the front. Due to frequent microearthquakes, shear failure tends to develop along the soft interlayer. However, the uplifted segment in the front supports the steeply dipping segment in the middle, prevents it from sliding, and limits its sliding space. For slopes affected by frequent microearthquakes, appropriate antislide measures can effectively limit slope displacement.
(3)
Under the action of frequent microearthquakes, impact of the priming effect on slope deformation is the greatest and microearthquakes’ impact on slope deformation is finitely convergent.
(4)
Under microearthquake action, slope deformation is passive deformation caused by thrust load. It starts from the steeply dipping segment in the middle and is controlled by the soft interlayer, as can be seen from translational sliding of the slope along the soft interlayer. The soft interlayer absorbs much more seismic wave energy than the upper rock mass, which enables it to serve as a filter.
(5)
Under the same probability of failure, as the number of microearthquake inputs increases, the safety factor of the slope gradually decreases.Although the slope deformation under microearthquake action tends to be convergent, reservoir-induced earthquake is usually in close association with the operation cycle of the reservoir and is characterized by continuous occurrence. The cumulative deformation effect thus caused should receive high attention. In particular, earthquake monitoring and prediction shall be properly carried out according to the dynamic response and deformation mechanism of the slope under seismic action.



---

*Source: 7583108-2016-07-31.xml*