Abstract

The instability of high rock slopes with tunnel structures increases under the effects of earthquakes and rainfall. Limited studies have been conducted on the long-term effect analysis of high rock slopes with tunnels. Field electrical measurements on rock slopes have been conducted to study the response of high rock slopes with large cross-section tunnels during rainfall and seismic loading. The material threshold was determined using the resistivity probability density statistical method, and a three-dimensional geological visualization model was obtained. A three-dimensional mesh reconstruction method consisting of material segmentation, cluster filtering, mesh generation, and material attribute mapping modules was proposed. The three-dimensional model of high-cutting rock slope was reconstructed, and the three-dimensional finite element model of the rock slope with tunnel was obtained. The response analysis of the slope was performed during rainfall and seismic loading by employing the El-Centro wave. The results demonstrated that plastic strain was generated in the fracture zone of the upper part of the tunnel during rainfall, and the strain value was 8.135 × 10−4. In the subsequent stages, it had a tendency to expand into the tunnel under continuous rainfall. The maximum principal strain was 3.324 × 10−2 under the effect of a strong earthquake. A large strain was generated in the fracture zone at the upper part of the tunnel. It had a tendency to expand under long-term load, which significantly affected the safety of the tunnel.

1. Introduction

The continuous increase in highway construction has led to the development of tunnels through high slopes with existing tunnel structures owing to the limitation of land use [1]. In a few cases, the tunnel is initially excavated, and the stress in the rock mass is allowed to gradually reach its equilibrium state untill the deformation of the rock mass becomes constant. The surface rock mass is excavated to form a slope. However, the slope stability and excavation response in such cases are different from those in cases without tunnel excavation in the rock mass because they affect the stress and deformation of the rock mass. The interaction between the slope and tunnel should be considered [2]. The presence of the tunnel might cause the collapse or instability of the excavated slope. Additionally, the excavation and disturbance of the slope affect the stability of the tunnel [3]. Cracks develop on the slope due to the construction disturbance, particularly on high slopes. These cracks expand during rainfall and earthquakes, which damage the tunnel and slope. For example, extensive cracking occurred at the top of the slope, and landslides occurred in the Nanshan Tunnel of the Wankai Highway in Chongqing, China, due to disturbances caused during the excavation of the roadbed and frequent precipitation [4].

The stability of rock slopes and tunnels is mainly affected by factors such as rainfall [5], earthquakes [6, 7], geological structures (such as folds, discontinuities, and faults) [8, 9]. The unloading and blasting performed during the excavation process directly reduce the integrity of the rock mass, which leads to fractures and cracks in the weak structural plane [10]. The weak structural surfaces significantly affect the stabilization of rock slopes [11]. Fractures and cracks expand during rainfall and earthquakes, and directly trigger the occurrence of landslides [12]. The tunnel structure and surrounding rock are vulnerable to damage in fracture zones under strong ground motion [13]. The tunnel-fault system can be divided into five main stages, strain localization, fracture initiation, crack acceleration, spontaneous crack growth and stabilization, which is a nonlinear failure process [14].

Several studies have been conducted on the stability of tunnels and slopes by various methods. Jiang et al. [15] proposed a modified BUS approach for the probabilistic back analysis of soil parameters and reliability updating of slopes in spatially variable soils, which facilitates the identification of the causes of slope failures and an understanding of the performance of in-service slopes. Shu et al. [16] analyzed the stability of 3D slopes characterized by spatially variable undrained shear strength based on the plastic limit analysis upper-bound theorem. Vo-Minh et al. [17] presented an isogeometric analysis (IGA) based on Bezier extraction to investigate the stability of a circular tunnel in cohesive-frictional soils subjected to uniform surcharge loading. Rasmussen et al. [18] developed the classic lattice spring model combined with the synthetic rock mass technique and Barton–Bandis joint constitutive model to analyze the probabilistic stability of rock tunnels in low in situ stress environments. Moreover, several studies have been conducted on the stability of high slopes with tunnel structure.

