Received: 12 November 2020 Revised: 2 April 2021 Accepted: 2 April 2021

DOI: 10.1002/eqe.3462

RESEARCH ARTICLE

A generalized seismic sliding model of slopes with multiple

slip surfaces

Jian Song1 Adrian Rodriguez-Marek2 Tugen Feng3 Jian Ji4

1 Key Laboratory of Ministry of Education

for Geomechanics and Embankment Abstract

Engineering, Hohai University, Nanjing, Earthquake-induced permanent displacements of slopes are generally evaluated
China
through the simplified sliding block analyses of a singular coherent mass on a
2Department of Civil and Environmental
Engineering, Virginia Tech, Blacksburg,
predefined failure plane. However, coseismic landslides may exhibit multiple
Virginia slip surfaces as a result of earthquake shaking. This study presents a general-
3 Key Laboratory of Ministry of Education ized seismic sliding analysis using a column of sliding blocks to model seismic
for Geomechanics and Embankment
slope failures with multiple slip surfaces. The formulation of the sliding mode
Engineering, Hohai University, Nanjing,
China equations of the block system is obtained by a rigorous theoretical derivation.
4 Key Laboratory of Ministry of Education Both a rigid soil column associated with an acceleration–time series that is con-
for Geomechanics and Embankment stant with depth, and a flexible soil column considering an acceleration–time
Engineering, Hohai University, Nanjing,
China
series that varies throughout the sliding mass is presented. The conditional yield
accelerations and the equivalent seismic accelerations are defined to character-
Correspondence ize the complex interplay of the sliding block assemblies. The approach is applied
Tugen Feng, Key Laboratory of Min-
istry of Education for Geomechanics to model the seismic sliding behavior of layered slopes in which failures would
and Embankment Engineering, Hohai likely occur at the interfaces between soil layers with different shear strengths.
University, Nanjing 210098, China.
The effects of input ground motions and soil deposits on the sliding patterns of
Email: ftghhu@163.com
slopes are investigated. The developed generalized model provides a promising
Funding information means for evaluating the seismic slope stability by combining a number of seis-
National Natural Science Foundation of
mic failure locations within slopes.
China, Grant/Award Numbers: 41630638,
51879091; Fundamental Research Funds
for the Central Universities in China, KEYWORDS
Grant/Award Number: B210202046; dynamic response, earthquakes, landslides, multiple slip surfaces, sliding block, slope failure
National Key Research and Development
Project of China, Grant/Award Number:
2018YFC1508604

1 INTRODUCTION

Earthquake-induced landslides typically occur in seismically active areas that commonly exhibit steep topography, and
have generally been considered the most damaging hazard in numerous catastrophic earthquakes.1 Current procedures

Notation: s̈ i , sliding acceleration of the block; A, equivalent seismic acceleration–time history used for determination of the sliding initiation; a,
acceleration–time history or average acceleration–time history acting on the sliding mass; ag , input acceleration–time history; c, cohesion of the soil;
FN , normal force from the inclined plane; K, conditional yield coefficient of sliding block; k, yield coefficient of the sliding block; l, length of the
interface; m, mass of the sliding block; N, total number of blocks shown in Figure; n, considered number of blocks shown in Figure; ns , number of
block that is sliding; s, sliding displacement of the block; t, thickness of the sliding mass (or thickness of soil layer); Vs, shear wave velocity of the soil;
γ, unit weight of the soil; φ, friction angle of the soil

Earthquake Engng Struct Dyn. 2021;1–18. wileyonlinelibrary.com/journal/eqe © 2021 John Wiley & Sons Ltd. 1
2 SONG et al.

for assessing earthquake-induced landslides and slope failures are commonly based on the computation of permanent,
downslope displacement of the potential sliding mass during earthquake shaking. To this end, the sliding block method
has been the most commonly used methodology. Alternative approaches, such as advanced numerical analyses with non-
linear constitutive models to represent soil behavior, use multiple parameters that can be difficult to constrain and are
computationally costly. Moreover, sliding block methods have been demonstrated to accurately predict the order of mag-
nitude of observed displacements in earthquakes,2,3 and their simplicity allows for the consideration of the influence of
ground motion variability in the computed displacements (e.g., Refs. 4 and 5).
The original sliding block method was proposed by Newmark in the 1965 Rankine Lecture.6 In this method, a rigid
block sliding on a plane is used as an analog to the earth mass sliding over a shear surface. The shaking accelerations
produce destabilizing inertial forces on the rigid sliding mass that can episodically exceed the yield resistance along the
slip surface and induce sliding. Once the sliding is initiated, the acceleration of the block is different from the acceleration
of the ground. The resulting relative sliding displacements continue until the relative sliding velocity equals zero and the
block and underlying ground move together. The displacements of all sliding episodes add up to represent the permanent
downslope displacement at the end of the earthquake-induced shaking.
After its original formulation, some modifications have been proposed to the Newmark sliding block model by relax-
ing the assumption of a rigid sliding mass and considering a flexible sliding mass. These modifications introduce the
dynamic responses of the sliding mass via a coupled or decoupled manner. The coupled approaches have been developed
by simultaneously computing the dynamic and sliding responses of the sliding mass.7–10 The coupled approach rigor-
ously accounts for the interaction of the soil deformability and the sliding based on the stick-slip model. Alternatively,
the decoupled approach, first introduced by Makdisi and Seed,11 is accomplished by separately performing the dynamic
response analysis and the displacement computation of the sliding mass. One- (1D) or two-dimensional (2D) dynamic
response analyses are commonly adopted to compute the average seismic acceleration experienced by the potential slide
mass.12 This average acceleration–time history including the effects of slope response is then used as the input motion in
the Newmark rigid block method to compute the sliding displacements. Although the simplifying assumption that the
dynamic response and sliding are mutually independent, the decoupled approach provides acceptable predictions in cer-
tain scenarios and is easy to apply with the current 1D seismic site response and numerical stress–deformation computer
codes, and hence, it is still widely used (e.g., Refs. 13–15).
The aforementioned Newmark-type deformation analyses are based on an assumption that displacements occur along a
single, well-defined slip surface. However, earthquake-induced landslides may exhibit multiple slip surfaces or distributed
shear deformations.16–18 This failure pattern, with the development of a set of sliding surfaces, could be significant for
slopes of heterogeneous soil layers with distinctly different shear strengths or slopes with some preexisting weak planes.
However, there are very few studies considering the presence of multiple shear movement.
In this paper, we present an analysis of the progressive sliding pattern associated with multiple slip surfaces under
earthquake loading, and propose a generalized seismic sliding model for slopes with the presence of multiple slip sur-
faces. Three aspects of significant improvements are introduced in this study compared to the previous studies, namely,
the possibility in the development of coseismic multiple slip surfaces of slopes is demonstrated from the numerical anal-
ysis, a generalized solution of different seismic sliding modes is innovatively derived for a generic number of sliding
blocks by a theoretical analysis, and this generalized model is extended to flexible sliding systems by introducing the
dynamic site responses. The complex interplay of the system of sliding masses is explicitly resolved. The paper is orga-
nized as follows. First, an analysis to show how multiple failure surfaces can develop is presented, followed by the devel-
opment of the proposed methodology. Thereafter, the theoretical solutions of multiblock seismic sliding analyses are ver-
ified with those obtained from numerical analyses. Finally, the methodology is applied to slopes with heterogeneous soil
layers.

2 DEVELOPMENT OF MULTIPLE FAILURE SURFACES

The multiple failure surfaces or regions of local shear deformations of earthquake-induced landslides have been observed,
albeit rarely, in the previous studies. For example, by performing the shaking table test of earthquake-induced deforma-
tions of earth slopes, Wartman et al17 found that the deformation mainly occurred along two or more localized slip surfaces
of the same general orientation for some clayey slope models. In our study, numerical simulations were first performed
to further demonstrate the development of multiple slip surfaces associated with coseismic slope failures. A slope model
consisting of a homogenous soil layer and another slope model with a soft overlay were considered. The slope geometry
SONG et al. 3

