Numerical Study on Dynamic Response of a Horizontal Layered-Structure Rock Slope under a Normally Incident Sv Wave

Several post-earthquake investigations have indicated that the slope structure plays a leading role in the stability of rock slopes under dynamic loads. In this paper, the dynamic response of a horizontal layered-structure rock slope under harmonic Sv wave is studied by making use of the Fast Lagrangian Analysis of Continua method (FLAC). The suitability of FLAC for studying wave transmission across rock joints is validated through comparison with analytical solutions. After parametric studies on Sv wave transmission across the horizontal layered-structure rock slope, it is found that the acceleration amplification coefficient η, which is defined as the ratio of the acceleration at the monitoring point to the value at the toe, wavily increases with an increase of the height along the slope surface. Meanwhile, the fluctuation weakens with normalized joint stiffness K increasing and enhances with normalized joint spacing ξ increasing. The acceleration amplification coefficient of the slope crest ηcrest does not monotonously increase with the increase of ξ, but decreases with the increase of K. Additionally, ηcrest is more sensitive to ξ compared to K. From the contour figures, it can also be found that the contour figures of η take on rhythm, and the effects of ξ on the acceleration amplification coefficient are more obvious compared to the effects on K.


Introduction
It is widely recognized that the dynamic response of rock slopes is one of the key issues in geotechnical engineering and earthquake engineering, and there is indisputable evidence of its significance even from the late 1960s [1,2]. On the other hand, a rock mass usually consists of multiple joints, known as joint sets, which govern the mechanical behavior of the rock mass [3]. The stability of rock slopes is often significantly influenced by the structural geology of the rock, especially discontinuities such as bedding planes, joints and faults [4]. Studies [5][6][7] have confirmed that the rock mass structure plays a leading role in the stability of rock slopes under dynamic loads. Therefore, research on wave propagation regularity in a layered-structure rock slope is important for assessing the stability and damage of rock slopes under dynamic loads.
Hence, many researchers have studied the dynamic response of rock slopes under stress and tried to simulate the complexity using analytical, numerical, model, and in situ test methods, and valuable conclusions about topography effects and dynamic response regularities have been drawn [8][9][10][11][12]. However, analytical solutions encounter serious limitations when the geometry is complicated and a great number of influencing parameters are to be incorporated [13]. The results obtained by model test and in situ test methods are limited and cannot be generalized. Due to complicated topography and slope structure, results from in situ test are difficult to obtain. Furthermore, the conclusion from model tests cannot reflect the natural slope. Meanwhile, numerical simulation is an economical and convenient alternative. Few works have been conducted to systemically study plane wave transmission through a layered-structure rock slope. As a result, numerical simulation of the dynamic response of a horizontal layered-structure rock slope is useful.
When a wave propagates through fractured rock masses, its attenuation and slowness are mainly induced by rock joints [14,15]. Analytical solutions and experimental data have proven that the incidence angle of waves, normalized fracture stiffness, number of fractures, and dimensionless fracture spacing have obvious effects on wave attenuation across parallel fractures [16,17].
Many characteristic parameters of slope and incident waves have effects on wave attenuation through a horizontal layered-structure rock slope, such as height, slope angle, acceleration amplitude, frequency, and so on. Several papers [18,19] have studied the effects of dynamic actions on the shear strength of rock mass discontinuities. In order to focus on studies of the parameters in which we are most interested, a homogeneous and elastic rock slope with normally incident Sv wave applied under the bottom of the slope is modeled, making full use of the Fast Lagrangian Analysis of Continua (FLAC) method (with the code of FLAC3D V5.0 issued by ITASCA Consulting Group, Inc., Minneapolis, MN, USA). FLAC is a three-dimensional explicit finite-difference program for engineering mechanics computation. This method follows the continuum hypothesis, which adopts the difference format. With continual change of the configuration, the coordinates update according to the solution of the time step integral. Based on the dynamic equation and dynamic method, it is able to simulate the dynamic problem well. One advantage of FLAC is that the FISH language is embedded; therefore, it is convenient to define the constitutive relation and to arrange the monitoring points. Before the simulation of the dynamic response of the slope, the verification of FLAC modeling on wave transmission across rock joints is executed. Through comparison with analytical solutions and FLAC results, the effects of joint mechanical and spatial properties on wave transmission in FLAC modeling are validated.
In this paper, as a simplified slope model and normally incident Sv wave are given, the incidence angle of waves and number of fractures are fixed, so the effects of normalized fracture stiffness and dimensionless fracture spacing to the dynamic response regularities are the main study factors. Subsequently, dynamic response and parametric studies on Sv wave propagation through horizontal layered-structure rock slope are performed. At the same time, dynamic response regularities are presented in detail.

