Numerical Investigation on Dynamic Response Characteristics and Deformation Mechanism of a Bedded Rock Mass Slope Subject to Earthquake Excitation

: In order to systematically reveal the dynamic response characteristics of rock mass slopes subject to seismic excitation, time-domain and frequency-domain analyses are used to investigate the dynamic response of a bedded rock slope from multiple perspectives, using the two-dimensional numerical dynamic analyses. Based on the numerical simulation results, the inﬂuence of the weak bedded structural planes on the propagation characteristics of seismic waves in the slope is analyzed. The time-domain analysis suggests that the topographic and geological conditions have an inﬂuence on the dynamic response of the slope. The effects of ground motion direction on the dynamic response characteristics of the slope are identiﬁed. In addition, according to the frequency-domain analysis, the impacts of slope surface, elevation, and structural plane on the seismic response characteristics of the slope are also clariﬁed. The intrinsic characteristics of the slope are investigated by using Fourier spectral analysis and modal analysis, and the deformation response characteristics of the slope are clariﬁed. The relationship between different natural frequencies of the slope, the predominant frequency of the seismic wave, and the dynamic response characteristics of the slope is discussed. Moreover, the dynamic failure mechanism of the slope is analyzed. This work provides a reference for the seismic analysis of this type of slope.


Introduction
Rock slope stability remains an important issue in geotechnical engineering, as the slope is the most common geotechnical infrastructure [1][2][3]. The failure of slopes is often involving in large deformations and results in many natural disasters, e.g., landslides, debris flow, and water flooding [1]. Recent earthquakes, including the 2008 Wenchuan earthquake in China, have demonstrated that the sliding and collapse of a slope during such large-scale events represent significant seismic hazards [4][5][6]. A surprising consistency has been observed between the global distributions of earthquake zones and landslide zones [7][8][9]; this finding sufficiently suggests that earthquakes have a substantial impact on the dynamic stability of slopes. The rock slopes are very prone to form landslides and collapses induced by external factors, such as subsequent earthquakes or rainfall, for example, the 2017 Xinmo village landslide in Sichuan Province, China [10,11]. In populated mountainous areas, earthquake-induced landslides and rockfalls are considered considerable natural threats to the daily lives of people [12][13][14].
Bedded rock slopes are the common type of geologic body in the southwest of China, whose dynamic stability greatly affects the surrounding infrastructure [15,16]. Under the long-term geological action, a large number of discontinuous joints are generated in rock mass [17][18][19]. The failure modes of bedded slopes containing discontinuities are as shown in Figure 1, including sliding-tension failure, sliding-shear failure, bending-tension failure, and sliding-bending failure [20,21]. Weak structural planes have an important influence on the dynamic characteristics and deformation and failure modes of rock slopes [22][23][24]. Due to the discontinuity of the structural planes in the rock mass, the interaction mechanism between the structural planes and seismic waves becomes very complex, thus making it difficult to fully understand the seismic response characteristics and failure modes of the rock slope [25,26]. Hence, the seismic stability of bedded slopes has become very complicated, and thus, focused research should be conducted on the seismic failure mode of bedded slopes. Generally, the seismic acceleration of the time history of slopes provides the basic information to evaluate their dynamic response [27,28]. The time-domain analysis can directly reflect the dynamic response law of a slope from the perspective of mechanics and deformation. The PGA is readily available and can be used to easily and efficiently analyze the changes in the amplification effect of slopes under earthquake excitation [28,29]. The existing research about the seismic response of rock slopes mostly use the time-domain analysis [22][23][24]. Many scholars have used PGA to study the seismic response of different types of rock slopes, such as bedded slope, anti-dip slope, slope with multiple lithologic combinations, slope with discontinuous joints, etc. [13][14][15][16][17]26]. It is worth noting that the existence of discontinuity in rock mass directly impacts the dynamic mechanics of slopes [30]; in particular, weak structural planes can result in the great change of frequency components of seismic waves in rock masses. However, the time-domain analysis using dynamic acceleration response cannot fully reflect the dynamic response characteristics of complex geological rock slopes [31,32].
In addition, the dynamic response of slope is the result of the interaction of different frequency components of seismic waves. The frequency domain analysis can study the influence of ground motion parameters, topography, and geological conditions on the dynamic response law of slope from the perspective of the frequency domain. Previous studies have shown that frequency domain analysis can reveal the relationship between the superior frequency of the seismic wave, the natural frequency of rock slope, and its dynamic response characteristics from a deep level [31,32]. At present, slope instability caused by the resonance between the superior frequency of the seismic wave and the natural frequency of rock and soil mass is, and has always been, a hot issue in the field of seismic engineering [33]. In addition, as an important part of frequency-domain analysis, modal analysis is mainly used to obtain the vibration mode and natural frequency of the engineering entity, which can effectively reflect the relationship between the natural frequency of the engineering entity and its deformation response characteristics [34][35][36]. However, modal analysis has not been widely used in geotechnical engineering. Therefore, the frequency domain research content of rock slope can be further improved by incorporating the modal analysis into the frequency domain research of rock slope and fully considering the mutual verification and supplement of Fourier spectrum and modal analysis results. Therefore, the frequency-domain analysis should be paid much more attention to further study the dynamic response characteristics of slopes.
In this work, the time-domain and frequency-domain analyses are fully considered to reveal the seismic response characteristics of a layered rock mass slope from multiple perspectives. Taking a bedded rock slope as an example, two-dimensional dynamic analyses are performed on the slope, using FLAC3D. Two aspects in this research were used to investigate the seismic response of the slope, including the time-domain analysis using dynamic acceleration response and frequency-domain analysis using modal analysis and Fourier spectrum analysis. The wave propagation characteristics through the rock mass slope and topographical dynamic amplification effect of the slope were analyzed. The impact of weak structural planes on the seismic response characteristics of the bedded slope was also identified. In particular, the relationship between the natural frequencies of the slope and its dynamic deformation characteristics was clarified under earthquakes. Moreover, the dynamic failure mechanism of the slope was discussed, based on the timeand frequency-domain analyses.

