Investigation of the seismic response characteristics of a rock mass slope containing weak structural planes under seismic excitation based on multi-domain coupling analysis

Abstract Because of the complexity of interactions between seismic waves and rock masses, the dynamic response of slopes has typical multi-domain characteristics. To fully consider the multi-domain characteristics of ground motion, the dynamic response of a rock slope containing weak structural planes was systematically investigated in the time domain, frequency domain, and time-frequency domain. Two-dimensional dynamic analyses were carried out on two numerical models, including the homogeneous slope- and slope-containing structural planes. The results show that slopes containing structural planes have dynamic topographic and geological amplification effects. Based on the frequency-domain analysis, the relationship between the frequency of waves and the dynamic response of slopes was further studied. The geological structure influences the high-frequency components of waves in slopes containing structural planes. The analyses of the marginal spectrum show that seismic energy is mainly concentrated in the frequency band of waves (7–10 Hz). In addition, the seismic failure mechanism of slopes containing structural planes was explored. Structural planes have an important effect on the seismic failure mode of slopes containing structural planes. Multi-domain coupling analysis can reveal the seismic response characteristics of slopes more comprehensively and provide a reference for the seismic design of rock slopes.


Introduction
The landscape of Western China is mainly mountainous, and earthquakes have induced many landslides (Zhou et al. 2017;Chen et al. 2020;. Field investigations after the 2008 Wenchuan earthquake show that seismic wave propagation in rock-soil bodies weakens the stability of slopes to a large extent (Yin et al. 2015). The spatial distribution of earthquake-induced landslides after the earthquake also reveals that seismic landslides have specific topographic and geological amplification effects (Wang et al. 2010;Huang et al. 2012). The influence of earthquakes on the stability of landslides has attracted the attention of many scholars for several years (Chen and Wu 2018;Chen, Song, Juliev, et al. 2021).
With long-term geological action, there are many discontinuous joints in the rock mass, including penetrating structural planes that cut the rock mass into discontinuous bodies (Wu et al. 2014;Du et al. 2020;Guo and Cui 2020;Cao et al. 2021). The discontinuous geological structure in the rock mass weakens its strength and stability (Liu et al. 2020;Song, Liu, Chen, et al. 2021). The failure modes of toppling slopes containing discontinuities mainly include bend toppling, block toppling, and massive bending and toppling (Alejano et al. 2010;Amini et al. 2012). The dynamic failure modes of rock slopes are affected by discontinuity distributions (Zhou et al. 2017;. Due to the randomness of earthquakes and the complexity of rock mass structural materials, the stability of rock slopes is difficult to fully understand. Hence, the dynamic stability of rock slopes needs to be further studied.
Many scholars have studied the seismic dynamic response characteristics of rock slopes (Liu et al. 2013;, Song, Liu, Li, et al. 2021. Time-domain analysis is the common method used to investigate the dynamic response of slopes (Fan et al. 2017;Li et al. 2018). The dynamic acceleration response is used to study the seismic response characteristics of various types of slopes, such as jointed slopes (Li et al. 2018), rock slopes with weak structural surfaces , and layered slopes ). In addition, frequency-domain analysis can explore the relationships between seismic wave frequency, natural frequency, and the dynamic response characteristics of slopes (Song, Liu, Li, et al. 2021). Moore et al. (2011) investigated the dynamic response of slopes using the Fourier spectrum response of large rock slopes. Li et al. (2019) investigated the dynamic response of jointed slopes and the damage characteristics of slopes by analysing the Fourier spectrum. He et al. (2020) investigated the seismic stability of a landslide by analysing the Fourier spectrum response of a strong earthquake seismograph array. The slope instability caused by the resonance between the high frequency of waves and the natural frequency of slopes has always been a hot issue in the field of seismic engineering. Yang et al. (2015) used the Hilbert-Huang transform (HHT) to study the dynamic response of a rock landslide and its triggering mechanism from an energy perspective, which showed that time-frequency domain analysis could better reflect the dynamic response of slopes. Fan et al. (2016Fan et al. ( , 2017Fan et al. ( , 2019 used the HHT marginal spectrum to study the seismic failure mode of slopes with structural planes, suggesting that the marginal spectrum could better reflect the seismic failure development process in slopes. The dynamic response of rock slopes with complex geological structures has typical multi-domain characteristics. Time domain, frequency domain, and time-frequency domain analyses have become important methods for analysing the seismic response of slopes.
Seismic waves are typical random waves containing complex frequency components . Due to the heterogeneity and nonlinearity characteristics of rock masses, the seismic response characteristics of rock masses are difficult to systematically understand. The dynamic characteristics of slopes should be further clarified from the perspective of the frequency domain. Moreover, due to the time-frequency-amplitude characteristics of seismic waves, it is necessary to further reveal the seismic response of slopes in the time-frequency domain. However, previous research focussed on the time-domain analysis of the seismic response of rock slopes but ignored the frequency domain and time-frequency domain characteristics. Therefore, to solve the above problems, it is necessary to establish a multi-domain coupling analysis method to systematically reveal the seismic response characteristics of complex rock slopes considering the three domains.
A multi-domain coupling analysis method is proposed to investigate the seismic response characteristics of a rock mass slope containing weak structural planes, as shown in Figure 1. Two-dimensional dynamic analyses are performed on two numerical models using FLAC 3 D (Fast Lagrangian Analysis of Continua). This work focuses on three aspects in this research: dynamic acceleration response (time-domain analysis), modal analysis and Fourier spectrum analysis (frequency-domain analysis), and marginal spectral analysis of the models (time-frequency domain analysis). The influence of topographic and geological conditions on the seismic response of the models was analysed. In addition, this work analyses the relationship between the natural frequencies and the dynamic characteristics of slopes. The dynamic failure mechanism of a slope containing structural planes is discussed as well.

Case study
The rock slope is located in the hilly region of Sichuan Province in Western China, where the river system, gullies, and valleys are developed. The landform in the region is dominated by mountains, hills, and valley plains ( Figure 2). The strata in the slope area are continuous, with gentle occurrence and no faults passing through. Under the influence of regional tectonic stress, deadweight stress, and unloading, rock joints in the area are very well developed. Several seismic fault zones are distributed near the slope. The slope body with weak structural planes is approximately 34 m long and 40 m high, with a gradient of approximately 70 . The development of the slope has multilayer weak structural planes, joints, and fissures. Weak structural planes are developed on the contact surface between silty mudstone layers. The structural planes are distributed in a zonally continuous manner on the slope. To obtain the physical and mechanical parameters of the rock mass in the slope, a series of laboratory tests, such as direct shear tests, triaxial tests, compression tests, and tensile tests, were performed (Table 1). According to the geological structure type of the slopes, the geological model of the slope is simplified, as shown in Figure 3.

Numerical model
The rock mass of the model is elastoplastic material, and the Mohr-Coulomb strength criterion is adopted. The calculated boundaries of the model are as follows: the length from the foot of the slope to the right boundary is twice the length of the slope, the length from the top of the slope to the left boundary is twice the height of the slope, and the height from the top to the bottom of the slope is twice the height of the slope. The boundary range meets the requirements of calculation accuracy under static and dynamic conditions (Xu et al. 2008). Model 1 (homogeneous slope) and Model 2 (slope containing weak structural planes) were established, and the two numerical mesh models, with dimensions of 80 Â 170 m, are shown in Figure 4. To reduce the difficulty of simulation, the structural planes were simplified into a weak band with a depth of 0.2 m. A quadrilateral grid with a side length of 0.5 m was used to investigate the characteristics of the model. FLAC 3 D was used as a solver for dynamic analysis via its transient dynamics simulation capability. Local damping is  used in the model. In dynamic analysis, model boundary processing is a key technology. An improper boundary leads to wave reflection and refraction, which seriously affect the results of dynamic analysis. In this model, the discontinuity is simulated as a soft material rigidly connected with the rock mass, and the selection of the thickness of the discontinuity is more critical. The range of the discontinuity should be greater than the wavelength, and the thickness should be less than the wavelength. For the dynamic response analysis of a semi-infinite space body such as a slope, the boundary problem that tends to infinity must be dealt with. FLAC provides a static boundary (viscous boundary) and free field boundary in the numerical simulation calculation. The research object in this work is a rock slope whose foundation modulus is large and can be considered a rigid foundation. Therefore, static boundary conditions are not required at the bottom of the model. The free field boundary is set around the model so that the side boundary of the main grid is coupled with the free field grid through dampers. The unbalanced force of the free field grid is applied to the boundary of the main grid. The boundary conditions are set as shown in Figure  3b. To study the dynamic response of the slope, some monitoring points are set up in the model. To eliminate the adverse effect of the deadweight of the slope, a stress balance should be carried out before the numerical simulation.
Two loading directions of Wenchuan earthquake waves (the z and x directions) were applied at the bottom boundary of the model, and the corresponding time history curve and Fourier spectrum are shown in Figure 5. The Wenchuan earthquake (WE) waves were recorded at Wudu in Gansu, China. The dominant frequency wave is 7.74 Hz, t ¼ 120 sec, and Dt ¼ 0.005 sec to create time-dependent acceleration inputs at the bottom of the slope during horizontal and vertical motions in the dynamic calculations.
3. Multi-domain coupling analysis of the dynamic response characteristics of slopes

Analysis of the time domain
To analyse the influence of structural planes on the wave propagation characteristics, taking the input wave (0.1 g) as an example, the acceleration distribution of the models at different times was extracted. The acceleration distributions of the two models are shown in Figures 6-9. Waves mainly propagate along the slope surface from the bottom to the slope crest. Comparing Model 1 and Model 2 reveals that the acceleration distribution characteristics of Model 2 have apparent changes. A noticeable change in the phase difference can be found on both sides of the weak planes. Furthermore, the local magnification amplification effect of the toppling slope can be identified, which indicates that structural planes affected the propagation characteristics of waves inside the slope. This effect occurs because the existence of structural planes causes refraction and reflection effects of waves when waves propagate in the rock mass so that waves appear superposed near weak planes, further resulting in a local amplification effect of waves. The peak ground acceleration (PGA) reflects the strongest response in the whole acceleration time history and reflects the maximum seismic inertial force at a certain position in the slope body. The M PGA is used to analyse the acceleration response of slopes. The ratio of PGA at any point on the slope to PGA at the foot of the slope (A1 point) was defined as M PGA , which represents acceleration magnification at a point of the slope. Four typical monitoring points (A1, A3, A15, and A7) of the slope surface are selected, and their acceleration-time histories and Fourier spectra curves  at the slope crest. The M PGAmax values are approximately 1.2 and 1.14 when the input waves are in the x-and z-directions, respectively. This shows that the acceleration amplification effect of the slope increases with elevation, and the slope has a dynamic elevation amplification effect. The M PGA of Model 2 increases with elevation as well, indicating that the slope also has a typical elevation dynamic amplification effect (Figure 12b). However, comparing Figure. 12a and b, the M PGA of Model 1 shows an increasing linear trend with elevation overall, while Model 2 shows a significant nonlinear increase, indicating that the structural planes influence the dynamic magnification effect of the model. This effect occurs because structural planes cause the discontinuity of the rock mass, and the refraction and reflection effects of waves appear near the weak structural plane, which leads to the nonlinear variation characteristics of the amplification effect. The M PGA at the slope surface is significantly larger than that inside the slope, suggesting that the slope has a significant amplification effect on the slope surface. Moreover, the M PGA of Model 2 is approximately 1.1 times that of Model 1, which suggests that the structural plane has a specific amplification effect on the slope. In addition, in Model 1 (Figure 12a), the M PGAmax under horizontal seismic force is approximately 1.05 times that under vertical seismic force. In Model 2 (Figure 12b), the M PGAmax under horizontal seismic force is approximately 1.15 times that under vertical seismic force. This finding indicates that the directions of seismic excitation influence the dynamic response of the slopes, and the acceleration amplification effect under horizontal earthquake excitation is greater. Fan et al. (2016) carried out a laboratory shaking table scaled model test on Model 2, and the comparison between the model test and the numerical test results is shown in Figure 13. The acceleration amplification factor obtained from the model test is similar to the numerical simulation result, indicating that the numerical result is reliable. In the numerical model, the size and physico-mechanical parameters of the   compared to the original slope. During the construction of the test model, there are some differences in the adhesion between the simulated rock mass and the structural plane, which is different from the original slope. These factors directly lead to the difference between the numerical results and the model test results. However, in general, the numerical results and experimental results can better reflect the dynamic response patterns of the slope, and the results of both analyses confirm each other.

Analysis of the frequency-time domain
The frequency-domain analysis mainly includes modal analysis and Fourier spectrum analysis. The principles of FFT and modal analysis are as follows. Modal analysis has become widely used type of dynamic frequency-domain analysis (Reale et al. 2016;Song, Liu, Li, et al. 2021). Modal analysis can reveal slopes' natural frequency and vibration mode and predict the dynamic response of slopes in the elastic domain. It is worth noting that the low-order vibration mode has a controlling effect on the dynamic characteristics of the engineering body. Therefore, in modal analysis, only the first few modes are considered. In addition, FFT is used to analyse the vibration response of rock and soil masses, which can identify different frequency components of signals. FFT can quickly identify the main components of the signal and can also be quickly filtered, which has become a common method to deal with seismic wave signals.
The modal analysis results of the models are shown in Figures 14-16. Figure 14 shows that the natural frequencies (f) of the slopes increase with the order. The first four f values of Model 1 are 2. 96, 7.84, 11.18, and 18.06 Hz, respectively. The corresponding f values of Model 2 are 2. 81, 7.34, 10.69, and 17.63 Hz, respectively. The natural frequency of Model 2 is less than that of Model 1 because Model 2 contains multiple structural planes, which decreases the stiffness of the rock mass. This finding suggests that structural planes affect the natural frequency of the models. In Figures  15 and 16, U refers to the relative displacement of the models under a certain mode of vibration. Figures 15 and 16 show that the first four modes of Model 1 are similar to Model 2 overall. However, U near the structural planes of Model 2 has a discontinuous distribution, with a certain phase shift, which indicates that the structural plane affects the dynamic response of the rock mass. Figures 15 and 16 show that the U max of Model 2 is larger than that of Model 1; in particular, the higher-order modal characteristics of the slopes are different. This difference arises because structural surfaces enlarge the slope deformation. The modal characteristics of different orders show that the U of the slope increases with the slope height, on the whole, suggesting that the slopes have a pronounced elevation amplification effect and are prone to damage on the slope crest. Moreover, Figures 15 and 16 also show that under the same elevation condition, the U of the slope surface is significantly larger than that inside the slope body, which is similar to the time-domain results.
In addition, the dynamic response of the slope is further analysed by using Fourier spectrum analysis. Figures 10 and 11 shows that Model 1 comprises three predominant frequencies (f 1 -f 3 ), and Model 2 has four predominant frequencies (f 1 -f 4 ). The values of the predominant frequencies are similar to those in the modal analysis. Fourier spectrum analysis and modal analysis further determined the f of the slopes. It is worth noting that Model 2 has f 4 , which indicates that structural planes have a noticeable magnification effect on the high-frequency components of waves. Figures  17 and 18 show that the PFSA of the slopes increases with elevation at different natural frequencies. Compared with the PFSA of the internal slope body, the slope surface is significantly larger. This indicates that the slope has prominent elevation and slope surface magnification effects, consistent with the analysis results in the time domain. Moreover, the PFSA of Model 2 is approximately 1.3 times that of Model 1.
The PFSA of f 1 and f 2 shows an increasing linear trend, while that of f 3 and f 4 shows a nonlinear increasing trend overall, suggesting that the structural plane affects the high-frequency components of waves. Therefore, modal analysis and Fourier spectrum analysis can identify the relationship between the frequency components of waves and the amplification effect of slopes.

Analysis of the time-frequency domain
The HHT method is mainly composed of the empirical mode decomposition (EMD) method and Hilbert spectrum analysis (HSA), which provides a feasible tool for identifying earthquake Hilbert energy distributions (Huang et al. 1998;Fan et al. 2017). Using the EMD, the seismic wave is decomposed into a series of intrinsic mode functions (IMFs). The flowchart for the EMD approach is shown in Figure 19 (Zhang 2006). The HHT marginal spectrum is defined as the integral of the Hilbert Huang can show the accumulation of IMF amplitude in the whole period. Every instantaneous frequency has a certain amount of energy. By accumulating the energy of these instantaneous frequencies, the total energy corresponding to the frequency in the original signal can be calculated, that is, the marginal spectrum amplitude (MSA). The marginal spectrum can be used to measure the contribution of each frequency value to the total Hilbert energy and corresponds to the Hilbert energy density at a certain frequency. The Hilbert marginal spectrum obtained by the HHT represents the distribution of the seismic signal energy amplitude along the frequency axis.
The marginal spectrum shows the energy distribution of the slope vibration, and the instantaneous spectrum can reflect the time-varying characteristics of the seismic wave signal. The propagation characteristics of seismic energy can be determined by  analysing the variation in marginal spectral amplitude. Taking the acceleration time history of point A1 in Model 2 as an example, the first four IMFs and instantaneous frequency are obtained by EMD, as shown in Figures 20 and 21. Then, the marginal spectrum of the IMFs can be obtained by using HHT. As shown in Figures 20 and  21, the amplitude of IMF2 is larger, and its frequency component is much richer, which makes it easier to identify. Hence, IMF2 was used to obtain the marginal spectrum. The marginal spectrum of typical measurement points in Model 2 is shown in Figure 22. The peak marginal spectrum amplitude (PMSA) is mainly located in the range of 7-10 Hz. This frequency band is similar to f 2 and f 3 of Model 2, so the seismic wave energy in this frequency band is more suitable for analysing the dynamic response characteristics of the slope.
The PMSAs of the models with elevation are shown in Figure 23. The PMSA of the slopes increases with elevation, and the PMSA of the slope surface is larger than that of the internal slope. This phenomenon indicates that the slope surface and elevation magnify the energy propagation of seismic waves to a certain degree. Moreover, Figure 23 shows that the PMSAz when the input wave in the z-direction is significantly less than the PMSAx under horizontal seismic loading. The PMSAx is approximately 1.1 times that of PMSAz in Model 1, and the PMSAx of Model 2 is approximately 1.15 times that of PMSAz. This indicates that the direction of the input wave influences the energy propagation of seismic waves and the geological amplification effect when the input wave in the x-direction is larger. The PMSA of Model 2 is approximately 1.2 times that of Model 1. This difference shows that the structural plane magnifies the propagation energy of seismic waves in the rock mass. Figure 19. Flowchart for EMD approach (Zhang 2006). Therefore, the dynamic response of slopes can be investigated from an energy-based perspective in the time-frequency domain.

Dynamic deformation characteristics of the slope
The dynamic deformation characteristics and evolution process of the slope can be investigated by analysing the distribution of shear strain increments of the slope under earthquakes. The contour of the shear strain increment of the slope under different ground motion intensities is shown in Figure 24. Figure 24 shows that the maximum value of the shear strain increment appears in the weak structural planes, while the value of the shear strain increment in the rock of the slope body is small, which indicates that weak structural planes are the most prone to shear failure subject to earthquake excitation; that is, the structural planes are the potential slip planes of the slope. Figure 24 shows that at the initial stage of ground motion, the shear strain increment of the topmost structural surface is large under small earthquakes ( Figure  24a, b), which indicates that the top of the slope is the area most prone to damage. With the increase in ground motion, the increment of shear strain gradually transfers to the lower part of the slope body (Figure 24c, d), which indicates that the sliding surface gradually expands to the bottom of the slope as the shear strain in the sliding range of the slope continues to expand with the increase in ground motion. These results of the damage process of the slope are consistent with the results of the shaking table test ( Figure 25) . Therefore, the structural planes greatly impact the failure mode of the slope, which are the potential sliding planes.

Conclusions
A numerical dynamic analysis method is used to study the dynamic response of the two models using the multi-domain coupling analysis method. Time-domain analysis can directly and preliminarily reflect the dynamic response characteristics of slopes. Frequency-domain analysis can clarify the relationship between the wave frequency component and slope dynamic response. Time-frequency analysis can reveal the dynamic response of slopes from an energy-based perspective. The following conclusions can be drawn.
1. The results of the time-frequency joint analysis show that Model 2 has a noticeable topographic and geological dynamic amplification effect. Structural planes have an impact on the wave propagation characteristics and amplification effect of the slope. The amplification effect of the slope is increased due to the distribution of structural planes. The M PGA , PFSA and PMSA of Model 2 are approximately 1.1, 1.3 and 1.2 times those of Model 1, respectively. Wave propagation directions have an impact on the dynamic response of the slope. The M PGAmax under horizontal seismic loading is approximately 1.15 times that under vertical seismic loading.
2. According to the frequency-domain analysis, the natural frequency impacts the dynamic response of the slope by using the Fourier spectrum and modal analysis. Structural planes have impacts on the dynamic deformation characteristics of the slope, which have an amplification effect on the high-frequency components of waves. The natural frequencies of Model 2 are smaller than those of Model 1.
3. According to the time-frequency joint analysis, in combination with the contour of the shear strain increment of Model 2, weak structural planes are potential slip surfaces and greatly affect the seismic mode of Model 2. The upper slope body is prone to deformation under small earthquakes. With the increase in earthquakes, the range of sliding bodies gradually extends to the lower slope. However, the energy method proposed in this work has some limitations in analysing the seismic response characteristics of rock slopes with discontinuities. The marginal spectrum is obtained based on a specific order of IMF without considering the distribution of seismic Hilbert energy in the time domain; hence, it is difficult to determine the specific time when the slope is damaged. Further research can apply the energy method to the evaluation of seismic dynamic response characteristics and failure evolution process of rock slopes with complex geological structures.