Numerical simulation, model testing, and field monitoring are mainly used to study the interaction between the tunnel and slope under excavation and earthquake [19]. He et al. [20] studied the effect of existing tunnels on the construction of intersecting tunnels in a shallow slope burial context through 3D numerical analysis and field monitoring. Causse et al. [21] focused on the mechanism of slope instability and damages that affect the structure of tunnels parallel to the slope. Kaya et al. [22] performed kinematic, limited equilibrium, and numerical stability analyses to investigate the mechanism of a failure of rock slope with a tunnel. Sun et al. [23] conducted a 1 : 25 large-scale shaking table model test on the bias loess tunnel, and the seismic responses and failure modes of the loess slope and bias tunnel were analyzed. The interaction between the tunnel structure and loess slope was studied. However, the previous studies conducted that the effect of weak structural surfaces on slopes with tunnels are inadequate, and deficiencies were observed in the analysis of the long-term effects of structures under the effects of rainfall and earthquakes.

Field drilling and geophysical exploration techniques are common methods that are used to determine the weak structural surfaces of rock slopes. Geological structures, including weak structural surfaces, can be identified through highly efficient geophysical exploration techniques [24]. In this study, field electrical measurements were conducted to determine the geological structures of the slope with tunnel structure and to study the response of a high-cutting slope with an existing large cross-section tunnel during rainfall and under seismic load. A three-dimensional mesh reconstruction method was proposed to develop a three-dimensional finite element model of the slope with tunnel structure. The long-term effect of the structure was analyzed by simulating rainfall and earthquakes.

2. Three-Dimensional Mesh Reconstruction

2.1. Material Segmentation

The three-dimensional mesh reconstruction method consisted of material segmentation, cluster filtering, mesh generation, and material attribute mapping. Material segmentation is a material processing method that is widely used in three-dimensional mesh model reconstruction. Material segmentation is used to classify the voxels in the three-dimensional visualization model, and each voxel corresponds to a material. Material segmentation can be divided into three types: threshold, edge, and region segmentation. Threshold segmentation is used to obtain the segmentation threshold through various algorithms, which can be considered the boundary of materials. Different material segmentation bases and segmentation methods result in different segmentation effects. The threshold (T) was determined using Otsu’s method [25]. The discrete data points in the three-dimensional visualization model F(i, j, k) were divided into different value range segments with different materials G(i, j, k).where i, j, k are the coordinate values in the three-dimensional visualization model, F(i, j, k) is the discrete data points in the three-dimensional visualization model, and G(i, j, k) is the value of data points after material segmentation.

2.2. Cluster Filtering

A cumulative marking method was proposed for a discrete data lattice to quickly and accurately determine the aggregation clusters in the three-dimensional data lattice. While accessing each data point, the marking of the current data point was determined, and the data points along the three axes, with the marking states of the previous adjacent points were combined. The data points were grouped according to the marking state to perform the identification of the scale of the data points in the aggregation cluster.

The position relationship of the data point (i, j, k) with six adjacent data points is shown in Figure 1. The six adjacent points were divided into three front and back points. The values G(i, j, k), G(i − 1, j, k), G(i, j − 1, k), and G(i, j, k − 1) of the current data point and the three front neighbors were recorded as GS, G1, G2, and G3, respectively.

The data lattice was sequentially marked by scanning along the positive directions of i, j, and k. If the data point of G(i, j, k) = 1 is the target to be marked first, the data point of G(i, j, k) = 0 was not marked temporarily. The cluster data points were combined, and eight previous states of the current element and the integration of the cluster are listed in Table 1, where m is the number of data points in the original cluster, and m1, m2, and m3 are the number of data points in the original cluster where the three front adjacent points (i − 1, j, k), (i, j − 1, k), and (i, j, k − 1) are located, respectively.

The abovementioned cumulative mark scanning was performed, and the data point with G(i, j, k) = 0 was considered the target to be marked and was scanned again. The data points were assigned to different clusters according to their previous states. The cluster with data points m less than that of a certain value was defined as an outlier cluster, and the material attribute value G(i, j, k) of the outlier cluster was converted, i.e., the cluster filtering was performed.

2.3. Mesh Generation

A three-dimensional mesh model was established that formed a mapping relationship with the data dot matrix. The three-dimensional mesh model consisted of eight node hexahedral elements, and each data point corresponded to a unit. The data point was located at the center of the unit, and the size of the unit along the three axes was similar to that of the three-axis spacing of the data point. The nodes of the grid model were numbered, as shown in Figure 2.