Validation of FLAC Modeling on Wave Transmission across Rock Joints
In the slope model, the normally incident wave is applied under the bottom of the slope. Before the transmitted wave arrives at the slope surface, just wave reflection and wave transmission exist. When the wave transmits to the slope surface, wave transformation as well as wave reflection occur. Subsequently, obliquely incident wave transmission across a joint set also occurs. Therefore, wave transmission and attenuation are very complicated due to wave transformation, wave reflection, and wave superposition.
In order to verify the capability of FLAC on wave transmission across rock joints, the modeling results are compared with analytical solutions with the propagator matrix method [16], adopting the displacement discontinuity model. In this section, FLAC results of normally incident P and Sv waves transmission across a single joint are compared with analytical solutions, and the effects of joint normal, shear stiffness, and incident wave types (P and Sv waves) on wave transmission in FLAC modeling are verified. Effects of joint spatial configuration and multiple wave reflections among joints on wave transmission are subsequently verified through comparison with theoretical solutions for normally incident P wave transmission across a joint set. Finally, obliquely incident P wave transmission across a single joint is studied to verify the dependence of wave transmission on the incident angle and the two-dimensional effect [15].

Normally Incident P Wave Transmission through a Single Joint
The FLAC model for P wave transmission through a single joint is illustrated in Figure 1. The origin of the X-Z coordinates is located at the left bottom of the model. The model size is 500 m × 500 m, and the joint is horizontally located at Z = 250 m. Brick elements is adopted as the mesh type, and all the mesh size is 1 m. Non-reflection viscous boundaries are placed at the bottom and upper boundaries to avoid wave reflections from the artificial boundaries. The left and right boundaries are fixed in X-direction. And a half-cycle sinusoidal incident stress pulse with amplitude 1.662276 × 10 7 Pa (to make sure that the velocity amplitude is 1.0 m/s) and frequency 50 Hz is normally applied at the bottom boundary. Appl. Sci. 2017, 7, 716 3 of 16

Normally Incident P Wave Transmission through a Single Joint
The FLAC model for P wave transmission through a single joint is illustrated in Figure 1. The origin of the X-Z coordinates is located at the left bottom of the model. The model size is 500 m × 500 m, and the joint is horizontally located at Z = 250 m. Brick elements is adopted as the mesh type, and all the mesh size is 1 m. Non-reflection viscous boundaries are placed at the bottom and upper boundaries to avoid wave reflections from the artificial boundaries. The left and right boundaries are fixed in X-direction. And a half-cycle sinusoidal incident stress pulse with amplitude 1.662276 × 10 7 Pa (to make sure that the velocity amplitude is 1.0 m/s) and frequency 50 Hz is normally applied at the bottom boundary. The rock mass is assumed to be purely elastic, and material damping is not included in this paper for two reasons [15]. First, when a wave propagates across jointed rock mass, the wave attenuation is mainly due to joints. Secondly, this study is focused on effects of rock joints on wave attenuation and thus the attenuation from material is ignored. The joint model adopts the Coulomb slip model. In order to ensure that the joint is linear and elastic, high values are assigned to the cohesion and tension of the joint. In the numerical simulation, cohesion is 10 GPa, the friction angle is 30°, and the tension strength is 10 5 GPa. The physico-mechanical properties of a rock matrix are listed in Table 1 [20].  The rock mass is assumed to be purely elastic, and material damping is not included in this paper for two reasons [15]. First, when a wave propagates across jointed rock mass, the wave attenuation is mainly due to joints. Secondly, this study is focused on effects of rock joints on wave attenuation and thus the attenuation from material is ignored. The joint model adopts the Coulomb slip model. In order to ensure that the joint is linear and elastic, high values are assigned to the cohesion and tension of the joint. In the numerical simulation, cohesion is 10 GPa, the friction angle is 30 • , and the tension strength is 10 5 GPa. The physico-mechanical properties of a rock matrix are listed in Table 1 [20]. In this section, only the transmission coefficients are studied. The FLAC results and analytical solutions are compared in Figure 2. In this model, the incident and transmitted measure point are located at (X, Z) = (250, 50) and (250, 450) of Figure 1, respectively. To focus on the study of the relation between the transmission coefficient and joints stiffness, the normal joint stiffness k n is equal to the shear joint stiffness k s ; for convenience, both of them are expressed as k. The X-coordinate is normalized joint stiffness K, which is defined as: where z is given by: where k is the normal or shear joint stiffness, ω is the angular frequency of the incident wave that is the product of wave frequency f and 2π, z is the wave impedance, ρ is the rock density, and C refers to C P or C S , which is the velocity of the P or Sv wave. The vertical coordinate is the magnitude of transmission coefficient across a single joint for P wave incidence T 1P , which is defined as [18]: where v IP is the particle velocity of the incident P wave and v TP is the particle velocity of the transmitted P wave. From Figure 2, it can be found that T 1P from the FLAC results agrees well with the analytical solutions. Therefore, FLAC is applicable to study P wave transmission. In this section, only the transmission coefficients are studied. The FLAC results and analytical solutions are compared in Figure 2. In this model, the incident and transmitted measure point are located at (X, Z) = (250, 50) and (250, 450) of Figure 1, respectively. To focus on the study of the relation between the transmission coefficient and joints stiffness, the normal joint stiffness kn is equal to the shear joint stiffness ks; for convenience, both of them are expressed as k. The X-coordinate is normalized joint stiffness K, which is defined as: where z is given by: where k is the normal or shear joint stiffness, ω is the angular frequency of the incident wave that is the product of wave frequency f and 2π, z is the wave impedance, ρ is the rock density, and C refers to CP or CS, which is the velocity of the P or Sv wave. The vertical coordinate is the magnitude of transmission coefficient across a single joint for P wave incidence T1P, which is defined as [18]: where vIP is the particle velocity of the incident P wave and vTP is the particle velocity of the transmitted P wave. From Figure 2, it can be found that T1P from the FLAC results agrees well with the analytical solutions. Therefore, FLAC is applicable to study P wave transmission.