Case Study
The rock slope is located in the hilly region of Sichuan Province in western China. The landform in the region is dominated by mountains, hills, and valley plains (Figure 2a). The overall topography of the slope area is high in the southwest and low in the northeast. There is a large amount of rainfall in the area, and the surface runoff has a significant influence on slope surface erosion. The strata in the slope area are continuous, with gentle occurrence and no fault passing through. Under the influence of regional tectonic stress, deadweight stress, and unloading, rock joints in the area are very developed. Several seismic fault zones are distributed near the slope. Slope instability in the study area is affected by discontinuity characteristics of the rock mass [2]. The geomechanical analysis was carried out at more than ten measuring stations in the study area, and the characteristics of nodes were collected by the scanning line method. Point loading tests were carried out on several representative rock samples at each measuring point in the laboratory. According to the method proposed by Hoek et al. (2013) [36], the geological strength index was used to estimate the rock mass properties at all survey points. The slope with a weak interlayer is about 70 m long, 76 m high. The microgeomorphology of the slope surface is characterized as alternately steep and gentle, and the upper platform of the slope crest is developed. The lithological composition of the slope mainly includes silty mudstone and argillaceous silty siltstone, and quaternary residual silty clay. The development of the slope has multilayer weak structural planes, joints, and fissures. The main component of the weak structural plane is clay, which mostly develops on the interface of silty mudstone and argillaceous silty siltstone. The weak structural plane is characterized by continuous zonal distribution, and its thickness is about 3-30 cm [23]. The thickness of the rock layer is about 5-8 m. A typical geological profile of the slope in the study area is shown in Figure 2b. The physical and mechanical parameters of rock are obtained through a series of experimental rock mechanics experiments. The cohesion and internal friction angle of rock are obtained by the direct shear test. The elastic modulus and Poisson's ratio of rock are obtained by the uniaxial compression test. The material parameters of slope rocks and weak structural surfaces are shown in Table 1 [13,23]. According to the geological structure type of the slope, the geological model of the bedded slope is simplified, as shown in Figure 3.