The coordinates of the Nth node of the model werewhere N is the number of node. x, y, and z are the coordinates of the Nth node in three directions. U, V, and W are the number of triaxial discrete data points. α, β, and λ are the values of the spacing of discrete data points along three axes. The 8 nodes in the cell were arranged counterclockwise and the number of the 8 nodes (N) in the nth element was determined using the following formula:where n is the number of element and r is the number of node.

2.4. Material Attribute Mapping

The material properties of the nth element in the model were determined using following formula:

The algorithms were developed according to the abovementioned steps to realize 3D mesh reconstruction. The 3D mesh reconstruction method can control the cell size of the 3D mesh model and improve the calculation efficiency to ensure the accuracy of the analysis.

3. Field Electrical Measurement and 3D Mesh Reconstruction of Slope with Tunnel

3.1. Project Overview

An interchange overpass was built on the Xiang’an Airport Highway, Xiamen, wherein a ramp crossed the existing Xiasha tunnel. The Xiasha tunnel was a separated shallow buried tunnel. The ramp was mainly excavated, and the maximum excavation depth was 19 m, as shown in Figure 3.

The results of geological drilling demonstrated the presence of gravelly and strongly weathered granite, fragmentary and strongly weathered granite, and moderately granite from top to bottom, as shown in Figure 4. The rock mass was disturbed, and the cracks expanded during the process of the rock slope excavation.

3.2. Field Test
3.2.1. Test Equipment

The electrical measurement system consisted of a dataset acquisition instrument, battery, booster, cable, and electrodes, as shown in Figure 5. The current and voltage were measured by the electrodes, and the signal was transmitted by the cable. The data acquisition instrument received signals and acquired data. The Wenner method was adopted for this high-density electrical exploration owing to the exploration depth and surface conditions. The electrode spacing was 2 m, and the power supply voltage was 480 V.

3.2.2. Survey Line Layout and Data Acquisition

Two survey lines were arranged to perform electrical measurements along the tunnel direction. Survey line 1 was arranged in tunnel 1, whereas survey line 2 was arranged in tunnel 2. Additionally, 120 electrodes were arranged on each survey line with a spacing of 2 m and a total length of 240 m. The maximum detection depth was 35 m. The layout of the survey line is shown in Figure 6.

3.2.3. Analysis of Test Results

The apparent resistivity profile obtained through field testing was inversely calculated in combination with the least square method, and the resistivity profiles of two survey lines were obtained, as shown in Figure 7. In the case of line 1, a low value of resistivity (Resistivity = 0–100000 Ω·m) was measured on the surface of the slope area (Distance = 60–240 m), with a thickness of approximately 5–10 m. The results were compared with the drilling data, and it was preliminarily determined that strongly weathered granite was present on the surface of the slope. The resistivity gradually increased with an increase in depth, and the resistivity was greater than 100000 Ω·m. This was preliminarily determined to be moderately weathered granite. The slope was in a state of semiexcavation. Therefore, the area (Distance = 0–60 m) was a temporary construction road with large compactness and high resistance in survey line 2. A strongly weathered granite layer with a thickness of approximately 5–10 m was present on the surface, and the resistivity was in the range of 0–100000 Ω·m. Additionally, the deeper layer was a moderately weathered granite layer. The resistivity increased due to the effect of the buried tunnel.

The resistivity was logarithmically converted, and the Gaussian distribution model was introduced, as shown in Figure 8. The resistivity exhibited an approximate normal distribution, and the data were processed according to the law of Gaussian distribution. The distribution exhibited multipeak characteristics, and the thresholds were 10000 and 100000 Ω·m. Combined with the geological survey data, it was determined that the gravelly strongly weathered granite had a resistivity of less than 10000 Ω·m, the fragment strongly weathered granite had a resistivity in the range of 10000–100000 Ω·m, and the moderately weathered granite had a resistivity greater than 100000 Ω·m.

The slope materials were determined in combination with the mathematical statistics method and the material profile of each survey line was obtained, as shown in Figure 9. According to the frequency statistical results, the materials were divided into two categories: strongly and moderately weathered granite. Compared with that of the information from BLK2 drilling points passed by survey line 1, the thickness of strongly weathered granite was approximately 5-6 m. This was in good agreement with the value (5.8 m) obtained from the drilling point. The deeper part was moderately weathered granite, which was further divided based on different joints and fissure states.

4. Long-Term Effect Analysis of the High Slope with Tunnel Structure

4.1. Construction of a Finite Element Model