F I G U R E 1 (A) Homogenous soil slope and (B) layered soil slope with a soft overlay (γ, Vs, c, and φ are the unit weight, shear wave
velocity, cohesion, and friction angle of soils, respectively)

F I G U R E 2 Shear–strain increment from dynamic analysis of the layered slope with a soft overlay for ground motions with different
linear scaling factors: (A) scaling factor = 0.25 and (B) scaling factor = 1.0

and the soil parameters are shown in Figure 1. The fully nonlinear dynamic analyses were conducted by using the 2D finite
difference code (FLAC19 ) with a constitutive model combining hysteretic and plastic behavior of soil. The shear modulus
reduction and damping versus shear strain curves were based on the standard one proposed by Sun et al20 for clay, and the
Mohr–Coulomb model was utilized to describe the plastic behavior. The width of the slope model was 400 m. The lateral
boundaries were set as the free-field condition, and a rigid bedrock was assumed at the bottom boundary. The ground
motion input was the 230◦ horizontal component of the El Centro Array #6 station recorded in the 1979 Imperial Valley-
06 earthquake. Element size dictates the range of input motion frequencies that can be captured in dynamic analyses. To
capture a broad band of frequencies, a maximum frequency of 20 Hz was considered, which corresponds to a minimum
wavelength of 20 and 10 m for Soil 1 and Soil 2, respectively. The element size was selected to be 2 and 1 m for Soil 1 and
Soil 2 (a tenth of the minimum wavelength).
The slip surfaces for the two slopes at the end of earthquake shaking are shown in Figures 2 and 3. The position of
these slip surfaces was determined based on the distribution of shear–strain increments. The results of dynamic analysis
for input motions linearly scaled to different amplitudes are also presented to exemplify the slope deformations with
the development of nonunique failure surfaces. Observe from Figure 2 that the initial failure of the layered slope with a
soft overlay is constrained within the surficial soil layer, resulting in local-failure mechanism. When the ground motion
becomes stronger, a deep failure is generated and the sliding displacement occurs along both slip surfaces during the
earthquake. For the homogenous soil slope (Figure 3), the failure surfaces pass through the toe for different ground motion
intensities, and a stronger earthquake shaking only produces larger shear–strain increments and thus larger permanent
displacements. The absence of the shallow failure after the deep failure in the homogenous soil slope is likely the result of

F I G U R E 3 Shear–strain increment from dynamic analysis of the slope with a uniform soil layer for ground motions with different linear
scaling factors: (A) scaling factor = 1.0 and (B) scaling factor = 1.6
4 SONG et al.

F I G U R E 4 Schematics of forces for the block at the limit equilibrium state and the sliding state under an earthquake motion parallel to
the slope (A and B) and under a horizontal earthquake motion (C and D)

the “isolation effects” of inertial forces transmitted to the upper slope bodies due to the sliding or permanent deformation
along the slip surface. This would be further discussed in the following presentation of block analysis.

3 METHODOLOGY

The numerical results presented above have demonstrated the possibility of the occurrence of coupled shallow and deep
failure modes associated with earthquake-induced landslides. A generalized seismic sliding model is built conceptually
on the Newmark method by considering multiple slip surfaces. The model discretizes the slope into a generic number of
sliding blocks, and the slip episodes of the block system are rigorously resolved for uniform and distributed acceleration–
time histories throughout the sliding masses under earthquake loading. The uniform acceleration–time histories represent
the rigid multiblock system, and the distributed acceleration–time histories describe the effects of dynamic response of
the flexible multiblock system. The proposed approach is more appropriate for translational failures, which is a common
failure mechanism to sliding block methodologies. This is a typical failure mode of infinite slope or dip slope conditions,21
and is different from the rotational sliding model of seismic slope failures that considers the geometry changes of the
sliding mass.22

3.1 Uniform acceleration–time series throughout the sliding masses

Before the model development, the yield acceleration (which is the yield coefficient multiplied by the acceleration of
gravity g) in different directions is first introduced. For a single block on an inclined plane, the yield acceleration is the
base acceleration that results in the block reaching a state of limit equilibrium (i.e., a factor of safety equal to 1). In the
previous studies, the base acceleration is usually considered to act parallel to the slope (Figure 4A). In this case, the yield
coefficient parallel to the slope (ks ) can be obtained as

𝑎𝑔 𝑐
𝑘𝑠 = − = + (𝑐𝑜𝑠𝛼 ⋅ 𝑡𝑎𝑛𝜑 − 𝑠𝑖𝑛𝛼), (1)
𝑔 𝛾⋅𝑡

where ag is the input acceleration; α is the inclination of the slope; c and φ are the cohesion and friction angle of the
interface; and γ and t are the unit weight and the thickness of the block. In the figure, m is the mass of the block, FN is
the normal force, and l is length of the interface. During the sliding condition under the input acceleration parallel to the
slope (Figure 4B), the relative sliding acceleration of the block along the slope direction can be obtained as

𝑠̈𝑠 = −𝑎𝑔 − 𝑘𝑠 ⋅ 𝑔. (2)

SONG et al. 5

FIGURE 5 Values of term (cosα+sinαtanφ)cosα for different slope angles (α) and friction angles of soil (φ)

FIGURE 6 Rigid block assemblies with N blocks when none of the block slides

Similarly, for the block under a horizontal inertial force (Figure 4C), the yield coefficient in the horizontal direction
(kh ) can also be obtained as

𝑎𝑔 𝑐
𝑘ℎ = − = + tan(𝜑 − 𝛼), (3)
𝑔 𝛾 ⋅ 𝑡 ⋅ (𝑐𝑜𝑠𝛼 + 𝑠𝑖𝑛𝛼 ⋅ 𝑡𝑎𝑛𝜑)

and the sliding acceleration of the block in the horizontal direction relative to the base (Figure 4D) can be obtained as
( )
𝑠̈ℎ = 𝑠̈𝑠 ⋅ cos 𝛼 = −𝑎𝑔 − 𝑘ℎ ⋅ 𝑔 ⋅ [(cos 𝛼 + sin 𝛼 ⋅ tan 𝜑) ⋅ cos 𝛼] . (4)

The value of term [(cos 𝛼 + sin 𝛼 ⋅ tan 𝜑) ⋅ cos 𝛼] is generally close to unity for a practical range of values of the param-
eters (α and φ), as illustrated in Figure 5. The value of this term is assumed equal to unity in this paper. Therefore, in the
following model development, the equations are presented in the unified form 𝑠̈ = −𝑎𝑔 − 𝑘 ⋅ 𝑔. When using these solu-
tions, if the input motion is assumed to be parallel to the slope, then the yield coefficient parallel to the slope should be
used, and the calculated result is the sliding displacement parallel to the slope. If the input motion is considered to be
horizontal, then the horizontal yield coefficient should also be used, and the calculated result would be the horizontal
sliding displacement.
Consider an infinite slope with N discrete sliding blocks. At the initial stage when none of the blocks slides (Figure 6),
each block may potentially move along the slip surface in the case that the input acceleration is larger than the correspond-
ing yield acceleration. This is the same as the traditional Newmark approach. The critical condition for the initiation of
sliding of block i is
[ 𝑖
] [ 𝑖
]
∑ ∑
− 𝑚𝑏 ⋅ 𝑎𝑔 = 𝑘𝑖 ⋅ 𝑚𝑏 ⋅ 𝑔, −𝑎𝑔 = 𝑘𝑖 ⋅ 𝑔, (5a)
𝑏=1 𝑏=1
6 SONG et al.