Numerical Calculation Model
Dynamic analyses were used to simulate the wave propagations through rock masses by using finite difference software FLAC3D. According to the geological survey data, the simplification numerical model is shown in Figure 3. The calculated boundary of the model is 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 [37,38]. The numerical mesh model was modeled as shown in Figure 4, and the model size is 80 × 170 m. In the numerical simulation of the rock mass stress field, the weak structural plane is often assumed to be a weak zone with weaker deformation than the surrounding rock, that is, the corresponding mechanical parameters are reduced to make it easier to deform. These parameters of weak structural planes can be determined by rock mechanics tests. In the numerical test, the structural planes were simplified into a weak band with a depth of 0.2 m. Their density, Poisson ratio, elastic Modulus, friction angle, and cohesion are 18.4 kN/m 3 , 0.35, 35 MPa, 23 • , and 14 kPa, respectively. In addition, a quadrilateral grid and square grid with a side length of 0.5 m were used to study some qualitative characteristics of behavior in the rock. The density, Poisson ratio, elastic Modulus, friction angle, and cohesion of rock are 2400 kN/m 3 , 0.16, 375 MPa, 35 • , and 1200 kPa, respectively. Next, 30,756 surface strain elements and 31,162 grid nodes were generated in the model. The damping in the numerical dynamic analysis module includes local damping, Rayleigh damping, and viscous damping. Local damping converges in the vibration cycle by increasing or decreasing mass at nodes or structural element nodes. Local damping was adopted in the model.

Boundary Condition Setting
For the dynamic response analysis of a semi-infinite space body, such as a slope, the boundary problem, which actually tends to infinity, must be dealt with [38]. FLAC3D provides the static boundary (viscous boundary) and free field boundary in the numerical dynamic analysis and solves the boundary problem in the dynamic analysis. 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 were not required at the bottom of the model. The free field boundary was set at both sides of the model so that the side boundary of the main grid was coupled with the free field grid through dampers. The unbalanced force of the free field grid was applied to the boundary of the main grid. Since the free field boundary provides the same effect as the infinite field, the upward surface wave does not distort. The boundary conditions were set as shown in Figure 3. To investigate the seismic response of the slope, some monitoring points were set up in the model. The rock mass is an elastoplastic material, and the Mohr-coulomb strength criterion was adopted. The boundary is the free field boundary. Local damping was adopted in the dynamic analyses, whose damping coefficient is 0.156. The static calculation was performed before the dynamic calculation [39]. To avoid the effect of rock gravity, the stress balance should be done before the dynamic analysis.

Seismic Wave Loading
The horizontal and vertical loading directions of the Wenchuan Earthquake (WE) waves were loaded at the bottom boundary of the model; the acceleration time history curves and the corresponding Fourier spectrum are shown in Figure 5. The WE waves were recorded by the Wudu seismic station in Gansu province, China. The dominant frequency wave is 7.74 Hz, t = 120 s, and ∆t = 0.005 s.

Analysis of Wave Propagation Characteristics in the Slope
To clarify the impact of the structural planes on the wave propagation characteristics through the slope, taking the input WE wave in xand z-directions (0.1 g) as examples, the acceleration distribution of the slope during the wave propagation process from the bottom slope to the slope crest are shown in Figures 6 and 7. The acceleration results from the calculations considering the initial stress due to gravity. Figure 6 shows an obvious layered characteristic of wave propagation in the bedrock mass, but an irregular propagation characteristic can be identified through the structural planes inside the slope. An obvious acceleration phase shift on both sides of the structural planes can be identified in Figure 6d-f, indicating that structural planes have a magnification effect on the wave propagation characteristics of the slope. Moreover, it is worth noting that, as shown in Figure 7, seismic waves in both the slope bedrock area and the slope body under the action of vertical seismic force show irregular propagation characteristics. In particular, at the junction of the slope surface and bedded structural planes, there is an obvious accelerated phase shift, which is manifested as an obvious local amplification effect. Moreover, the acceleration increases with the elevation during the wave propagation process, and reaches the maximum at the slope crest, indicating that the slope owns an elevation magnification effect subject to earthquake excitation.