The inversion was performed using the least square method after format conversion, noise removal, smoothing, and interpolation. The three-dimensional spatial distribution of resistivity was obtained based on the two-dimensional electrical profiles spread in space through block Kriging interpolation. The resistivity of rock is dependent upon various factors such as rock properties, joints and fissures, and moisture content. Different types of rock have significant differences in electrical conductivity, while the resistivity of the same material is concentrated in a certain range. The three-dimensional resistivity model was combined with mathematical statistical methods such as frequency analysis and cluster analysis and divided into several materials with different conductive properties to generate a three-dimensional geological model.

Combined with the 3D geological model, two 3D grid models were obtained through the three-dimensional mesh reconstruction method, which were affected by rainfall and earthquakes, respectively. The 3D grid model was composed of eight-node hexahedral units, and the size of the unit along the three axes was similar, which was 1 m × 1 m × 1 m. The nodes of the mesh model were numbered, as shown in Figure 10. The length and width of the model were 130 and 45 m, respectively. The height above and below the surface was 31 and 46 m, respectively. The surface layer was strongly weathered granite, and the lower part was moderately weathered granite. The model affected by rainfall had 269,747 elements, while the model affected by the earthquake had 292,900 elements. For both models affected by rainfall and earthquake, horizontal displacement was constrained on lateral surfaces, while horizontal and vertical displacement was constrained on the bottom surface. To reduce edge effects on the model affected by the earthquake, infinite elements were adopted on the boundary of the model. Combined with the geological exploration data, the mapping relationship between grid elements and material properties was established, and different material parameters such as elastic modulus and internal friction angle were assigned to the materials defined by the three-dimensional models.

4.2. Material Parameters

Weathered granite is a rock material with weak nonlinear properties, and the Mohr–Coulomb criterion can fulfill the analysis accuracy, whose parameter is simple and has a physical meaning. Weathered granite was sampled by drilling on-site. The unit weight (γ), cohesion (c′), internal friction angle (φ′), Young’s modulus (E), and Poisson’s ratio () were determined by mechanical experiments. The material parameters are shown in Table 2.

4.3. Load and Cases

Rainfall and seismic force were selected as loads, and the cases are shown in Table 3. The change of strongly weathered rock water content was performed for case 1 under a gravity load to simulate rainfall. Based on the relationship between strongly weathered rock mechanical parameters and water content obtained by relevant literature and mechanical property tests [26], the reduction coefficient of 0.7 was determined when the soil reaches an approximately saturated state. And the stress-strain state of the slope and tunnel rock was observed. A dynamic calculation was performed for cases 2–6 to simulate an earthquake. A horizontal ground motion load was exerted at the bottom of the model, and the El-Centro wave was employed. The seismic time interval was 0.02 s and the duration was 30 s, as shown in Figure 11. The dynamic responses of the slope and tunnel, including acceleration, stress, strain, and earth pressure, were analyzed.

4.4. Responses of Slope and Tunnel under Rainfall

The stress and strain distribution of the slope and tunnel, considering the reduction of soil parameters to simulate rainfall, is shown in Figures 12(a)12(c). A stress concentration area was observed in the tunnel area, and the maximum Mises stress value was 3.576 MPa and the maximum principal strain was 3.268 × 10−5. The plastic strain appeared in the fracture zone in the upper part of the tunnel, and the strain was 8.135 × 10−4. A tendency to expand the tunnel was observed under continuous rainfall. The displacement of the slope and tunnel was small, and the maximum displacement was 9.579 × 10−6 m, as shown in Figure 12(d).

4.5. Responses of Slope and Tunnel under Earthquake

The acceleration response of the slope and tunnel is shown in Figure 13. The maximum accelerations at VI and VIII degrees were 0.57 and 2.97 m/s2, respectively, which had an amplification effect on the slope.

The infinite element boundary was adopted and the seismic load was continuously applied at the bottom of the model. The stress-strain state of the slope and tunnel is shown in Table 4. The maximum Mises stress value of the slope and tunnel under different seismic intensities at peak seismic wave acceleration (t = 5 s) is shown in Figure 14. The maximum principal strain under different seismic intensities is shown in Figure 15. The maximum Mises stress initially increased under the same seismic intensity and then subsequently decreased. The stress rapidly increased with an increase in seismic intensity. Stress concentration areas were observed in the tunnel area. The maximum Mises stress values at t = 0, 5, and 10 s were 1.694, 11.640, and 4.367 MPa when the seismic intensity was IX. A large strain was observed in the upper fracture zone of the tunnel with a tendency to expand into the tunnel under the long-term load, which significantly affected the safety of the tunnel.