FIGURE 7 Rigid block assemblies with N blocks when one of the blocks slides: (A) block i2 slides first and (B) block i1 slides first

and after sliding, the sliding acceleration of block i can be obtained based on the equilibrium at the slip surface as
[ 𝑖
] [ 𝑖 ]
∑ ( ) ∑
− 𝑚𝑏 ⋅ 𝑎𝑔 + 𝑠̈𝑖 = 𝑘𝑖 ⋅ 𝑚𝑏 ⋅ 𝑔, 𝑠̈𝑖 = −𝑎𝑔 − 𝑘𝑖 ⋅ 𝑔, (5b)
𝑏=1 𝑏=1

where ki and s̈ i represent the yield coefficient and sliding acceleration of block i, respectively; m is the mass of each block.
For the case when the acceleration–time series are uniform throughout the sliding masses, the initial sliding occurs on
the block with the minimum value of yield coefficient (k).
After initial sliding of the block with the minimum k, the state of the block system is shown in Figure 7. For the upper
block i1 and the lower block i2 , if block i2 slides first (Figure 7A, i.e., ki2 < ki1 ), the governing equation of equilibrium on
the slip surface below block i2 and the sliding acceleration of block i2 are given by
[ 𝑖2
] [𝑖 ]
∑ ( ) ∑2
Block# 𝑖2 ∶− 𝑚𝑏 ⋅ 𝑎𝑔 + 𝑠̈𝑖2 = 𝑘𝑖2 ⋅ 𝑚𝑏 ⋅ 𝑔, 𝑠̈𝑖2 = −𝑎𝑔 − 𝑘𝑖2 ⋅ 𝑔. (6)
𝑏=1 𝑏=1

Then, the total inertial force above the slip surface below block i1 can be obtained as
[ 𝑖1
] [𝑖 ] [𝑖 ]
∑ ( ) ∑1 ∑1
− 𝑚𝑏 ⋅ 𝑎𝑔 + 𝑠̈𝑖2 = 𝑘𝑖2 ⋅ 𝑚𝑏 ⋅ 𝑔 < 𝑘𝑖1 ⋅ 𝑚𝑏 ⋅ 𝑔. (7)
𝑏=1 𝑏=1 𝑏=1

Therefore, relative sliding would not occur for block i1 .

If block i1 slides first (Figure 7B, i.e., ki1 < ki2 ), the sliding acceleration of block i1 and critical condition for the sliding
initiation of block i2 are given by
[ 𝑖1
] [𝑖 ]
∑ ( ) ∑1
Block# 𝑖1 ∶ − 𝑚𝑏 ⋅ 𝑎𝑔 + 𝑠̈𝑖1 = 𝑘𝑖1 ⋅ 𝑚𝑏 ⋅ 𝑔, 𝑠̈𝑖1 = −𝑎𝑔 − 𝑘𝑖1 ⋅ 𝑔, (8a)
𝑏=1 𝑏=1
[ 𝑖2
] [ 𝑖1
] [ 𝑖2
]
∑ ∑ ∑
Block# 𝑖2 ∶ − 𝑚𝑏 ⋅ 𝑎𝑔 + 𝑘𝑖1 ⋅ 𝑚𝑏 ⋅ 𝑔 = 𝑘𝑖2 ⋅ 𝑚𝑏 ⋅ 𝑔, −𝑎𝑔 = 𝐾𝑖2 −𝑖1 ⋅ 𝑔, (8b)
𝑏=𝑖1 +1 𝑏=1 𝑏=1

where Ki2 −i1 represents the conditional yield coefficient of block i2 under the premise of the sliding occurrence of block
i1 . This is defined as

∑𝑖1
𝑚𝑏 ( )
𝐾𝑖2 −𝑖1 =𝑘𝑖2 + ∑𝑖 𝑏=1 ⋅ 𝑘𝑖2 − 𝑘𝑖1 . (9)
2
𝑏=𝑖1 +1 𝑚𝑏

A more complex and generic state of the block system is to consider n of the N blocks, where n−1 blocks slide simultane-
ously, as shown in Figure 8. The three subplots in the figure describe the different locations of the block, whose sliding has
SONG et al. 7

FIGURE 8 Rigid block assemblies with N blocks when n−1 of the blocks slide

not yet begun, in the n blocks. For the case shown in Figure 8A, the governing equations of equilibrium on the different
slip surfaces are given by
[ 𝑖2
] [ 𝑖2
]
∑ ( ) ∑
Block# 𝑖2 ∶ − 𝑚𝑏 ⋅ 𝑎𝑔 + 𝑠̈𝑖2 + +𝑠̈𝑖𝑛 = 𝑘𝑖2 ⋅ 𝑚𝑏 ⋅ 𝑔, (10a)
𝑏=1 𝑏=1
[ 𝑖3
] [𝑖 ] [𝑖 ]
∑ ( ) ∑2 ∑3
Block# 𝑖3 ∶− 𝑚𝑏 ⋅ 𝑎𝑔 + 𝑠̈𝑖3 + +𝑠̈𝑖𝑛 + 𝑘𝑖2 ⋅ 𝑚𝑏 ⋅ 𝑔 = 𝑘𝑖3 ⋅ 𝑚𝑏 ⋅ 𝑔, (10b)
𝑏=𝑖2 +1 𝑏=1 𝑏=1

⎡ ∑ 𝑖𝑗 ⎤ ( ) ⎡∑𝑖𝑗−1 ⎤ ⎡∑𝑖𝑗 ⎤
Block# 𝑖𝑗 ∶ − ⎢ 𝑚𝑏 ⎥ ⋅ 𝑎𝑔 + 𝑠̈𝑖𝑗 + +𝑠̈𝑖𝑛 + 𝑘𝑖𝑗−1 ⋅ ⎢ 𝑚𝑏 ⎥ ⋅ 𝑔 = 𝑘𝑖𝑗 ⋅ ⎢ 𝑚𝑏 ⎥ ⋅ 𝑔, (10c)
⎢𝑏=𝑖 +1 ⎥ ⎢𝑏=1 ⎥ ⎢𝑏=1 ⎥
⎣ 𝑗−1 ⎦ ⎣ ⎦ ⎣ ⎦
⎡ ∑
𝑖𝑗+1 ⎤ ( ) ⎡∑𝑖𝑗 ⎤ ⎡∑𝑖𝑗+1 ⎤
Block# 𝑖𝑗 +1 ∶ − ⎢ 𝑚𝑏 ⎥ ⋅ 𝑎𝑔 + 𝑠̈𝑖𝑗+1 + +𝑠̈𝑖𝑛 + 𝑘𝑖𝑗 ⋅ ⎢ 𝑚𝑏 ⎥ ⋅ 𝑔 = 𝑘𝑖𝑗+1 ⋅ ⎢ 𝑚𝑏 ⎥ ⋅ 𝑔, (10d)
⎢𝑏=𝑖 +1 ⎥ ⎢𝑏=1 ⎥ ⎢𝑏=1 ⎥
⎣ 𝑗 ⎦ ⎣ ⎦ ⎣ ⎦
[ 𝑖 ] [𝑖 ] [𝑖 ]
∑
𝑛−1
( ) ∑
𝑛−2 ∑𝑛−1
#
Block 𝑖𝑛 −1 ∶ − 𝑚𝑏 ⋅ 𝑎𝑔 + 𝑠̈𝑖𝑛−1 + 𝑠̈𝑖𝑛 + 𝑘𝑖𝑛−2 ⋅ 𝑚𝑏 ⋅ 𝑔 = 𝑘𝑖𝑛−1 ⋅ 𝑚𝑏 ⋅ 𝑔, (10e)
𝑏=𝑖𝑛−2 +1 𝑏=1 𝑏=1
[ 𝑖𝑛
] [𝑖 ] [𝑖 ]
∑ ( ) ∑𝑛−1 ∑𝑛
Block# 𝑖𝑛 ∶ − 𝑚𝑏 ⋅ 𝑎𝑔 + 𝑠̈𝑖𝑛 + 𝑘𝑖𝑛−1 ⋅ 𝑚𝑏 ⋅ 𝑔 = 𝑘𝑖𝑛 ⋅ 𝑚𝑏 ⋅ 𝑔. (10f)
𝑏=𝑖𝑛−1 +1 𝑏=1 𝑏=1