Analysis of Acceleration Magnification Effect of the Slope
To investigate the seismic response characteristics of the slope, taking the WE wave (0.1 g) when input in x-direction as an example, the acceleration-time history of the points is shown in Figures 8 and 9. According to Figure 8, the PGA (peak ground acceleration) of the slope surface shows an increasing trend with the elevation. Figure 9 also shows that the PGA increases with the elevation inside the slope. This phenomenon suggests that the elevation is positively correlated with the PGA. To further investigate the seismic response of the slope, M PGA is introduced to define the ratio of PGA at any point in the slope and PGA at the slope toe (point A1) as M PGA , which represents acceleration magnification at certain a point of the slope. The change rules of M PGA with the slope relative elevation (h/H) are shown in Figure 10  The dynamic acceleration magnification effect of the slope has an increased effect with the elevation, indicating that the slope has an elevation amplification effect. However, the M PGA shows an obvious nonlinear variation characteristic, which is caused by the discontinuity of the rock mass, and wave refraction and reflection effects appear near the structural planes, resulting in nonlinear variation characteristics of the amplification effect. By comparing the M PGA of the internal slope with that of the slope surface, it can be found that the M PGA of the slope surface is significantly larger than that of the internal slope. For example, when the elevation is about 0.6 m, M PGA of the points A5 and A6 are 1.27 and 1.02, respectively, which indicates that the slope has an obvious slope surface amplification effect. The slope surface is the free surface, which results in a great amplification effect on the seismic energy and indicates that the surface slope amplification effect is much larger and damage is more severe near the slope surface. Figure 10b shows that the M PGA variation rule of the bedded slope under vertical seismic force is similar to Figure 10a, and the M PGA reaches the maximum at the slope crest (point A9), approximately 1.19. Moreover, the direction of ground motion has an impact on the dynamic magnification effect of the slope. Figure 10 shows that M PGAx under horizontal seismic loading is about 1.10-1.20 times larger than M PGAz under vertical seismic loading. For example, the M PGAx and M PGAz of A4 are 1.36 and 1.12, respectively. The horizontal seismic force has a greater effect on the dynamic amplification and deformation of the slope. When waves propagate through joints, refraction and reflection phenomena occur on the joint surfaces, which impacts the seismic energy propagation within the slope.

Modal Analysis
Modal analysis has been a common method for the dynamic analysis of structures in the frequency domain [33][34][35][36]. The modal analysis of slopes can clarify the natural frequencies and corresponding vibration modes. The dynamic control equation of the modal analysis is as follows [40] Abaqus/frequency is adopted as the solver for the modal analysis. The modal analysis results are shown in Figures 11 and 12. Figure 11 shows that the first four natural frequencies are 2.53 Hz, 7.21 Hz, 10.64 Hz, and 17.47 Hz, respectively, and the natural frequencies of the bedded slope increase gradually with the order number of the modal mode. The first-order modal shows that the U at the slope crest is significantly greater than other areas (Figure 12a). The U increases with the elevation from the bottom to the slope crest, suggesting that the slope has a typical elevation amplification effect. The U near the slope surface is significantly greater, showing that the slope has a surface amplification effect. The U max appears in the top area, which indicates that the slope crest is more prone to be unstable under the action of the earthquake, which is consistent with the analysis results, using PGA (Figures 8-10). Moreover, the second-order mode shows that U max is mainly concentrated at the top of the slope (Figure 12b), which indicates that the surface area is prone to large deformation. Figure 12c,d show that the U max of the third-and fourth-order modes are mainly concentrated in the local area of the slope crest. The low-order mode suggests that the low-order natural frequency (2.53 Hz) mainly induces the overall sliding deformation of the slope crest area, while the high-order natural frequency mainly induces the local deformation of the slope crest.   According to the modal analysis, the dynamic failure mechanism of the slope can be summarized as follows: under the action of high-order natural frequency, local failure deformation appears at the surface slope that is above the topmost structural plane (potential slip surface); under the action of low-order natural frequency, the surface slope appears as overall shear sliding failure along the potential slip surface; with the increase in seismic load, under the action of higher and lower natural frequencies, the slipping mass gradually expands downward, and further shear sliding failure continues to occur along the lower structural plane.