Normally Incident SV Wave Transmission through a Single Joint
In order to verify the capability of FLAC in modeling Sv wave transmission, in this section, normally incident Sv wave transmission across a single joint is studied. The model size is also 500 m × 500 m, and the joint is horizontally located at Z = 250 m. Non-reflection viscous boundaries are placed at the upper and lower boundaries to avoid wave reflections from the artificial boundaries. The left and right boundaries are fixed in the Z direction. A half-cycle sinusoidal incident stress pulse with amplitude 1.01772 × 10 7 Pa (to make sure that the velocity amplitude is 1.0 m/s) and frequency 50 Hz is normally applied at the bottom boundary. The physico-mechanical properties of the rock matrix are the same as those in Table 1. In this model, the incident and transmitted measure point are located at (X, Z) = (250, 50) and (250, 450) of Figure 1, respectively, as well as in Section 2.1.

Normally Incident S V Wave Transmission through a Single Joint
In order to verify the capability of FLAC in modeling Sv wave transmission, in this section, normally incident Sv wave transmission across a single joint is studied. The model size is also 500 m × 500 m, and the joint is horizontally located at Z = 250 m. Non-reflection viscous boundaries are placed at the upper and lower boundaries to avoid wave reflections from the artificial boundaries. The left and right boundaries are fixed in the Z direction. A half-cycle sinusoidal incident stress pulse with amplitude 1.01772 × 10 7 Pa (to make sure that the velocity amplitude is 1.0 m/s) and frequency 50 Hz is normally applied at the bottom boundary. The physico-mechanical properties of the rock matrix are the same as those in Table 1. In this model, the incident and transmitted measure point are located at (X, Z) = (250, 50) and (250, 450) of Figure 1, respectively, as well as in Section 2.1. The FLAC results and analytical solutions are compared in Figure 2. The horizontal coordinate is normalized joint stiffness K. The vertical coordinate is the magnitude of transmission coefficient across a single joint for Sv wave incidence T 1S , which is defined as: where v IS is the particle velocity of incident P wave, v TS is the particle velocity of transmitted P wave. From Figure 3, it can be found that T 1S from the FLAC results agrees well with the analytical solutions. Therefore, FLAC is applicable to study Sv wave transmission.
Appl. Sci. 2017, 7, 716 5 of 16 The FLAC results and analytical solutions are compared in Figure 2. The horizontal coordinate is normalized joint stiffness K. The vertical coordinate is the magnitude of transmission coefficient across a single joint for Sv wave incidence T1S, which is defined as: where vIS is the particle velocity of incident P wave, vTS is the particle velocity of transmitted P wave. From Figure 3, it can be found that T1S from the FLAC results agrees well with the analytical solutions. Therefore, FLAC is applicable to study Sv wave transmission.

Normally Incident P Wave Transmission through a Joint Set
When a wave propagates through a rock mass with multiple parallel joints, multiple reflections occur between joints, which makes the problem more complex [21]. The joint set involved is assumed to have the same mechanical parameters and physico-mechanical properties of the rock matrix as those in Table 1. The FLAC model is also the same as in Figure 1 except that multiple joints exist with different spacing. In this model, the incident and transmitted measure point are located at (X, Z) = (250, 50) and (250, 456), respectively.
The FLAC results and analytical solutions are compared in Figure 4. The horizontal coordinate is the dimensionless fracture spacing ξ, which is defined as the ratio of joint spacing d to incident wavelength λ, where K = 1. The vertical coordinate is the magnitude of transmission coefficient across a joint set for P wave incidence TN, which is defined as: where vIP is the particle velocity of incident P wave and vTP is the particle velocity of transmitted P wave through a joint set. From Figure 4, it can be found that TN from the FLAC results agrees well with the analytical solutions. Therefore, FLAC is applicable to study wave propagating through a rock mass with multiple parallel joints.

Normally Incident P Wave Transmission through a Joint Set
When a wave propagates through a rock mass with multiple parallel joints, multiple reflections occur between joints, which makes the problem more complex [21]. The joint set involved is assumed to have the same mechanical parameters and physico-mechanical properties of the rock matrix as those in Table 1. The FLAC model is also the same as in Figure 1 except that multiple joints exist with different spacing. In this model, the incident and transmitted measure point are located at (X, Z) = (250, 50) and (250, 456), respectively.
The FLAC results and analytical solutions are compared in Figure 4. The horizontal coordinate is the dimensionless fracture spacing ξ, which is defined as the ratio of joint spacing d to incident wavelength λ, where K = 1. The vertical coordinate is the magnitude of transmission coefficient across a joint set for P wave incidence T N , which is defined as: where v IP is the particle velocity of incident P wave and v TP is the particle velocity of transmitted P wave through a joint set. From Figure 4, it can be found that T N from the FLAC results agrees well with the analytical solutions. Therefore, FLAC is applicable to study wave propagating through a rock mass with multiple parallel joints.