By solving Equations (10a)–(10f), the sliding accelerations of blocks are obtained as

𝑠̈𝑖2 = 𝐾𝑖3 − 𝑖 ⋅ 𝑔 − 𝑘𝑖2 ⋅ 𝑔, (11a)

2

𝑠̈𝑖𝑗 = 𝐾𝑖𝑗+1 − ⋅ 𝑔 − 𝐾𝑖𝑗 −𝑖𝑗−1 ⋅ 𝑔, (11b)

𝑖𝑗

𝑠̈𝑖𝑛−1 = 𝐾𝑖𝑛 − 𝑖 ⋅ 𝑔 − 𝐾𝑖𝑛−1 −𝑖𝑛−2 ⋅ 𝑔, (11c)

𝑛−1

𝑠̈𝑖𝑛 = −𝑎𝑔 − 𝐾𝑖𝑛 − 𝑖 ⋅ 𝑔, (11d)

𝑛−1

where K is the conditional yield coefficient of different blocks, and is defined similar to Equation (9). Also, block i1 would
be locked and never slide because the total inertial force above the slip surface below block i1 is consistently less than the
yield resistance, as shown in Equation (7).
8 SONG et al.

T A B L E 1 Conditional yield coefficient and equivalent seismic acceleration used for determination of the sliding initiation under
uniform acceleration–time series
Equivalent seismic
ns Conditional yield coefficient acceleration
ns = 0 For block i, Ki = ki 𝐴 = −𝑎𝑔
ns > 0 For block i1 , i2 ,. . . , ij-1 , ij , ij+1 ,. . . , in-1 , in (i1 < i2 < . . . < in ), 𝐴 = −𝑎𝑔
𝐾𝑖𝑗 −𝑖2 …𝑖𝑛 = +∞, 𝑖𝑗 = 𝑖1
𝐾𝑖𝑗 −𝑖1 𝑖𝑛 = +∞, 𝑖1 < 𝑖𝑗 < 𝑖𝑛
𝑖2 …𝑖𝑗−1 𝑖𝑗+1 …
∑𝑖𝑛−1
𝑚𝑏
𝐾𝑖𝑗 −𝑖1 𝑖2 …𝑖𝑛−1 = 𝐾𝑖𝑗 − = 𝑘𝑖𝑛 + ∑𝑖 𝑛
𝑏=1
⋅ (𝑘𝑖𝑛 − 𝑘𝑖𝑛−1 ), 𝑖𝑗 = 𝑖𝑛
𝑖𝑛−1 𝑏=𝑖𝑛−1 +1
𝑚𝑏

Note: ns is the number of block that is sliding; the block number is in ascending order from the top to the bottom of the block system.

TA B L E 2 Sliding acceleration of blocks in different conditions under uniform acceleration–time series

ns Sliding acceleration
ns = 1 For block i, 𝑠̈𝑖 = −𝑎𝑔 − 𝑘𝑖 ⋅ 𝑔
ns > 1 For block i1 , i2 ,. . . , ij-1 , ij , ij+1 ,. . . , in-1 , in (i1 < i2 < . . . < in )
𝑠̈𝑖1 = 𝐾𝑖2 −𝑖1 ⋅ 𝑔 − 𝑘𝑖1 ⋅ 𝑔
𝑠̈𝑖𝑗 = 𝐾𝑖𝑗+1 −𝑖𝑗 ⋅ 𝑔 − 𝐾𝑖𝑗 −𝑖𝑗−1 ⋅ 𝑔
𝑠̈𝑖𝑛 = −𝑎𝑔 − 𝐾𝑖𝑛 −𝑖𝑛−1 ⋅ 𝑔

Similar derivation can be performed for the cases shown in Figure 8B and C. The sliding accelerations of blocks shown
in Figure 8B are obtained as

𝑠̈𝑖1 = 𝐾𝑖2 − 𝑖 ⋅ 𝑔 − 𝑘𝑖1 ⋅ 𝑔, (12a)

1

𝑠̈𝑖2 = 𝐾𝑖3 − 𝑖 ⋅ 𝑔 − 𝐾𝑖2 −𝑖1 ⋅ 𝑔, (12b)

2

𝑠̈𝑖𝑗 ≡ 0, (12c)

𝑠̈𝑖𝑛−1 = 𝐾𝑖𝑛 − 𝑖 ⋅ 𝑔 − 𝐾𝑖𝑛−1 −𝑖𝑛−2 ⋅ 𝑔, (12d)

𝑛−1

𝑠̈𝑖𝑛 = −𝑎𝑔 − 𝐾𝑖𝑛 − 𝑖 ⋅ 𝑔, (12e)

𝑛−1

and the sliding accelerations of blocks and the condition for the sliding initiation of block in shown in Figure 8C are
obtained as

𝑠̈𝑖1 = 𝐾𝑖2 − 𝑖 ⋅ 𝑔 − 𝑘𝑖1 ⋅ 𝑔, (13a)

1

𝑠̈𝑖2 = 𝐾𝑖3 − 𝑖 ⋅ 𝑔 − 𝐾𝑖2 −𝑖1 ⋅ 𝑔, (13b)

2

𝑠̈𝑖𝑗 = 𝐾𝑖𝑗+1 − ⋅ 𝑔 − 𝐾𝑖𝑗 −𝑖𝑗−1 ⋅ 𝑔, (13c)

𝑖𝑗
𝑠̈𝑖𝑛−1 = −𝑎𝑔 − 𝐾𝑖𝑛−1 − 𝑖 ⋅ 𝑔, (13d)
𝑛−2

Block# 𝑖𝑛 ∶ −𝑎𝑔 = 𝐾𝑖𝑛 − 𝑖 ⋅ 𝑔. (13e)

𝑛−1

As is evident from the previous discussion, the sliding process of the block system is much more complex than the
conventional single-block sliding. The conditions of sliding initiation and the sliding acceleration of each block interact
because of the “isolation effects” on inertial forces. For the uniform acceleration–time series (i.e., constant acceleration
throughout the sliding masses), the sliding of the lower block (i.e., a smaller yield acceleration) results in a constant
force transmitted to the upper block (with a larger yield acceleration), which is consistently less than the yield resistance,
and thus the upper block would be locked. In addition, an initial sliding of the upper block would restrict the force
transmitted to the lower block to its yield resistance, and thus, the larger inertial forces are isolated by the sliding of the
upper block. In spite of the complexity of the sliding block system, the slip episodes can be resolved with the explicit
description of the initiation of sliding and the sliding accelerations in different cases. These conditions, which follow
from the above derivations, are summarized in Tables 1 and 2. The initiation of sliding corresponds to conditions when
SONG et al. 9

the equivalent seismic acceleration (A) becomes larger than the corresponding conditional yield acceleration (K⋅g). For
rigid-block system above, the equivalent seismic acceleration (A) is the input acceleration, but it would be different for
the following flexible-block system.

3.2 Distributed acceleration–time series throughout the sliding masses