Dynamic Response of the Slope Using Fourier Spectrum Analysis
To study the amplification effect of the slope, taking the horizontal seismic load (0.1 g) as examples, the Fourier spectra of the points are shown in Figures 13 and 14. There are four predominant frequencies (f 1 , f 2 , f 3 , and f 4 ) in Figures 13 and 14, including the f 1 (2.50 Hz), the f 2 (7.30 Hz), the f 3 (10.50 Hz), and the f 4 (17.30 Hz). The change of topographic and geological conditions directly leads to the change of the magnification effect of the slope [23,24]. Figures 13 and 14 show that the Fourier spectrum amplitude of the first three natural frequencies is much larger, and the PFSA (peak Fourier spectrum amplitude) of the f 1 , f 2 , and f 3 increases gradually with the slope elevation on the whole, which indicates that f 1 , f 2 , and f 3 are closely related to the magnification effect of the bedded slope. In combination with the results of the modal analysis of the slope (Figure 12), the f 1 , f 2 , and f 3 of the slope have an important impact on the modal characteristic of the slope. In particular, the PFSA of the f 1 is the largest of all, and increases rapidly with the slope elevation, indicating that f 1 has a controlling effect on the deformation of the bedded slope. Moreover, when the frequency is >20 Hz, the Fourier spectrum amplitude has no obvious change with the increasing elevation, suggesting that the amplification effect is not obvious in the frequency component range. Therefore, certain frequency components (f 1 , f 2 , and f 3 ) have magnification effects on the dynamic response of the slope, and their PFSA was selected to clarify the dynamic response of the slope during earthquakes. To investigate the dynamic response of the slope, taking the horizontal seismic load (0.1 g) as an example, the PFSA of f 1 , f 2 , and f 3 is shown in Figure 15. The PFSA has an increase with the elevation of the slope, and reaches the maximum at the top slope, suggesting that the slope elevation has an impact on the seismic response of the slope. Figure 15a shows that the PFSA increases linearly overall firstly, but a rapid increase in the PFSA can be identified when the h/H is > 0.62 on the slope surface and the h/H is > 0.82 inside the slope, which suggests that an obvious amplification effect can be found in the surface slope that is above the topmost structural plane, and the difference of magnification effect between the surface slope and the lower rock mass, resulting in the different dynamic response of the slip mass and slip bed. The phenomenon indicates that the structural surface has a magnification effect on the dynamic response of the slope, due to the following reason: when the seismic waves pass through structural planes, the reflection and refraction of waves occur, which results in the wave superposition in some frequency bands, and a great difference of physical properties between structural surfaces and rock masses directly leads to the redistribution of energy [41], directly magnifying the dynamic response of slopes. Figure 15 also shows that the PFSA reaches the maximum at the slope crest, indicating that the dynamic amplification effect is mainly concentrated in the surface slope, where the mechanical properties are very unstable. Moreover, it can be seen that when the slope has multiple structural planes, the topmost structural plane has the most obvious magnification effect on the slope. When the surface slope begins to be damaged, the topmost structural plane is a potential sliding surface. With the increase in seismic load, after the formation of the primary slip mass (surface slope), the potential sliding surface continues to expand in the lower structural plane, leading to further promoting the formation of the larger-scale sliding body.

Analysis of Dynamic Failure Mechanism of the Slope
Weak structural planes have an important effect on the seismic failure evolution process of the bedded slope. Taking the WE wave (0.1 g) as an example, the contour of the shear strain increment of the slope under different earthquake intensities is shown in Figure 16. The shear strain increment in the rock mass is small, while that in the weak structural planes is much large (Figure 16), which indicates that the weak structural planes are the potential sliding surfaces. According to the results of the numerical simulations (Figure 16), the initiation and dynamic failure evolution processes of the landslide are investigated by analyzing the dynamic change characteristics of the shear strain increment of the weak structural surfaces: When the earthquake intensity is 0.1 g, the shear strain increment of the topmost structural surface is relatively large, while that of the lower structural surface is relatively small, suggesting that shear failure is initiated to develop at the slope crest (Figure 16a). Ever-larger shear strain increments occur in the structural planes when the earthquake intensity is 0.3 g and expand to the whole topmost structural plane (Figure 16b), indicating that the surface slope begins to slide along the structural plane, which is consistent with the shaking table tests [23,24]. With the increase of seismic loads, under the excitations of 0.4 g and 0.6 g (Figure 16c,d), the shear strain increment continues to increase and extends to the lower structural planes, which indicates that the dynamic failure area of the slope extends to the lower part of the slope body, which is consistent with the experimental results [23,24].
Therefore, structural planes have important influences on the sliding development process and mode of the bedded landslide, which are the potential sliding surfaces. With the increase in seismic force, the sliding surface extends from the upper structural plane to the lower structural plane. The seismic failure mechanism of the bedded slope is as follows: crack first appears in the topmost structural surface that is the primary slip surface, and the primary slip mass (surface slope) gradually forms. Under the action of the earthquake, the cracks gradually spread in the structural plane, and promote the formation of the sliding surface of the topmost structural plane. When the earthquake force continues to reach a certain value, the surface slope produces shear failure along the bedded structural surface. With the increase in seismic force, the sliding body gradually extends to the lower rock mass, and the structural plane of the lower part continues to act as the sliding plane, which promotes the sliding body to slide on a larger scale.