Obliquely Incident P Wave Transmission through a Single Joint
In this slope model, the inclined angle is 30°. For an oblique slope surface, obliquely incident wave transmission across joints occurs after being reflected from the oblique surface. In order to verify the dependence of wave transmission on the incident angle and the two-dimensional effect in FLAC, an obliquely incident P wave through a horizontal joint is given, and the wave transformation, wave reflection as well as wave transmission are shown in Figure 5. As the wave transmission is two-dimensional, in order to eliminate the effects of wave superposition of the reflected wave from the upper and side boundaries on the transmission coefficient, the point (X, Z) = (250 m, 251 m) used to monitor the transmitted wave is near the joint. Through this method, the transmitted wave across the joint and the reflected wave from the upper and side boundaries can easily be separated because of the time difference in the arrival of the P wave and the Sv wave.

Obliquely Incident P Wave Transmission through a Single Joint
In this slope model, the inclined angle is 30 • . For an oblique slope surface, obliquely incident wave transmission across joints occurs after being reflected from the oblique surface. In order to verify the dependence of wave transmission on the incident angle and the two-dimensional effect in FLAC, an obliquely incident P wave through a horizontal joint is given, and the wave transformation, wave reflection as well as wave transmission are shown in Figure 5. As the wave transmission is two-dimensional, in order to eliminate the effects of wave superposition of the reflected wave from the upper and side boundaries on the transmission coefficient, the point (X, Z) = (250 m, 251 m) used to monitor the transmitted wave is near the joint. Through this method, the transmitted wave across the joint and the reflected wave from the upper and side boundaries can easily be separated because of the time difference in the arrival of the P wave and the Sv wave.

Obliquely Incident P Wave Transmission through a Single Joint
In this slope model, the inclined angle is 30°. For an oblique slope surface, obliquely incident wave transmission across joints occurs after being reflected from the oblique surface. In order to verify the dependence of wave transmission on the incident angle and the two-dimensional effect in FLAC, an obliquely incident P wave through a horizontal joint is given, and the wave transformation, wave reflection as well as wave transmission are shown in Figure 5. As the wave transmission is two-dimensional, in order to eliminate the effects of wave superposition of the reflected wave from the upper and side boundaries on the transmission coefficient, the point (X, Z) = (250 m, 251 m) used to monitor the transmitted wave is near the joint. Through this method, the transmitted wave across the joint and the reflected wave from the upper and side boundaries can easily be separated because of the time difference in the arrival of the P wave and the Sv wave.  The FLAC results and analytical solutions are compared in Figure 6. The horizontal coordinate is the incident angle θ. The vertical coordinate is the magnitude of transmission coefficient across a single joint for P wave oblique incidence T 1P , which is defined as Equation (3). From Figure 3, it can be found that T 1P from the FLAC results agrees well with the analytical solutions. Therefore, FLAC is applicable to study the dependence of wave transmission on the incident angle and the two-dimensional effect.
Appl. Sci. 2017, 7, 716 7 of 16 The FLAC results and analytical solutions are compared in Figure 6. The horizontal coordinate is the incident angle θ. The vertical coordinate is the magnitude of transmission coefficient across a single joint for P wave oblique incidence T1P, which is defined as Equation (3). From Figure 3, it can be found that T1P from the FLAC results agrees well with the analytical solutions. Therefore, FLAC is applicable to study the dependence of wave transmission on the incident angle and the two-dimensional effect.

Dynamic Response Analysis Method
As many factors could influence the dynamic response of the slope, a homogeneous and elastic rock slope with horizontal layered-structure, slope angle of 30°, and height of 30 m is specified to focus on the study of the parameters in which we are most interested (normalized fracture stiffness and dimensionless fracture spacing). As Figure 7 shows, the number of the joint in the slope model (N) is 5; the equal space of the joint (d) is 10 m. The joint model adopts the Coulomb slip model, which is commonly used in FLAC. High values are assigned to the cohesion and tension of the joint to make sure that the joint is linear and elastic. As in Section 2.1, cohesion is 10 GPa, the friction angle is 30°, and the tension strength is 10 5 GPa.