For the compliant slope in which the dynamic response is significant during earthquake shaking, the effects of dynamic
response of slopes on the inertial force acting on the sliding block system should be considered to provide more realistic
input motions for the sliding block analysis. To this end, distributed acceleration–time histories throughout the block
system are implemented in the multiblock seismic sliding model. These distributed acceleration–time histories describe
the average dynamic response of each sliding mass and the inertial force at the potential failure surfaces. Both 2D and
simplified 1D analyses can be employed to produce the time-varying distributed accelerations.
Under the distributed inertial forces, each sliding mass may potentially move along the slip surface. At the initial stage
when none of the blocks slide (Figure 6), the critical condition for the initiation of sliding and the sliding acceleration of
block i are:

𝑖
[ 𝑖
]
∑ ∑
− (𝑚𝑏 ⋅ 𝑎𝑏 ) = 𝑘𝑖 ⋅ 𝑚𝑏 ⋅ 𝑔, −𝑎1𝑖 = 𝑘𝑖 ⋅ 𝑔, (14a)
𝑏=1 𝑏=1
𝑖
[ 𝑖
]
∑ ∑
− [𝑚𝑏 ⋅ (𝑎𝑏 + 𝑠̈𝑖 )] = 𝑘𝑖 ⋅ 𝑚𝑏 ⋅ 𝑔, 𝑠̈𝑖 = −𝑎1…𝑖 − 𝑘𝑖 ⋅ 𝑔, (14b)
𝑏=1 𝑏=1

where a is the acceleration–time history of each block; a1. . . i is the average acceleration–time history acting on the sliding
mass (from block 1 to block i) and is defined as
∑𝑖
𝑏=1 (𝑚𝑏 ⋅ 𝑎𝑏 )
𝑎1⋅⋅⋅𝑖 = ∑𝑖 . (15)
𝑚
𝑏=1 𝑏

Contrary to the condition of uniform accelerations, the block with the minimum yield acceleration is not necessarily
the first to slide under the distributed acceleration–time series.
After the occurrence of initial sliding of a block (Figure 7), the governing equations of equilibrium on the slip surface
below block i2 in Figure 7A are given by

𝑖2
[𝑖 ]
∑ [ ( )] ∑2
Block# 𝑖2 ∶− 𝑚𝑏 ⋅ 𝑎𝑏 + 𝑠̈𝑖2 = 𝑘𝑖2 ⋅ 𝑚𝑏 ⋅ 𝑔. (16)
𝑏=1 𝑏=1

The sliding acceleration of block i2 and the condition of sliding along the potential slip surface below block i1 are derived
as

𝑠̈𝑖2 = −𝑎1…𝑖2 − 𝑘𝑖2 ⋅ 𝑔, (17a)

𝑖1
[𝑖 ]
∑ [ ( )] ∑1
Block# 𝑖1 ∶ − 𝑚𝑏 ⋅ 𝑎𝑏 + 𝑠̈𝑖2 = 𝑘𝑖1 ⋅ 𝑚𝑏 ⋅ 𝑔, −𝑎1…𝑖1 − (−𝑎1…𝑖2 ) = 𝑘𝑖1 ⋅ 𝑔 − 𝑘𝑖2 ⋅ 𝑔. (17b)
𝑏=1 𝑏=1

The sliding acceleration of block i1 and critical condition for the sliding along the potential slip surface below block i2
shown in Figure 7B are given by

𝑖1
[𝑖 ]
∑ [ ( )] ∑1
Block# 𝑖1 ∶− 𝑚𝑏 ⋅ 𝑎𝑏 + 𝑠̈𝑖1 = 𝑘𝑖1 ⋅ 𝑚𝑏 ⋅ 𝑔, 𝑠̈𝑖1 = −𝑎1…𝑖1 − 𝑘𝑖1 ⋅ 𝑔, (18a)
𝑏=1 𝑏=1
10 SONG et al.

T A B L E 3 Conditional yield coefficient and equivalent seismic acceleration used for determination of the sliding initiation under
distributed acceleration–time series
ns Conditional yield coefficient Equivalent seismic acceleration
ns = 0 For block i, Ki = ki 𝐴𝑖 = −𝑎1…𝑖
ns > 0 For block i1 , i2 ,. . . , ij-1 , ij , ij+1 ,. . . , in-1 , in 𝐴𝑖𝑗 −𝑖2 …𝑖𝑛 = −𝑎1…𝑖1 − (−𝑎1…𝑖2 ), 𝑖𝑗 = 𝑖1
(i1 < i2 < . . . < in ), 𝐴𝑖𝑗 −𝑖1 𝑖2 …𝑖𝑗−1 𝑖𝑗+1 … 𝑖𝑛 = −𝑎(𝑖𝑗−1 +1)…𝑖𝑗 − [−𝑎(𝑖𝑗−1 +1)…𝑖𝑗+1 ],
𝐾𝑖𝑗 −𝑖2 …𝑖𝑛 = 𝐾𝑖𝑗 −𝑖2 = 𝑘𝑖1 − 𝑘𝑖2 , 𝑖𝑗 = 𝑖1 𝑖1 < 𝑖 𝑗 < 𝑖 𝑛
𝐾𝑖𝑗 −𝑖1 𝑖2 …𝑖𝑗−1 𝑖𝑗+1 …𝑖𝑛 = 𝐾𝑖𝑗 −𝑖𝑗−1 − 𝐾𝑖𝑗+1 −𝑖𝑗−1 , 𝑖1 < 𝑖𝑗 < 𝑖𝑛 𝐴𝑖𝑗 −𝑖1 𝑖2 … 𝑖𝑛−1 = −𝑎(𝑖𝑛−1 +1)…𝑖𝑛 , 𝑖𝑗 = 𝑖𝑛
∑𝑖𝑛−1 ∑𝑖
𝑚𝑏 (𝑚𝑏 ⋅𝑎𝑏 )
𝐾𝑖𝑗 −𝑖1 𝑖2 …𝑖𝑛−1 = 𝐾𝑖𝑗 − = 𝑘𝑖𝑛 + ∑𝑖 𝑛
𝑏=1
⋅ (𝑘𝑖𝑛 − 𝑘𝑖𝑛−1 ), 𝑖𝑗 = 𝑖𝑛 where 𝑎1⋅⋅⋅𝑖 = 𝑏=1
∑𝑖
𝑖𝑛−1 𝑏=𝑖𝑛−1 +1
𝑚𝑏 𝑏=1
𝑚𝑏
∑𝑖 𝑛 ∑𝑖1 −1
[ 𝑏=1 𝑚𝑏 ]⋅𝑎1…𝑖𝑛 −[ 𝑏=1 (𝑚𝑏 )]⋅𝑎1…(𝑖1 −1)
and 𝑎𝑖1 …𝑖𝑛 = ∑𝑖 𝑛
𝑏=𝑖1
𝑚𝑏

TA B L E 4 Sliding acceleration of blocks in different conditions under distributed acceleration–time series

ns Sliding acceleration
ns = 1 For block i, 𝑠̈𝑖 = −𝑎1…𝑖 − 𝑘𝑖 ⋅ 𝑔
ns > 1 For block i1 , i2 ,. . . , ij-1 , ij , ij+1 ,. . . , in-1 , in (i1 < i2 < . . . < in )
𝑠̈𝑖1 = −𝑎1…𝑖1 − [−𝑎(𝑖1 +1) …𝑖 ] + (𝐾𝑖2 −𝑖1 ⋅ 𝑔 − 𝑘𝑖1 ⋅ 𝑔)
2
𝑠̈𝑖𝑗 = −𝑎( 𝑖 +1)…𝑖 − [−𝑎(𝑖𝑗 +1) ] + (𝐾𝑖𝑗+1 −𝑖𝑗 ⋅ 𝑔 − 𝐾𝑖𝑗 −𝑖𝑗−1 ⋅ 𝑔)
𝑗−1 𝑗 …𝑖𝑗+1
𝑠̈𝑖𝑛 = −𝑎(𝑖𝑛−1 +1)…𝑖𝑛 − 𝐾𝑖𝑛 −𝑖𝑛−1 ⋅ 𝑔