Compared with other monitoring points, the strain at D3 is higher due to the existence of a fracture zone, which is about 10.76 × 10−4. Four settlement monitoring points were arranged at the vault of the tunnel on-site to obtain the change in the tunnel settlement during excavation, as shown in Figure 16. The monitored settlement data along the tunnel is 0, −2.3 mm, −4.2 mm, and −3.5 mm. The trend is consistent with the strain calculation results, and the settlement at D3 is higher due to the existence of a fracture zone.

5. Discussion

5.1. Responses of the High Slope with Tunnel Structure under the Combination of Rainfall and Earthquake

It is inevitable that a tunnel will be close to a noncausative fracture zone during the line selection design of a tunnel. Several studies have been conducted on the seismic dynamic interaction between a noncausative fracture zone and lined tunnel. The stress amplification effect of the lined tunnel increased when the lined tunnel was in the foot wall of the fracture zone [27]. However, the disturbance on the fracture zone in the upper part of the tunnel was not considered, and the distribution of the fracture zone was mostly artificially assumed. The slope with tunnel model was developed using the three-dimensional mesh reconstruction method, and the distribution area of the fracture zone was accurately obtained. Fracture growth was observed with a release of high geo-energies due to the excavation of the high slope [28]. Additionally, the thickness of the overlying soil layer of the tunnel decreased, and the rock in the fracture zone was further disturbed due to rainfall. The degradation of the mechanical parameters of the rock under the effect of rainfall was observed owing to the presence of the fracture zone. This severely affected the ground motion of earthquakes due to the multiple reflections of waves in the fault fracture zone [29, 30]. The concentration of stress and plastic strain appeared in the fracture zone. There was a tendency to expand to the tunnel under the long-term load. It is crucial to study the failure mechanism of a high slope with tunnel structure under rainfall and earthquake.

5.2. Failure Mechanism Analysis of the High Slope with Tunnel Structure

The combined finite-discrete element method (FDEM) can be employed to investigate the failure mechanism and failure process of rock masses and the effect of different factors on the failure mode [31]. The fracture zone and tunnel can be divided into discontinuous elements, and other areas can be set as continuous elements to obtain the three-dimensional numerical model [32]. Combined with the damage evolution of rock under rainfall and earthquake, the damage effects of rainfall and earthquake can be coupled into the constitutive model of rock based on continuous damage mechanics and statistical theory [33]. An elastoplastic constitutive model considering damage deterioration can be established, and the mechanical parameters of continuous elements can be determined [34]. Combined with the quantitative correlation between the characteristic parameters of rock microstructure and macroparameters under rainfall and earthquake, the mechanical parameters of discontinuous elements can be defined. The interaction relationship between discontinuous elements can be determined based on the contact constitutive law, and the boundary conditions of the coupling domain can be determined [35]. Research can be performed on the deformation and failure processes of a high slope with the tunnel structure considering the effects of rainfall and earthquake, and its failure patterns and mechanisms can be determined.

6. Conclusions

(1)A field test was performed on the cutting slope with an existing large cross-section tunnel through electrical measurement. To obtain the 3D geological visual model, a mathematical statistics method was used to determine the material thresholds of strongly weathered and moderately weathered granite, which were 10000 and 100000 Ω·m. The three-dimensional mesh reconstruction method was used to reconstruct the 3D mesh model of the slope with the tunnel for the long-term effect analysis.(2)The rainfall and seismic cases were simulated based on the reconstructed finite element models. For the effect of the rainfall, plastic strain was observed in the fracture zone of the upper part of the tunnel. The strain value was 8.135 × 10−4 with a tendency to expand in the tunnel. For the effect of the earthquake, the maximum principal strain was 3.324 × 10−2. A large strain was observed in the fracture zone at the upper part of the tunnel with a tendency to expand into the tunnel, which significantly affected the safety of the tunnel.

Data Availability

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

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the Natural Science Foundation of Fujian Province (2022J011252), and the Science and Technology Project of Xiamen Construction Bureau (XJK-2021-9).