Dynamic Response Analysis Method
As many factors could influence the dynamic response of the slope, a homogeneous and elastic rock slope with horizontal layered-structure, slope angle of 30 • , and height of 30 m is specified to focus on the study of the parameters in which we are most interested (normalized fracture stiffness and dimensionless fracture spacing). As Figure 7 shows, the number of the joint in the slope model (N) is 5; the equal space of the joint (d) is 10 m. The joint model adopts the Coulomb slip model, which is commonly used in FLAC. High values are assigned to the cohesion and tension of the joint to make sure that the joint is linear and elastic. As in Section 2.1, cohesion is 10 GPa, the friction angle is 30 • , and the tension strength is 10 5 GPa. The FLAC results and analytical solutions are compared in Figure 6. The horizontal coordinate is the incident angle θ. The vertical coordinate is the magnitude of transmission coefficient across a single joint for P wave oblique incidence T1P, which is defined as Equation (3). From Figure 3, it can be found that T1P from the FLAC results agrees well with the analytical solutions. Therefore, FLAC is applicable to study the dependence of wave transmission on the incident angle and the two-dimensional effect.

Dynamic Response Analysis Method
As many factors could influence the dynamic response of the slope, a homogeneous and elastic rock slope with horizontal layered-structure, slope angle of 30°, and height of 30 m is specified to focus on the study of the parameters in which we are most interested (normalized fracture stiffness and dimensionless fracture spacing). As Figure 7 shows, the number of the joint in the slope model (N) is 5; the equal space of the joint (d) is 10 m. The joint model adopts the Coulomb slip model, which is commonly used in FLAC. High values are assigned to the cohesion and tension of the joint to make sure that the joint is linear and elastic. As in Section 2.1, cohesion is 10 GPa, the friction angle is 30°, and the tension strength is 10 5 GPa.  Through using the FISH language embedded in FLAC, the monitoring points are laid out at equal intervals in the upper slope in which we are most interested. The duration of Sv wave applied under the bottom of the slope was 1.0 s, allowing the incident wave to transmit between the bottom and top boundaries of the slope several times.
For a slope with multiple horizontally layered structures, analytical solution of wave transmission through the horizontal layered-structure rock slope is impossible to obtain and experiments as well as in situ tests are difficult to conduct. Therefore, numerical simulation is the most appropriate approach to study the complexity of this problem.
A FLAC model with five horizontally equal interval joints (d = 10 m), slope angle of 30 • , and height of 30 m is shown in Figure 7. The viscous boundaries are placed at the left side, right side, and bottom of the model. The incident harmonic Sv wave is applied at the bottom boundary of the slope model. The monitoring points are evenly distributed between the joints in the slope. The origin of the X-Z coordinates is located at the toe of the slope.
It is known that the horizontal acceleration is greatly amplified at the crest of slopes [12]. In this section, the horizontal acceleration (i.e., X-acceleration) of monitoring points in the slope shown in Figure 7 is the concerned factor. For convenience, the horizontal acceleration is described as the acceleration, and the acceleration amplification coefficient of monitoring point is defined as η, which can be expressed as: where A 0 is the acceleration amplitude of the reference point near the toe of the slope and A n is the acceleration amplitude of any monitoring point in the slope shown in Figure 7.
In order to study the dynamic response regularities of the slope, two parameters are considered. They are the normalized fracture stiffness (K) and the dimensionless fracture spacing (ξ), respectively. Through the monitoring points placed in the slope profile and Equation (6), the acceleration amplification coefficients of monitoring points in the whole slope can be obtained. It should be noted that the mechanical properties, spatial properties, and incident wave properties of the actual slope are very complicated. In this article, with the rock slope model specified, many factors are fixed. Therefore, the main study can focus on the most important parameters.

Computational Model
In static analyses, fixed or elastic boundaries can be realistically placed at some distance from the region of interest. In dynamic problems, however, such boundary conditions cause the reflection of outward propagating waves back into the model and do not allow the necessary energy radiation. The use of a larger model can minimize the problem, since material damping will absorb most of the energy in the waves reflected from distant boundaries. However, this solution leads to a large computational burden [22]. At the same time, the slope simulated here is not a structure, but a part of the geological body, and how large a model is needed to accurately model the behavior of the slope is not clear yet [23]. Thus, the alternative is to use quiet (or absorbing) boundaries. The viscous boundaries developed by Kuhlemeyer and Lysmer [24] are used in our simulations; they are placed at the left side, right side, and bottom of the model in the normal and shear directions (Figure 7). The dashpots provide viscous normal and shear tractions [22]: where v n and v s are the normal and shear components of the velocity at the boundary, ρ is the mass density, and C P and C S are the P and Sv wave velocities. C P is given by: Appl. Sci. 2017, 7, 716 9 of 16 and C S is given by: where E is the dynamic elastic modulus, G is the dynamic shear modulus, and ν is the Poisson's ratio. The slope is meshed into brick elements [22]. The size of elements has a great influence on the computational results. Kuhlemeyer and Lysmer [24] show that, for accurate representation of wave transmission through a model, the spatial element size must be smaller than approximately one-tenth to one-eighth of the wavelength associated with the highest frequency component of the input wave. In this paper, to balance between accuracy and efficiency, the ratio of the mesh size to the wavelength is equal to 1/12. In addition, the deformational behavior of joints is assumed to be linearly elastic both in the normal and shear directions, and the rock material is assumed to be elastic, isotropic, and homogeneous.
To the viscous boundary, the dynamic input must be applied in a stress or force form, or the effect of the viscous boundary would be nullified. To input seismic motion at a quiet boundary, a stress boundary condition is used (i.e., a velocity or acceleration record is transformed into a stress record and applied to a quiet boundary). A velocity wave may be converted to an applied stress using the formula: or σ s = 2(ρC s ) v s (12) where σ n is the applied normal stress and σ S is the applied shear stress; the factor of two in Equations (11) and (12) accounts for the fact that the applied stress must be double that observed in an infinite medium, since half the input energy is absorbed by the viscous boundary.
In this study, one shear wave a = 0.5gcos(2πt/T) is applied under the bottom of the slope. According to the character described above, the acceleration time history is converted to a stress-time history: σ s = (gT/2π)(Eρ/(2(1 + ν))] 0.5 sin(2πt/T) (13) where T is the period of the wave and g is the gravity acceleration.