𝑖2
[ 𝑖1
] [ 𝑖2
]
∑ ∑ ∑
Block# 𝑖2 ∶ − (𝑚𝑏 ⋅ 𝑎𝑏 ) + 𝑘𝑖1 ⋅ 𝑚𝑏 ⋅ 𝑔 = 𝑘𝑖2 ⋅ 𝑚𝑏 ⋅ 𝑔, −𝑎(𝑖1 +1)𝑖2 = 𝐾𝑖2 −𝑖1 ⋅ 𝑔, (18b)
𝑏=𝑖1 +1 𝑏=1 𝑏=1

where a(i1 +1)…i2 is the average acceleration–time history from block (i1 + 1) to block i2 , and is defined as
[∑ ] [∑ ]
∑𝑖2 ∑𝑖2 ∑𝑖1 𝑖2 𝑖1
𝑚 ⋅ 𝑎 − (𝑚 ) ⋅ 𝑎1…𝑖1
𝑏=𝑖1 +1 (𝑚𝑏 ⋅ 𝑎𝑏 ) 𝑏=1 (𝑚𝑏 ⋅ 𝑎𝑏 ) − 𝑏=1 (𝑚𝑏 ⋅ 𝑎𝑏 ) 𝑏=1 𝑏 1…𝑖2 𝑏=1 𝑏
𝑎(𝑖1 +1)…𝑖2 = ∑𝑖2 = ∑𝑖2 = ∑𝑖2 . (19)
𝑏=𝑖1 +1 𝑚𝑏 𝑏=𝑖1 +1 𝑚𝑏 𝑏=𝑖1 +1 𝑚𝑏

Similarly, the slip episodes of a generalized sliding block system with distributed acceleration–time series can be
resolved with the explicit descriptions of the initial condition of sliding and the sliding accelerations for different cases,
which are summarized in Tables 3 and 4.
In the application of the multiblock sliding model, the first step is to divide the slope into a series of critical failure
masses and to obtain the corresponding yield coefficient. The discretization can be performed at potentially critical slip
surfaces (e.g., at the interfaces between soil layers with different shear strengths and zones of weakness, or at geosynthetic
interfaces) to accurately describe the multiple slip failure or the dispersed shear movement that may occur during earth-
quakes. The distributed acceleration–time histories of sliding masses above each potential slip surface are then determined
from the seismic response analysis.

4 MODEL VERIFICATION

4.1 Verification of a multiblock model

To evaluate the multiblock seismic sliding model, a comparison with numerical results from the FLAC was performed,
with reference to a three-block system, as shown in Figure 9. The upper three layers represent the rigid blocks and the
lower layer is the ground base. The inclination of the plane is 11.4◦ . The rigid blocks and the base were modeled as elastic
material using the typical parameters of structural steel (the Young’s modulus is 206 GPa, the Poisson’s ratio is 0.25, and
the density is 7.8 g/cm3 ). Interface elements were applied between the three blocks and between the lower block and
the ground to model preexisting potentially sliding surfaces. The interfaces were characterized by Coulomb sliding with a
friction angle of 12.6◦ , 17.1◦ , and 19.9◦ for the upper, middle, and lower ones, respectively. A harmonic sinusoidal excitation
SONG et al. 11

F I G U R E 9 Finite difference model of an

inclined three-block system

F I G U R E 1 0 Acceleration of (A) the upper block, (B) the middle block, and (C) the lower block from the seismic sliding analysis
compared to the numerical results