Discussion
The seismic response of rock mass slopes is a multi-domain scientific problem. Investigation on the dynamic response of rock slopes based on the time domain and frequency domain is helpful to reveal the seismic response characteristics and dynamic failure evolution mechanism of rock slopes with complex geological structures from multiple perspectives. The time-domain analysis can well explore the influence of ground motion parameters and topographic and geological factors on the dynamic response characteristics of rock slopes. The frequency-domain analysis can reveal the relationship between seismic wave frequency, the natural frequency of rock slopes, and dynamic response characteristics of slopes, which can be further verified and supplemented with the results of a time-domain analysis. In addition, the previous frequency analysis content does not fully consider the modal analysis results; the modal analysis can directly reflect the natural frequency and modal characteristics of rock slope. In this work, the frequency research content is extended to the spectral characteristics analysis of waves and modal analysis of slopes, and the mutual verification and supplementation of the analysis results of the two aspects can be fully considered to further improve the frequency domain research content. In this work, the dynamic response characteristics of the layered slope are systematically studied, especially from the perspective of modal analysis and Fourier spectrum analysis, and the influence of seismic wave frequency on the dynamic response characteristics of the slope is revealed. The multi-domain characteristics of the rock slope are fully considered in this work, which can provide a new idea for seismic reinforcement analysis of the actual slopes. However, this work also has some limitations. The time-frequency domain analysis fully considers the characteristics of the time frequency and amplitude of ground motion. The characteristics of the seismic dynamic response of the rock slope with complex geological structures still need to be studied in the time-frequency domain.

Conclusions
According to the dynamic analyses of the bedded slope, the seismic response characteristics of the slope were investigated. Some conclusions can be drawn:

1.
Structural planes have a significant effect on the wave propagation characteristics through the slope. Irregular wave propagation characteristics can be identified through structural planes. To analyze the MPGA, the slope has obvious elevation and slope surface amplification effects. Structural planes have an impact on the slope dynamic response. The MPGA increases gradually with the slope elevation and reaches the maximum at the slope crest. The MPGA shows obvious nonlinear variation characteristics in the slope; in particular, a rapid increase of MPGA above the topmost structural plane can be found. The seismic direction of ground motion is closely related to the magnification effect of the slope. The MPGA under horizontal seismic load is approximately 1.10-1.20 times as much as the MPGA under vertical seismic force.

2.
According to the frequency domain analysis, the first four natural frequencies of the bedding slope are approximately 2.53 Hz, 7.21 Hz, 10.64 Hz, and 17.47 Hz, respectively. The Fourier spectrum and modal analyses show that the slope owns the elevation, slope surface, and structural plane amplification effects. The natural frequency is closely related to the dynamic response characteristics of the slope. The low-order and high-order natural frequencies mainly cause the overall sliding deformation and the local deformation of the top slope.

3.
Structural planes have an important influence on the sliding development process and mode. With the increase in seismic force, the sliding surface extends from the upper structural plane to the lower structural plane. The dynamic failure mechanism of the bedding slope is as follows: a crack first appears in the topmost structural plane, and the primary slip mass (surface slope) gradually forms. Cracks gradually extend in structural planes and promote the formation of the sliding surface; then, the surface slope produces shear failure along the topmost structural plane. With the increase in seismic force, the sliding body gradually extends to the lower slope body, and the structural plane of the lower part continues to form the sliding plane, which promotes the sliding body to slide on a larger scale. Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data used to support the findings of this study are available from the corresponding author upon request.