Model Properties
The horizontal layered-structure slope is made up of homogeneous, isotropic, and elastic marble [10], with the dynamic modulus listed in Table 2. The material damping is not included in this paper for the same reasons listed in Section 2.1. In order to consider the effect of normalized fracture stiffness (K) and dimensionless fracture spacing (ξ) on the dynamic response of the horizontal layered-structure rock slope, the scheme for the properties of joints and applied stress (σ S ) under the bottom of the slope is listed as Table 3.
Based on the values available in Tables 2 and 3, numerous numerical simulations were conducted to study the effect of normalized fracture stiffness and dimensionless fracture spacing on the dynamic response of the horizontal layered-structure rock slope.

Parametric Study on Normalized Fracture Stiffness K
In order to study the effect of normalized fracture stiffness K on dynamic response of the horizontal layered-structure rock slope, the dimensionless fracture spacing ξ is fixed at ξ = 0.4. Other parameters of the slope are the same as described above. The numerical simulation results of different parts of the slope are described below.
In Figure 8, the horizontal coordinate is the normalized fracture stiffness K, the vertical coordinate is the acceleration amplification coefficient of the slope crest η crest , which is defined by Equation (6). It can be found that η crest decreases with the increase of K, which may be caused by the increase of k s . Figure 9 shows η versus H, where η is the acceleration amplification coefficient of monitoring points on the slope surface and H is the vertical distance of monitoring points on the slope surface to the horizontal coordinate. It can be found that η wavily increases with the increase of H, and the peak values occur at H = 7.5 m, 17.5 m, and 27.5 m in turn, whereas the points near the bottom of joints (the valley values) are at H = 12.5 m and 22.5 m respectively, shown in Figure 9. In order to highlight the area we are most concerned about, only the dynamic response result of the part of the slope above its foot is given. Figure 10 shows the distribution regularity of acceleration amplification coefficient of the slope with different K (K = 0.1, 0.5, 2, and 4, respectively). 3.93 × 10 9 7.86 × 10 9 1.179 × 10 10 1.572 × 10 10 K = 4 5.24 × 10 8 6.56 × 10 8 1.05 × 10 9 2.62 × 10 9 5.24 × 10 9 1.05 × 10 10 1.572 × 10 10 2.096 × 10 10 K = 5 6.55 × 10 8 8.20 × 10 8 1.31 × 10 9 3.28 × 10 9 6.55 × 10 9 1.31 × 10 10 1.965 × 10 10 2.62 × 10 10 σS (Pa) 3.98 × 10 6 3.18 × 10 6 1.99 × 10 6 7.93 × 10 5 4.02 × 10 5 1.95 × 10 5 1.376 × 10 5 1.032 × 10 5 Based on the values available in Tables 2 and 3, numerous numerical simulations were conducted to study the effect of normalized fracture stiffness and dimensionless fracture spacing on the dynamic response of the horizontal layered-structure rock slope.

Parametric Study on Normalized Fracture Stiffness K
In order to study the effect of normalized fracture stiffness K on dynamic response of the horizontal layered-structure rock slope, the dimensionless fracture spacing ξ is fixed at ξ = 0.4. Other parameters of the slope are the same as described above. The numerical simulation results of different parts of the slope are described below.
In Figure 8, the horizontal coordinate is the normalized fracture stiffness K, the vertical coordinate is the acceleration amplification coefficient of the slope crest ηcrest, which is defined by Equation (6). It can be found that ηcrest decreases with the increase of K, which may be caused by the increase of ks. Figure 9 shows η versus H, where η is the acceleration amplification coefficient of monitoring points on the slope surface and H is the vertical distance of monitoring points on the slope surface to the horizontal coordinate. It can be found that η wavily increases with the increase of H, and the peak values occur at H = 7.5 m, 17.5 m, and 27.5 m in turn, whereas the points near the bottom of joints (the valley values) are at H = 12.5 m and 22.5 m respectively, shown in Figure 9. In order to highlight the area we are most concerned about, only the dynamic response result of the part of the slope above its foot is given. Figure 10 shows the distribution regularity of acceleration amplification coefficient of the slope with different K (K = 0.1, 0.5, 2, and 4, respectively).