with a frequency of 1.5 Hz and an amplitude of 0.3 g was applied as an input acceleration to the base of the model. A small
amount of stiffness-proportional Rayleigh damping (0.5%) was added to avoid numerical noise at high frequencies.
As the input motion was exerted horizontally, the horizontal yield coefficient was calculated based on Equation (3), and
the resulting horizontal yield coefficients of the upper, middle, and lower interfaces are 0.02, 0.1, and 0.15, respectively,
based on the friction angles and the inclination. Figure 10 shows the horizontal acceleration–time histories of the three
blocks from the multiblock seismic sliding analysis and the numerical results. Observe that for the upper block, its accel-
eration reaches the yield acceleration (0.02 g) when the relative sliding occurs, while for the middle and lower blocks,
their accelerations in the sliding condition (0.14 g for Block # 2 and 0.2 g for Block # 3) are larger than the corresponding
yield accelerations. These larger accelerations coincide with the conditional yield accelerations (K2-1 g and K3-2 g) computed
from Table 1. This observation demonstrates the influence of interactions among the sliding scenario of each block in the
multiblock system, and the developed analytical method can well capture the interactions and explain the mechanism.
Figure 11 presents the horizontal sliding displacement of the blocks predicted by the theoretical and numerical analyses.
12 SONG et al.

F I G U R E 1 1 Relative sliding displacement of (A) the upper block, (B) the middle block, and (C) the lower block from the seismic sliding
analysis compared to the numerical results

The results show a good agreement for all blocks. The computed sliding displacements from the traditional Newmark
analysis that treats the system as independent rigid blocks deviate from the values that consider the interactions among
the sliding blocks.

4.2 Verification of a simple slope model

Additional verification analyses were also performed for the numerical results of the simple slope model for the case
of Figure 2B in Section 2. The horizontal yield accelerations associated with the slip surfaces were first obtained as the
horizontal component of gravitational acceleration of the numerical model that results in a factor of safety of the slope
equal to 1 (0.05 g for the shallow slip surface and 0.06 g for the deep slip surface). Average acceleration–time histories
of sliding mass above the shallow slip surface and above the deep slip surface were then derived from the dynamic site
response analysis without consideration of soil plasticity. The displacement in the horizontal direction of the shallow and
the deep sliding masses relative to the deep slip surface are shown in Figure 12 for both rigid and decoupled multiblock
sliding analysis (i.e., uniform and distributed acceleration–time series, respectively), along with the displacement from
the traditional decoupled single-block analysis and the numerical results. Note that the displacement for the shallow
sliding mass here is thus the sum of the relative horizontal sliding displacement that occurs along the two slip surfaces.
This is different from the displacement presented in other figures of this paper, which is the relative horizontal sliding
displacement that occurs only along each of the slip surface. Because the displacements are different throughout the
sliding mass for the numerical model, the displacement range and the average displacement of the sliding mass are shown
in the figure.
In general, the decoupled multiblock model can better describe the development of deformation for the numerical
model than the rigid multiblock and the traditional decoupled single-block models. For the shallow sliding mass, the
displacements of all the sliding block methods deviate from the displacement range of numerical model, but the decou-
pled multiblock method provides results that are closer to the maximum displacement of the numerical model. The rigid
multiblock and the decoupled single-block analyses result in displacement estimates that are significantly unconserva-
tive and conservative. The unconservative results from the rigid multiblock are due to the neglection of the seismic site
SONG et al. 13

F I G U R E 1 2 Displacement in the horizontal direction relative to the deep slip surface from the rigid and decoupled multiblock sliding
analysis, and the decoupled single-block analysis compared to the numerical results (the range and average displacement): (A) shallow
sliding mass and (B) deep sliding mass

TA B L E 5 Properties of the four soil layers

Soil name t (m) γ (kN/m3 ) c (kPa) φ (◦ ) Vs (m/s)
Soil-1 10 17.5 20 18 180
Soil-2 15 18 30 21 240
Soil-3 15 19 42 25 320
Soil-4 20 20 50 27 460
Note: t is the soil-layer thickness.

amplification, and the conservative results from traditional decoupled single-block analysis are due to the neglection of
the interaction of sliding along both slip surfaces. For the deep sliding mass, the displacement from the decoupled multi-
block method is closer to the average displacement of the numerical model although all the three sliding block methods
provide displacements within the range of those predicted by the numerical model.
Therefore, although being conservative, the decoupled multiblock model generally provides acceptable results com-
pared with the numerical ones. The conservative estimation of the decoupled multiblock method may be explained by
three considerations. First, previous studies for the single sliding block analysis have demonstrated that the decoupled
approximation could be conservative (e.g., Refs. 7 and 9). Second, the failure mode from the numerical model is close to a
rotational mode instead of the translational mode. The change in geometry may result in an increasing yield acceleration
during sliding, and thus a smaller permanent displacement. Finally, the mechanisms of the displacement development
from the sliding block and the numerical analyses are not exactly the same.

5 SAMPLE APPLICATION AND DISCUSSION

5.1 Application of a four-layer soil deposit

To illustrate the application of the proposed approach to a realistic case, a hypothetical infinite-layered slope is considered.
The slope angle is 19◦ , and the slope is a four-layer deposit (referred to, from softest to strongest, as Soil-1, Soil-2, Soil-3,
and Soil-4). The soil properties of these layers are given in Table 5. These soil properties are selected using typical ranges of
fine-grained soils, with the intent of generating sites with heterogeneous soil layers. Different locations of Soil-1 (i.e., the
softest layer) are considered, resulting in three different sites (i.e., Site-1, Site-2, and Site-3; see Figure 13). The analysis in
this subsection will focus on Site-1, in which the shear strength of soil increases with depth (i.e., the soil layer is Soil-1, Soil-
2, Soil-3, and Soil-4 from the top to the bottom layer, respectively), and the other sites will be discussed later. Figure 13A
shows that the horizontal yield coefficient profile of this site was determined from pseudostatic analysis (Equation (3))
by considering infinite slopes with a varying sliding surface. Observe that there are discontinuities occurred in the yield
coefficient at the soil-layer interfaces. This may result in a susceptibility to sliding at these locations during the earthquake
14 SONG et al.

FIGURE 13 Yield coefficient profile for different soil deposits (different locations of the softest layer Soil-1)

TA B L E 6 Yield coefficient at different locations of slip surfaces for different layered soil deposits
Yield coefficient k
First-slip Fourth-slip
surface Second-slip Third-slip surface
(topmost) surface surface (bottommost)
Site-1 0.09 0.1 0.16 0.18
Site-2 0.13 0.03 0.16 0.18
Site-3 0.13 0.17 0.01 0.18

shaking. Therefore, the infinite slope was discretized into four potential sliding surfaces and the corresponding yield
coefficients for Site-1 are shown in Table 6.
The 1D equivalent-linear seismic site response analysis code Strata23 was employed to compute the average acceleration
time histories above the potential slip surfaces (referred to as k-acc here). The k-acc can be calculated from the ratio of
the shear stress to the total vertical stress at the slip surface in the 1D case. The nonlinear soil properties were modeled
using the curves of Darendeli and Stokoe24 with a plasticity index of 20 and appropriate values of confining pressure
for soils at different depths. The input earthquake ground motion was also the El Centro Array #6 recorded in the 1979
Imperial Valley-06 earthquake. The peak values of the average acceleration–time and velocity–time histories above the
potential slip surfaces at different depths of Site-1 are presented in Table 7. The average velocity–time history (referred
to as k-vel) was generated through numerical integration of the average acceleration–time history. Observe that the peak
values of average acceleration–time histories (k-accmax ) are smaller than the input PGA (PGA = 0.45 g). However, the peak
values of average velocity–time histories (k-velmax ) are larger than the value of input PGV (114 cm/s) due to the increase
in long-period contents. These average acceleration–time histories of sliding masses above different slip surfaces provide
the appropriate seismic loading for the multiblock analysis.
Based on the yield coefficients associated with the slip surfaces and the average acceleration responses of the sliding
mass system, the seismic displacement can be determined through the multiblock sliding approach, and are shown in
Figure 14 for Site-1, where both the results of rigid (i.e., uniform dynamic responses) and flexible (i.e., distributed dynamic
responses) sliding systems are presented. Observe that the sliding is triggered along the two shallower slip surfaces both in
the rigid and flexible models, but the displacements are larger for the flexible sliding masses due to the dynamic response of

T A B L E 7 Peak values of the average acceleration–time and velocity–time histories above different locations of slip surfaces for different
layered soil deposits
k-accmax (g) k-velmax (cm/s)
First Second Third Fourth First Second Third Fourth
layer layer layer layer layer layer layer layer
Site-1 0.40 0.32 0.26 0.26 143 135 126 119
Site-2 0.29 0.28 0.27 0.23 164 153 124 112
Site-3 0.28 0.27 0.24 0.27 166 157 143 113
SONG et al. 15

FIGURE 14 Relative sliding displacement at different slip surfaces of (A) rigid and (B) flexible sliding systems

F I G U R E 1 5 Relative sliding displacement at different slip surfaces of flexible sliding system of Site-1 subject to ground motions with
different amplitudes: (A) scaling factor = 0.5 and (B) scaling factor = 2.0

the site. The type of progressive sliding on multiple potential slip surfaces is a unique failure mode under the time-varying
earthquake loading, which may not be expected in static slope stability problems where the failure is constrained on the
single-slip surface with the minimum factor of safety. The seismic response analysis in the application of the proposed
multiblock sliding model was based on simplified 1D analysis. Extension to consider topographic effects of slopes may be
implemented using the same process presented herein.

5.2 Various amplitudes of earthquake ground motion

To further clarify the development of multiple failure surfaces associated with coseismic landslides, the ground motion
was linearly scaled to different amplitudes, and the sliding displacements of the block system in Site-1 including the
dynamic site response are shown in Figure 15. As can be seen from the figure, the sliding is only triggered for the block
with the minimum yield acceleration (i.e., Block # 1) for the ground motion with a scaling factor of 0.5. As the intensity of
the input ground motion increases, subsequent sliding along the deeper slip surfaces occurs (i.e., Block # 1 and Block # 2
for nonscaled motion, and all the four blocks for a scaling factor of 2.0). This provides a means to refine the assessment of
earthquake-induced landslide hazards by including both the occurrence and the consequence of different failure modes.
In the traditional displacement-based seismic hazard analysis of slopes, higher ground motion intensities only result in
higher probabilities of slope failure at the predefined sliding surface. A multilevel seismic landslide hazard assessment
combining a number of failure modes and also their potential effects may be performed based on the developed approach.

5.3 Alternative soil layering

The influence of soil deposits is examined by considering different locations of the softest layer Soil-1. The yield coefficient
profile for different soil deposits is shown in Figure 13, and the yield coefficients at different locations of discretized slip
surfaces are shown in Table 6. The change of the location of the softest layer Soil-1 results in the shift of the slip surface with
16 SONG et al.

FIGURE 16 Relative sliding displacement of flexible blocks of (A) Site-2 and (B) Site-3

F I G U R E 1 7 Equivalent seismic loading and yield coefficients conditioned on the sliding of Block # 2 for Site-2: (A) Block # 1, (B) Block
#
3, and (C) Block # 4

the minimum yield coefficient. Note that the soil deposit with Soil-1 at the bottom soil layer is not present herein because
the minimum yield coefficient would be negative (i.e., the slope is initially unstable before earthquakes). Seismic site
response analysis was performed for Site-2 and Site-3, and the peak values of the average acceleration–time and velocity–
time histories above different locations of slip surfaces are presented in Table 7. Observe that the presence of a relative soft
soil layer at deeper locations of the strata would further decrease the high-frequency contents (smaller values of k-accmax )
and increase the low-frequency contents (larger values of k-velmax ) of the response of sliding masses above.
The computed sliding displacements of the flexible blocks at different sites are shown in Figure 16. It can be seen that
both the sliding of Block # 1 (11.8 cm) and Block # 2 (480.5 cm) is triggered at Site-2, and only the sliding associated with a
very large displacement of Block # 3 (1187.8 cm) is triggered at Site-3. Note that the yield coefficient of the topmost sliding
surface (0.13) is much larger than the minimum yield coefficient (0.03) for Site-2, but a sliding displacement of 11.8 cm
occurs at the topmost sliding surface, which would not be this case for rigid blocks. This can be explained from Figure 17,
in which the equivalent seismic loading and conditional yield accelerations on the sliding of Block # 2 are illustrated. After
the sliding of Block # 2 is initiated, there are some moments of time when the equivalent seismic acceleration exceeds the
conditional yield acceleration of Block # 1 due to the dynamic site response. However, this is not the case for the two deeper
slip surfaces of Site-2 and the slip surfaces of Site-3 (the results are similar and are not presented), and hence the sliding
does not occur.
SONG et al. 17

6 CONCLUSIONS

A Newmark-type multiblock model was developed to describe the complex sliding scenarios of slopes exhibiting multiple
failure surfaces, which may occur in heterogeneous soils during time-varying earthquake shaking. The model discretizes
the infinite slope into a generic number of sliding blocks, and the slip episodes of the block system are rigorously resolved
for uniform and distributed acceleration–time histories throughout the sliding masses under earthquake loading. The dis-
cretization can be performed reproducing the site-specific conditions (e.g., at the interfaces of the different soil layers or at
the preexisting potentially sliding surfaces) to accurately describe the multiple slip failures that may occur during earth-
quakes. The distributed acceleration–time histories of sliding masses above each potential slip surface can be determined
from seismic response analyses.
The conditional yield accelerations and the equivalent seismic accelerations were defined to characterize the complex
interactions of the block system, and to determine the initial condition triggering sliding and the sliding accelerations
for different cases (i.e., cases in which different slip surfaces initiate sliding). The development of multiple failure
surfaces was investigated for sites with different soil layering. The ground motion was applied uniformly throughout the
profile and was allowed to vary with depth, respectively. It is found that for the uniform acceleration–time series, the
initial sliding would occur along the failure surface with the minimum yield acceleration. In this case, the blocks above
this failure surface are locked and the blocks below the failure surface would have larger yield accelerations due to the
isolation effects that result from sliding. This would be different when considering the seismic site response due to varying
acceleration–time series throughout the sliding masses. The solution of the multiblock model allows the possibility of
combining the different failure modes and their potential effects in the assessment of the earthquake-induced landslide
hazard.
Although the proposed approach introduces the compliant response of slopes by using a decoupled manner, it is a
significant enhancement of the rigid-block system. Future work may incorporate coupled dynamic response and sliding
analyses into the multiblock seismic sliding model. This would result in the most rigorous procedure to fully describe
the interactions of “stick-slip” and “slip-slip” process of sliding mass systems, but not without a significant increase in
complexity.

AC K N OW L E D G M E N T S
This research has been supported by the National Natural Science Foundation of China (Grant No. 41630638 and 51879091),
the Fundamental Research Funds for the Central Universities in China (Grant No. B210202046), and the National Key
Research and Development Project of China (Grant No. 2018YFC1508604).

D A T A AVA I L A B I L I T Y S T A T E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.

ORCID
Tugen Feng https://orcid.org/0000-0002-6638-5818

REFERENCES
1. Cui SH, Wang GH, Pei XJ, Huang RQ, Kamai T. On the initiation and movement mechanisms of a catastrophic landslide triggered by the
2008 Wenchuan (M-s 8.0) earthquake in the epicenter area. Landslides. 2017;14(3):805-819.
2. Bray JD, Travasarou T. Simplified procedure for estimating earthquake-induced deviatoric slope displacements. J Geotech Geoenviron Eng
ASCE. 2007;133(4):381-392.
3. Jibson RW. Methods for assessing the stability of slopes during earthquakes—A retrospective. Eng Geol. 2011;122:43-50.
4. Ambraseys NN, Menu JM. Earthquake-induced ground displacements. Earthq Eng Struct Dyn. 1988;16(7):985-1006.
5. Du W, Wang G. Fully probabilistic seismic displacement analysis of spatially distributed slopes using spatially correlated vector intensity
measures. Earthq Eng Struct Dyn. 2014;43(5):661-679.
6. Newmark NM. Effects of earthquakes on dams and embankments. Géotechnique. 1965;15(2):139-160.
7. Lin JS, Whitman RV. Decoupling approximation to the evaluation of earthquake-induced plastic slip in dams. Earthq Eng Struct Dyn.
1983;11(5):667-678.
8. Kramer SL, Smith MW. Modified Newmark model for seismic displacements of compliant slopes. J Geotech Geoenviron Eng.
1997;123(7):635-644.
9. Rathje EM, Bray JD. Nonlinear coupled seismic sliding analysis of earth structures. J Geotech Geoenviron Eng. 2000;126(11):1002-1014.
18 SONG et al.

10. Tropeano G, Chiaradonna A, D’Onofrio A, Silvestri F. An innovative computer code for 1D seismic response analysis including shear
strength of soils. Géotechnique. 2016;66(2):95-105.
11. Makdisi FI, Seed HB. Simplified procedure for estimation dam and embankment earthquake induced deformations. J Geotech Eng Div
ASCE. 1978;104(GT7):849-867.
12. Chopra AK. Earthquake response of earth dams. J Soil Mech Found Div ASCE. 1967;93(2):65-81.
13. Rathje EM, Antonakos G. A unified model for predicting earthquake-induced sliding displacements of rigid and flexible slopes. Eng Geol.
2011;122(1-2):51-60.
14. Song J, Gao YF, Feng TG, Xu GZ. Effect of site condition below slip surface on prediction of equivalent seismic loading parameters and
sliding displacement. Eng Geol. 2018;242:169-183.
15. Tsai CC, Lin CH. Prediction of earthquake-induced slope displacements considering 2D topographic amplification and flexible sliding
mass. Soil Dyn Earthq Eng. 2018;113:25-34.
16. Stewart JP, Bray JD, McMahon DJ, Smith PM, Kropp AL. Seismic performance of hillside fills. J Geotech Geoenviron Eng. 2001;127(11):905-
919.
17. Wartman J, Seed RB, Bray JD. Shaking table modeling of seismically induced deformations in slopes. J Geotech Geoenviron Eng.
2005;131(5):610-622.
18. Song J, Fan QQ, Feng TG, Chen ZQ, Chen J, Gao YF. A multi-block sliding approach to calculate the permanent seismic displacement of
slopes. Eng Geol. 2019;255:48-58.
19. Itasca Consulting Group. FLAC: Fast Lagrangian analysis of continua, 2005.
20. Sun JI, Golesorkhi R, Seed HB. Dynamic Moduli and Damping Ratios for Cohesive Soils. Report No. UCB/EERC-88/15, Earthquake Engi-
neering Research Center, University of California, Berkeley; 1988.
21. Chiu CC, Weng MC. DEM simulation of planar sliding using a particulate interface model considering velocity-dependent friction. Comput
Geotech. 2019;112:51-59.
22. Stamatopoulos CA. Sliding system predicting large permanent co-seismic movements of slopes. Earthq Eng Struct Dyn. 1996;25(10):1075-
1093.
23. Kottke EM, Rathje EM. Technical Manual for Strata. PEER Report 2008/10, Pacific Earthquake Engineering Research Center, University
of California at Berkeley; 2008:84.
24. Darendeli MB. Development of a New Family of Normalized Modulus Reduction and Material Damping Curves. Geotechnical Engineering
Report GD01-1, University of Texas, Austin, TX; 2001.

How to cite this article: Song J, Rodriguez-Marek A, Feng T, Ji J. A generalized seismic sliding model of slopes
with multiple slip surfaces. Earthquake Engng Struct Dyn. 2021;1–18. https://doi.org/10.1002/eqe.3462