Parametric Study on Dimensionless Fracture Spacing ξ
In order to study the effects of dimensionless fracture spacing ξ on the dynamic response of the horizontal layered-structure rock slope, the normalized fracture stiffness K is fixed at K = 0.5. The other parameters of the slope are the same as described above. The numerical simulation results of different parts of the slope are shown below.
In Figure 11, the horizontal coordinate is the dimensionless fracture spacing ξ; the vertical coordinate is acceleration amplification coefficients of the slope crest ηcrest. It can be found that the η crest does not monotonously increase with an increase in ξ; it sometimes increases, and sometimes decreases. Figure 12 shows η versus H, where η is the acceleration amplification coefficient of monitoring points on the slope surface. It can be found that η wavily increases with the increase of H, and when H increases to the slope crest, η reaches the maximum. It also can be seen that the fluctuation enhances with ξ increasing, which is contrary to Figure 9.

Parametric Study on Dimensionless Fracture Spacing ξ
In order to study the effects of dimensionless fracture spacing ξ on the dynamic response of the horizontal layered-structure rock slope, the normalized fracture stiffness K is fixed at K = 0.5. The other parameters of the slope are the same as described above. The numerical simulation results of different parts of the slope are shown below.
In Figure 11, the horizontal coordinate is the dimensionless fracture spacing ξ; the vertical coordinate is acceleration amplification coefficients of the slope crest η crest . It can be found that the η crest does not monotonously increase with an increase in ξ; it sometimes increases, and sometimes decreases. Figure 12 shows η versus H, where η is the acceleration amplification coefficient of monitoring points on the slope surface. It can be found that η wavily increases with the increase of H, and when H increases to the slope crest, η reaches the maximum. It also can be seen that the fluctuation enhances with ξ increasing, which is contrary to Figure 9.

Parametric Study on Dimensionless Fracture Spacing ξ
In order to study the effects of dimensionless fracture spacing ξ on the dynamic response of the horizontal layered-structure rock slope, the normalized fracture stiffness K is fixed at K = 0.5. The other parameters of the slope are the same as described above. The numerical simulation results of different parts of the slope are shown below.
In Figure 11, the horizontal coordinate is the dimensionless fracture spacing ξ; the vertical coordinate is acceleration amplification coefficients of the slope crest ηcrest. It can be found that the η crest does not monotonously increase with an increase in ξ; it sometimes increases, and sometimes decreases. Figure 12 shows η versus H, where η is the acceleration amplification coefficient of monitoring points on the slope surface. It can be found that η wavily increases with the increase of H, and when H increases to the slope crest, η reaches the maximum. It also can be seen that the fluctuation enhances with ξ increasing, which is contrary to Figure 9.    Figure 13 shows the distribution regularity of acceleration amplification coefficient of the slope with different ξ (ξ = 0.05, 0.1, 0.2, 0.3, and 0.4, respectively). It can be found that the effects of ξ are obvious to the distribution regularity of acceleration amplification coefficient of the slope. With ξ increasing, the rhythm of contours η is more and more apparent. The rhythm mainly distributes along the joint, as Figure 13 shows. Before the elevation reaches the slope shoulder, the contours are approximately parallel to the oblique slope surface. When they reach the slope shoulder, contours begin to increase along the surface, which is identical to the conclusion drawn by Qi et al. [23]. In addition, the acceleration amplification regularities of the slope in the vertical direction presented in Figure 13 are the same as the conclusions acquired by Wang and Wang [25], and He and Lu [26].  Figure 13 shows the distribution regularity of acceleration amplification coefficient of the slope with different ξ (ξ = 0.05, 0.1, 0.2, 0.3, and 0.4, respectively). It can be found that the effects of ξ are obvious to the distribution regularity of acceleration amplification coefficient of the slope. With ξ increasing, the rhythm of contours η is more and more apparent. The rhythm mainly distributes along the joint, as Figure 13 shows. Before the elevation reaches the slope shoulder, the contours are approximately parallel to the oblique slope surface. When they reach the slope shoulder, contours begin to increase along the surface, which is identical to the conclusion drawn by Qi et al. [23]. In addition, the acceleration amplification regularities of the slope in the vertical direction presented in Figure 13 are the same as the conclusions acquired by Wang and Wang [25], and He and Lu [26].  Figure 13 shows the distribution regularity of acceleration amplification coefficient of the slope with different ξ (ξ = 0.05, 0.1, 0.2, 0.3, and 0.4, respectively). It can be found that the effects of ξ are obvious to the distribution regularity of acceleration amplification coefficient of the slope. With ξ increasing, the rhythm of contours η is more and more apparent. The rhythm mainly distributes along the joint, as Figure 13 shows. Before the elevation reaches the slope shoulder, the contours are approximately parallel to the oblique slope surface. When they reach the slope shoulder, contours begin to increase along the surface, which is identical to the conclusion drawn by Qi et al. [23]. In addition, the acceleration amplification regularities of the slope in the vertical direction presented in Figure 13 are the same as the conclusions acquired by Wang and Wang [25], and He and Lu [26].

Discussion of Results
From the results in Figures 2-4 and 6, we can find that the FLAC results of normally incident P and Sv waves' transmission across a single joint, normal incident P wave transmission through a joint set, and obliquely incident P wave transmission through a single joint are comparable with the analytical solutions. The FLAC results agree well with the analytical solutions. Therefore, FLAC is applicable to study the dependence of wave transmission on the incident angle and the two-dimensional effect.
Then parametric studies on Sv wave transmission through the horizontal layered-structure rock slope were conducted. From Figure 8, ηcrest decreases with the increase of K, which may be caused by the increase of ks. From Figure 3, it can be seen that with the increasing ks the gradient of transmission coefficient T1S decreases. According to the definition of transmission coefficient T1S, the gradient is the acceleration. This can explain why increasing K could lead to ηcrest decreasing. There is a clear change of trend at K = 2, which may be caused by the interaction of joints and stress wave; more studies need to be done. From Figure 9, η presents strong wavy increase. Because of the attenuation effects of joints on wave propagation, once waves transmit across the joint, the value of η will drop to a valley value, and then the topography effects make the values of η increase again. Moreover, the fluctuation weakens as K increases, which is the same as ηcrest versus K. In addition, when H increases to the slope crest, η reaches a maximum, which is consistent with the dynamic response regularities of homogeneous slopes acquired by Qi et al. [10]. There is also a sudden

Discussion of Results
From the results in Figures 2-4 and 6, we can find that the FLAC results of normally incident P and Sv waves' transmission across a single joint, normal incident P wave transmission through a joint set, and obliquely incident P wave transmission through a single joint are comparable with the analytical solutions. The FLAC results agree well with the analytical solutions. Therefore, FLAC is applicable to study the dependence of wave transmission on the incident angle and the two-dimensional effect.
Then parametric studies on Sv wave transmission through the horizontal layered-structure rock slope were conducted. From Figure 8, η crest decreases with the increase of K, which may be caused by the increase of k s . From Figure 3, it can be seen that with the increasing k s the gradient of transmission coefficient T 1S decreases. According to the definition of transmission coefficient T 1S , the gradient is the acceleration. This can explain why increasing K could lead to η crest decreasing. There is a clear change of trend at K = 2, which may be caused by the interaction of joints and stress wave; more studies need to be done. From Figure 9, η presents strong wavy increase. Because of the attenuation effects of joints on wave propagation, once waves transmit across the joint, the value of η will drop to a valley value, and then the topography effects make the values of η increase again. Moreover, the fluctuation weakens as K increases, which is the same as η crest versus K. In addition, when H increases to the slope crest, η reaches a maximum, which is consistent with the dynamic response regularities of homogeneous slopes acquired by Qi et al. [10]. There is also a sudden increase at about H = 27.5 m, caused by the strong amplification effect at the slope crest. From Figure 10, it can be found that the effects of K are not apparent from the distribution regularity of acceleration amplification coefficient of the slope, and the contour of η takes on rhythm, which is similar to the dynamic response of the homogenous slope made by Qi et al. [27]. From Figure 11, we can see that η crest is more sensitive to ξ, compared to K. This phenomenon may be caused by wave frequency, as can be seen in Table 3, but further study is needed. From Figure 12, it can also be seen that the fluctuation enhances as ξ increases, which is contrary to Figure 9. The enhancement of fluctuation may be caused by the same factors described in Section 4.1, which is a result of the combined action from attenuation effects of joints and topography effects.

Conclusions
For a rock slope with multiple joints, it is complex to describe its dynamic response regularities under stress wave loading. One reason is that many factors can influence the dynamic response of the slope, such as slope angle, slope height, properties of joints, properties of applied stress wave, and so on. In this paper, the simplified configuration of a rock slope with a horizontal layered-structure is given. Therefore, this method makes it possible to focus on the study of the effect of the parameters in which we are most interested (K and ξ). Meanwhile, other parameters are fixed.
This study focused on two issues: one is the verification of FLAC modeling on wave transmission across rock joints, and the other is the study of parametric studies and dynamic response on Sv wave propagation through a horizontal layered-structure rock slope. The existence of different forms of joints in a natural rock slope usually makes the problem more complicated. However, the given numerical slope model simplifies this issue. Before the simulation of the dynamic response of the slope, a verification of FLAC modeling on wave transmission across rock joints is executed. Through comparison with analytical solutions and FLAC results, the effects of joint mechanical and spatial properties on wave transmission in FLAC modeling are validated. Through the numerical simulation study of the dynamic response of a horizontal layered-structure rock slope using FLAC, some conclusions about the attenuation effects of joints and topography effects can be drawn as follows:

•
The acceleration amplification coefficient η decreases with the increase of the normalized fracture stiffness K, and the decrease rate is obviously reduced once K is larger than 2; • The acceleration amplification coefficient η does not increase with the increase of the dimensionless fracture spacing ξ, and there is a decrease range due to the superposition of the reflection and the refraction from the joints; • The effect of the dimensionless fracture spacing ξ is more obvious on the acceleration amplification coefficient η than that of the normalized fracture stiffness K.