Numerical Modeling of 3D Slopes with Weak Zones by Random Field and Finite Elements

: This work investigates an analysis method for the stability of a three-dimensional (3D) slope with weak zones considering spatial variability on the basis of two-phase random media and the ﬁnite element method. By controlling the volume fractions of rock and weak zones, two-phase random media are incorporated into the 3D slope model to simulate the random distribution of rock and weak zones. Then, a rotation of a Gaussian random ﬁeld is performed to account for the inclination of the weak zones. The validity of the proposed model for use in the analysis of the stability of 3D slopes with weak zones was veriﬁed by comparing it to existing results and analytical solutions. The failure mechanism of the slope is considered by examining the plastic failure zone at incipient slope failure. The safety factor is sensitive to the inclination of the weak zones, but it is predictable. Parametric studies on the inclination of the layer of weak zones demonstrate that when the rotation angle of the weak zones is approximately parallel to the slope inclination, the slope is prone to slippage along the weak zones, resulting in a signiﬁcant reduction in the safety factor. The ﬁndings of this research can serve as the foundation for further research on the stability of slopes with weak zones.


Introduction
The prediction and control of landslides caused by the instability of slopes with weak zones is of great practical significance in the risk assessment of earth disasters. Slope stability is closely related to the inclination, position and shape of the weak zones inside the slopes. Due to the strong degree of randomness in the spatial distribution of the weak zones within the slope, these factors are, to a certain extent, uncertain, so it is difficult to quantitatively evaluate the stability and failure consequences of the slope [1]. Stability analysis of slopes with weak zones is an essential element of safety design in practice.
The stability of a slope will be strongly affected by the distribution of its internal weak zones. For a slope composed of multiple zones, the damage is usually dominated by the spatial distribution of the weak zones [2]. To date, numerous approaches have been proposed to analyze the stability of slopes with weak zones, such as the rigid body limit equilibrium method, the model test method, and the finite element method [3][4][5][6][7][8][9][10][11].
However, natural rock and weak zones exhibit obvious spatial variability [12][13][14][15]. In some geotechnical engineering practices, due to long-term physical, chemical or geological effects, investigations that do not consider the spatial variability of rock and weak zones can lead to unconservative results [16]. The spatial variability of rock and weak zones can cause randomness and uncertainty in the distribution of weak zones inside the slope. Since the rock and weak zones of slopes are the products of a shared sedimentary history, structure and experience of human activities, the spatial distribution of the weak zones at a Finite element strength reduction method Numerical Reflect the influence of the weak zones on the stability of the slope, but cannot reflect the spatial heterogeneity of the weak zones Liu et al., 2018 [33] Two-phase random media Numerical Reflect the characteristics of the rock-soil slope, but limited to a two-dimensional model, and cannot reflect the spatial variability of the slope shape Huang et al., 2020 [38] Discrete element method Numerical Simulate the nonlinear large deformation characteristics, but have low solution efficiency for complex models Wang et al., 2021 [28] Random finite element method Numerical Characterize the stratification of the soil inside the slope, but cannot reflect the slope of the soil-rock mixture.
This research aims to estimate the stability of a 3D corner slope with weak zones, introducing a two-phase random media model to invert the spatial distribution of rock and weak zones. Based on the basic theory of random fields, a mathematical theoretical model for characterizing rock weak zones using two-phase random media and a simulation method for the spatial inclination of weak zones are proposed. In this study, using the USDFILD platform in ABAQUS software, a Fortran program corresponding to the proposed model is developed, and 3D finite element modeling of slopes with weak zones is realized. The two-phase random media model and the finite element method are combined to quantitatively evaluate the safety factor of the slope, offering technical support for stability studies of slopes with weak zones.

Method for Generating 3D Random Fields
The 3D random field is generated using the modified linear estimation method. The basic principle of this method is to generate a stationary Gaussian random field, where the mean is zero, and the unit variance and the spatial correlation length are √ π [12]. Spatially correlated variables can be expressed by a linear equation combining fixed effects and random effects as follows: where t represents the spatial variation of rock material properties at different locations, and fixed effects and random effects are represented by the mean value µ and the residual e, respectively. e is often assumed to be a constant, with deviations from this having a mean value u. The covariance matrix V of the residual e can be descirbed as follows: where σ is the standard deviation of the random variable at each node, and C is the spatial autocorrelation matrix. To avoid generating negative values for the strength parameters of rock and weak zones, the strength parameters of some rock and weak zones usually adopt log-normal distribution. In this case, the random field can be modeled by the following method conversion: where ε is the vector of correlated random variables; µ ln z is the mean of the logarithm of the material parameters; and σ ln z is the standard deviation of the logarithm of the rock property. The modified linear estimation method needs to generate an n-dimensional m-vector attribute field with spatial cross-correlation matrix C. Then, vector attribute field stretch modeling is used to generate a random field with any specified spatial correlation length. The specific steps are as follows: First, the n-dimensional space containing position vector s is discretized into ndimensional grids with unit grid spacing. In addition, the grid nodes are filled using a random vector r with m random numbers following the standard Gaussian distribution.
Second, Cholesky decomposition is performed on the spatial cross-correlation matrix C according to Equation (4).
where L is the lower triangular matrix. Third, calculate the position of y in s according to Equation (5).
where y is the attribute field containing the position vector; s is the random field containing the position vector; J is the direction vector and ε is a translation vector with n components. J and s are independent random variables, and they have specific values. In this study, the components of J and ε are uniformly distributed in the range of [0, π/2] and [0,1], respectively. According to Equation (2), the attribute field y will be rotated and translated to obtain the random field s. Finally, the continuous attribute field of any coordinates can be calculated by Equation (6), as shown below.
where f j (y) and f j i,k are the jth component of the attribute field containing the position vector y and the ith node of element k, respectively; N i (y) is a form function with 2 n nodal elements, which is taken from the library of form functions established by the finite element method.

Method of Generating 3D Random Field
The stabilities of slope projects are strongly affected by the distribution of its internal weak zones. In the destruction of a slope composed of multiple zones, landslides are usually dominated by the spatial distribution of the weak zones. However, due to the constraints of project budget and construction schedules, it is almost impossible to obtain a complete distribution of zones within the slope through on-site surveys. In typical slope projects, only limited survey work is carried out. Therefore, the distribution characteristics of weak zones are only known in the survey operation area, and such information cannot be obtained at other locations.
The spatial distribution of weak zones at the same site and the spatial distribution of the mechanical characteristics of the rock and weak zones should have similar characteristics. In the random field, the spatial correlation of the layer of weak zones in the layer of weak zones between the two-slope rock and weak zones elements can be expressed using the square exponential autocorrelation function R(x, y, z), as shown below: where R (x, y, z) represents the spatial correlation of the material unit; θ x , θ y and θ z are the correlation lengths along x, y and z directions, respectively. In this study, the anisotropic random field is used to characterize the inhomogeneity of the spatial distribution of the weak zones. When the spatial correlation of each unit in the simulated weak zone is relatively high, it can be regarded as the same weak zone. In addition, since the weak zones are usually continuously layered inside the slope, the simulation results obtained using this method are closer to the ground truth.
For a 3D stationary Gaussian random field, a two-phase random media can be reconstructed for the rock and weak zones using the following formula: where B(x, y, z) is the random media of weak zones and rock; and p 0 is the volume fraction of weak zones. Reconstructions of typical weak zones and rock random media are shown in Figure 1, where 0 represents the volume fraction of weak zones, and 1 represents the volume fraction of rock [33].
the weak zone of the slope according to the actual geological survey situation, and can truly reflect the random distribution characteristics of the layer of weak zones in the slope. On the basis of the reconstruction of the typical two-phase random media shown in Figure  1, two typical samples of the two-phase random media are generated according to Equation (8), as shown in Figure 2, where the volume fraction of the weak zones is as assumed to be 20%. The two-phase random media simulation method based on Gaussian random field is adopted in this research. This method can flexibly control the shape and orientation of the weak zone of the slope according to the actual geological survey situation, and can truly reflect the random distribution characteristics of the layer of weak zones in the slope. On the basis of the reconstruction of the typical two-phase random media shown in Figure 1, two typical samples of the two-phase random media are generated according to Equation (8), as shown in Figure 2, where the volume fraction of the weak zones is as assumed to be 20%. The two-phase random media simulation method based on Gaussian random field is adopted in this research. This method can flexibly control the shape and orientation of the weak zone of the slope according to the actual geological survey situation, and can truly reflect the random distribution characteristics of the layer of weak zones in the slope. On the basis of the reconstruction of the typical two-phase random media shown in Figure  1, two typical samples of the two-phase random media are generated according to Equation (8), as shown in Figure 2, where the volume fraction of the weak zones is as assumed to be 20%.

Method for Characterization of Anisotropic Correlation Structure
To make the simulation more similar to the actual slope, a 3D model was established, taking into account rotational anisotropy. The coordinates of the Gaussian random field were rotated around the X-axis in accordance with the given inclination angle. The rotation process can be expressed as in Equation (9).
where x', y' and z' are the corresponding coordinates of x, y and z after rotation of the α angle about the X-axis. The distance after rotation can be obtained using Equations (10)- (12), as shown below. ∆x = ∆x (10) The expression of the squared exponential autocorrelation function of rotational anisotropy can be obtained as: Figure 3 illustrates the rotationally anisotropic after rotation of the α angle about the X-axis, the transversely anisotropic weak zone structure, and the correlation lengths. It can be clearly found from the figure that transverse anisotropy is a special case of rotational anisotropy. In addition, the autocorrelation function can well characterize the weak zones inside the slope.

Method for Characterization of Anisotropic Correlation Structure
To make the simulation more similar to the actual slope, a 3D model was estab taking into account rotational anisotropy. The coordinates of the Gaussian random were rotated around the X-axis in accordance with the given inclination angle. Th tion process can be expressed as in Equation (9).
where x', y' and z' are the corresponding coordinates of x, y and z after rotation o angle about the X-axis. The distance after rotation can be obtained using Equation (12), as shown below.
The expression of the squared exponential autocorrelation function of rotatio isotropy can be obtained as: Figure 3 illustrates the rotationally anisotropic after rotation of the α angle abo X-axis, the transversely anisotropic weak zone structure, and the correlation leng can be clearly found from the figure that transverse anisotropy is a special case o tional anisotropy. In addition, the autocorrelation function can well characterize the zones inside the slope. For rotated anisotropy around the Y-axis, the correlation distance after rotation can be obtained using Equations (14)- (16), as shown below.
The expression of the squared exponential autocorrelation function of rotational anisotropy can be obtained as: Figure 4 illustrates the rotationally anisotropic after rotation of α angle about the Y-axis, the transversely anisotropic weak zone structure, and the corresponding correlation lengths.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 7 of 1 For rotated anisotropy around the Y-axis, the correlation distance after rotation ca be obtained using Equations (14)- (16), as shown below.
The expression of the squared exponential autocorrelation function of rotational an isotropy can be obtained as:  For rotated anisotropy around the Z-axis, the correlation distance after rotation ca be obtained by Equations (18)-(20), as shown below.
The expression of the squared exponential autocorrelation function of rotational an isotropy can be obtained as: Figure 5 illustrates the rotationally anisotropic after rotation of α angle about the Z axis, the transversely anisotropic weak zone structure, and the corresponding correlatio lengths. For rotated anisotropy around the Z-axis, the correlation distance after rotation can be obtained by Equations (18)-(20), as shown below.
The expression of the squared exponential autocorrelation function of rotational anisotropy can be obtained as: Figure 5 illustrates the rotationally anisotropic after rotation of α angle about the Z-axis, the transversely anisotropic weak zone structure, and the corresponding correlation lengths.  On the basis of the above method, it is possible to characterize the slope where the weak zones are rotated by 0-60° around the X, Y, and Z-axes, respectively. Taking the soil volume fraction of the weak zones as 20%, the two-phase random media B (x, y, z) of the slope were generated, and the soil layer of the weak zones with different inclination angles was inverted. Some simulated samples of the slope are shown in Figure 6.

Simulation and Model Verification of Weak Zones in Slopes
Taking into account the high degree of randomness the parameters such as shape inclination angle, and the position of the weak zones, in order to perform verification, we used the Monte-Carlo simulation method to statistically analyze the results, generating a large number of random field samples through the modified linear estimation method The airport was mapped to a finite element grid in order to construct a 3D model of a slope with weak zones. The stability of each slope in 100 examples of Monte-Carlo simulation was calculated. The safety factor was evaluated and compared with the content of the previous research [33] (Liu et al., 2018) for purposes of verification. As shown in Figure  7, the specific steps of combining USDFLD subroutine random field with the software Abaqus were as follows: (1) Use Abaqus to build a 3D finite element slope model, including defining the material parameters of the rocks and weak zones, and setting the grid and boundary conditions of the model in the analysis step. Then output an inp file (a file executed by the Abaqus software). (2) Associate the yield strength in the material parameters of the slope with the userdefined field variable User Defined Field, and set the state variable Depvar at the same time. On the basis of the above method, it is possible to characterize the slope where the weak zones are rotated by 0-60 • around the X, Y, and Z-axes, respectively. Taking the soil volume fraction of the weak zones as 20%, the two-phase random media B (x, y, z) of the slope were generated, and the soil layer of the weak zones with different inclination angles was inverted. Some simulated samples of the slope are shown in Figure 6. On the basis of the above method, it is possible to characterize the slope where the weak zones are rotated by 0-60° around the X, Y, and Z-axes, respectively. Taking the soil volume fraction of the weak zones as 20%, the two-phase random media B (x, y, z) of the slope were generated, and the soil layer of the weak zones with different inclination angles was inverted. Some simulated samples of the slope are shown in Figure 6.

Simulation and Model Verification of Weak Zones in Slopes
Taking into account the high degree of randomness the parameters such as shape, inclination angle, and the position of the weak zones, in order to perform verification, we used the Monte-Carlo simulation method to statistically analyze the results, generating a large number of random field samples through the modified linear estimation method. The airport was mapped to a finite element grid in order to construct a 3D model of a slope with weak zones. The stability of each slope in 100 examples of Monte-Carlo simulation was calculated. The safety factor was evaluated and compared with the content of the previous research [33] (Liu et al., 2018) for purposes of verification. As shown in Figure  7, the specific steps of combining USDFLD subroutine random field with the software Abaqus were as follows: (1) Use Abaqus to build a 3D finite element slope model, including defining the material parameters of the rocks and weak zones, and setting the grid and boundary conditions of the model in the analysis step. Then output an inp file (a file executed by the Abaqus software). (2) Associate the yield strength in the material parameters of the slope with the userdefined field variable User Defined Field, and set the state variable Depvar at the same time.

Simulation and Model Verification of Weak Zones in Slopes
Taking into account the high degree of randomness the parameters such as shape, inclination angle, and the position of the weak zones, in order to perform verification, we used the Monte-Carlo simulation method to statistically analyze the results, generating a large number of random field samples through the modified linear estimation method. The airport was mapped to a finite element grid in order to construct a 3D model of a slope with weak zones. The stability of each slope in 100 examples of Monte-Carlo simulation was calculated. The safety factor was evaluated and compared with the content of the previous research [33] (Liu et al., 2018) for purposes of verification. As shown in Figure 7, the specific steps of combining USDFLD subroutine random field with the software Abaqus were as follows: (1) Use Abaqus to build a 3D finite element slope model, including defining the material parameters of the rocks and weak zones, and setting the grid and boundary conditions of the model in the analysis step. Then output an inp file (a file executed by the Abaqus software). (2) Associate the yield strength in the material parameters of the slope with the userdefined field variable User Defined Field, and set the state variable Depvar at the same time. (3) Call the random field generated by the USDFLD subroutine to replace the yield strength in the original slope model. (4) Carry out the stability calculation of the slope model with weak zones constructed in (3), and obtain the safety factor of the slope model. (5) According to the above steps, perform 100 Monte-Carlo simulations on the established model, and perform statistical analysis on the calculation results. (3) Call the random field generated by the USDFLD subroutine to replace the yield strength in the original slope model. (4) Carry out the stability calculation of the slope model with weak zones constructed in (3), and obtain the safety factor of the slope model.  To verify the rationality of the two-phase random media characterization method, a typical 3D slope was established as shown in Figure 8. The model contains 213,240 C3D8 elements, and two-phase random media was used to characterize the weak zone and rock mass (using the ideal elastoplastic model); the two-dimensional slope model corresponding to the 3D model was in agreement with the content found in the literature, and the boundary conditions of the 3D model were consistent with the two-dimensional slope model found in the literature [33]. The material parameters used in this section are shown in Table 2. To verify the rationality of the two-phase random media characterization method, a typical 3D slope was established as shown in Figure 8. The model contains 213,240 C3D8 elements, and two-phase random media was used to characterize the weak zone and rock mass (using the ideal elastoplastic model); the two-dimensional slope model corresponding to the 3D model was in agreement with the content found in the literature, and the boundary conditions of the 3D model were consistent with the two-dimensional slope model found in the literature [33]. The material parameters used in this section are shown in Table 2.  Note: -means that this parameter has no unit.
On the basis of the parameters shown in Table 2, a Monte-Carlo simulation (100 Figure 8. A typical realization of two-phase random media for a slope with a weak zone fraction equal to 0.5 (the black area represents weak zones, and the gray area represents rock). (a) Schematic diagram of a typical 3D slope with weak zones (3D model used in this research); (b) schematic diagram of a typical two-dimensional slope with weak zones. Note: -means that this parameter has no unit.
On the basis of the parameters shown in Table 2, a Monte-Carlo simulation (100 times) of the model in Figure 8 was performed, and the results obtained for the average safety factor are shown in Figure 9. When the volume fraction of the weak zones was 50%, the safety factor of the 3D model was 1.3% higher than that of the two-dimensional model   [33]. Considering that the number of Monte-Carlo is 100, the deviation of the mean is basically within the acceptable range. Compared with existing results, it can be verified that the 3D slope method representing the weak zones on the basis of two-phase random media is reasonable. Meanwhile, due to the spatial variability of the material parameters of the slope with weak zones, there is no straightforward rule for the failure mode of the slope and the irregular critical slip surface. The results obtained in this study suggest that the characterization of a 3D slope with weak zones on the basis of two-phase random media, combined with the finite element method, is able to automatically locate the critical slip surface where the internal strength of the slope is insufficient to resist the shear stress.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 11 of 19 be verified that the 3D slope method representing the weak zones on the basis of twophase random media is reasonable. Meanwhile, due to the spatial variability of the material parameters of the slope with weak zones, there is no straightforward rule for the failure mode of the slope and the irregular critical slip surface. The results obtained in this study suggest that the characterization of a 3D slope with weak zones on the basis of twophase random media, combined with the finite element method, is able to automatically locate the critical slip surface where the internal strength of the slope is insufficient to resist the shear stress. At the same time, it can be seen from Figure 9 that when the slope contains weak zones, the safety factor is 0.64, and when the slope is composed of rocks, the safety factor is 3.2 (the simulation results are consistent with the analytical solution). The median safety factor of slopes with weak zones and rocks is 1.92, which is greater than the safety factor of 1.14 obtained when the volume fraction of weak zones is 50%, meaning that the failure mechanism of slopes with weak zones may be dominated by weak zones. To study the failure principle of slopes with weak zones, Figure 10 shows the plastic strain diagrams of homogeneous slopes and slopes with weak zones. It can be seen from the figure that the homogeneous slope forms an obvious overall shear failure zone. In the slope with weak zones, due to the existence of the weak zones, a series of small shear failure zones appear in the failure mode of the slope along the laminar tendency. These small plastic strain regions reduce the safety factor of the slope.  At the same time, it can be seen from Figure 9 that when the slope contains weak zones, the safety factor is 0.64, and when the slope is composed of rocks, the safety factor is 3.2 (the simulation results are consistent with the analytical solution). The median safety factor of slopes with weak zones and rocks is 1.92, which is greater than the safety factor of 1.14 obtained when the volume fraction of weak zones is 50%, meaning that the failure mechanism of slopes with weak zones may be dominated by weak zones. To study the failure principle of slopes with weak zones, Figure 10 shows the plastic strain diagrams of homogeneous slopes and slopes with weak zones. It can be seen from the figure that the homogeneous slope forms an obvious overall shear failure zone. In the slope with weak zones, due to the existence of the weak zones, a series of small shear failure zones appear in the failure mode of the slope along the laminar tendency. These small plastic strain regions reduce the safety factor of the slope.

Slope Stability Analysis
Taking the Gangtou tunnel slope as the object of analysis, as shown in Figure 11, the position and shape of each section of the slope are quite different, and the distribution of weak zones inside the slope has spatial variability. The traditional two-dimensional model can only represent the rotation angle of the layer of weak zones of the weak zones in a certain section. However, it cannot reflect the real spatial distribution of the weak zones inside the slope.

Slope Stability Analysis
Taking the Gangtou tunnel slope as the object of analysis, as shown in Figure 11, the position and shape of each section of the slope are quite different, and the distribution of weak zones inside the slope has spatial variability. The traditional two-dimensional model can only represent the rotation angle of the layer of weak zones of the weak zones in a certain section. However, it cannot reflect the real spatial distribution of the weak zones inside the slope. To illustrate the influence of the rotation angle of the layer of weak zones inside the slope on the stability and failure mechanism of the slope, weak zones with a rotation of 0-60° relative to the X-axis, Y-axis and Z-axis were generated in the 3D slope. The model was about 34 m long, 33 m high and 49 m wide. The number of finite element mesh elements was 399,611 and the number of nodes was 73669. The constraint boundaries around the model were normal bearing constraints, and the bottom of the slope was a complete constraint condition as is shown in Figure 12. Both weak zones and weathered flint strips dolomite are elastoplastic materials, and the Mohr-Coulomb yield criterion was adopted. The rock in the slope shown in Figure 11 is weathered flint-striped dolomite, and the soil in the weak zones is residual cohesive soil. The parameters are shown in Table 3.
(a) (b) (c) Figure 11. Typical failure mode of Gangtou tunnel slope engineering (the rock in the slope is weathered flint-striped dolomite, and the soil in the weak zones is residual cohesive soil).
To illustrate the influence of the rotation angle of the layer of weak zones inside the slope on the stability and failure mechanism of the slope, weak zones with a rotation of 0-60 • relative to the X-axis, Y-axis and Z-axis were generated in the 3D slope. The model was about 34 m long, 33 m high and 49 m wide. The number of finite element mesh elements was 399,611 and the number of nodes was 73669. The constraint boundaries around the model were normal bearing constraints, and the bottom of the slope was a complete constraint condition as is shown in Figure 12. Both weak zones and weathered flint strips dolomite are elastoplastic materials, and the Mohr-Coulomb yield criterion was adopted. The rock in the slope shown in Figure 11 is weathered flint-striped dolomite, and the soil in the weak zones is residual cohesive soil. The parameters are shown in Table 3.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 13 of 19 Figure 11. Typical failure mode of Gangtou tunnel slope engineering (the rock in the slope is weathered flint-striped dolomite, and the soil in the weak zones is residual cohesive soil).
To illustrate the influence of the rotation angle of the layer of weak zones inside the slope on the stability and failure mechanism of the slope, weak zones with a rotation of 0-60° relative to the X-axis, Y-axis and Z-axis were generated in the 3D slope. The model was about 34 m long, 33 m high and 49 m wide. The number of finite element mesh elements was 399,611 and the number of nodes was 73669. The constraint boundaries around the model were normal bearing constraints, and the bottom of the slope was a complete constraint condition as is shown in Figure 12. Both weak zones and weathered flint strips dolomite are elastoplastic materials, and the Mohr-Coulomb yield criterion was adopted. The rock in the slope shown in Figure 11 is weathered flint-striped dolomite, and the soil in the weak zones is residual cohesive soil. The parameters are shown in Table 3.

Results and Discussion
The stability calculation of the slope is shown in Figure 13. The rotation angle of the weak zone layer inside the slope is one of the main controlling factors for the safety factor of the slope. In the figure, when the weak zone is rotated 50° around the X-axis and 40° around the Y-axis, the slope and the inclination of the weak zone are approximately parallel, while simultaneously resulting in the safety factor of the slope being greatly reduced. Therefore, when the inclination of the weak zone is approximately balanced with the slope, the safety factor of the slope will decrease, and the slope is more prone to landslides.

Results and Discussion
The stability calculation of the slope is shown in Figure 13. The rotation angle of the weak zone layer inside the slope is one of the main controlling factors for the safety factor of the slope. In the figure, when the weak zone is rotated 50 • around the X-axis and 40 • around the Y-axis, the slope and the inclination of the weak zone are approximately parallel, while simultaneously resulting in the safety factor of the slope being greatly reduced. Therefore, when the inclination of the weak zone is approximately balanced with the slope, the safety factor of the slope will decrease, and the slope is more prone to landslides. The plastic strain diagram of the slope is shown in Figures 14-16. It can be seen from the figures that when the slope is in a limit equilibrium state, yield occurs first at the weak zones, and the deformation of the weak zones develops continuously, resulting in an increase in the plastic zone. Therefore, when the weak zone layer is rotated at an angle of 50° around the X-axis and a rotation angle of 40° around the Y-axis, the safety factor of the slope reaches its lowest value. The figures suggest that the maximum plastic strain zone is likely to be parallel with the slope and the weak structural surfaces, thus decreasing the stability of the slope. The plastic strain diagram of the slope is shown in Figures 14-16. It can be seen from the figures that when the slope is in a limit equilibrium state, yield occurs first at the weak zones, and the deformation of the weak zones develops continuously, resulting in an increase in the plastic zone. Therefore, when the weak zone layer is rotated at an angle of 50 • around the X-axis and a rotation angle of 40 • around the Y-axis, the safety factor of the slope reaches its lowest value. The figures suggest that the maximum plastic strain zone is likely to be parallel with the slope and the weak structural surfaces, thus decreasing the stability of the slope. The plastic strain diagram of the slope is shown in Figures 14-16. It can be seen from the figures that when the slope is in a limit equilibrium state, yield occurs first at the weak zones, and the deformation of the weak zones develops continuously, resulting in an increase in the plastic zone. Therefore, when the weak zone layer is rotated at an angle of 50° around the X-axis and a rotation angle of 40° around the Y-axis, the safety factor of the slope reaches its lowest value. The figures suggest that the maximum plastic strain zone is likely to be parallel with the slope and the weak structural surfaces, thus decreasing the stability of the slope. The plastic strain diagram of the slope is shown in Figures 14-16. It can be seen from the figures that when the slope is in a limit equilibrium state, yield occurs first at the weak zones, and the deformation of the weak zones develops continuously, resulting in an increase in the plastic zone. Therefore, when the weak zone layer is rotated at an angle of 50° around the X-axis and a rotation angle of 40° around the Y-axis, the safety factor of the slope reaches its lowest value. The figures suggest that the maximum plastic strain zone is likely to be parallel with the slope and the weak structural surfaces, thus decreasing the stability of the slope.

Conclusions
Aiming at the uneven spatial distribution of weak zones in slopes, this work investigated a mathematical theoretical model based on two-phase random media to characterize a 3D slope with weak zones. Based on the USDFILD platform in ABAQUS software, a Fortran program corresponding to the proposed model was developed, and the 3D finite element modeling of the slope with weak zones was performed. By setting the coordinate rotation angle of the anisotropic correlation structure, the inclination angle of the weak zones in the slope model was controlled. The validity of the proposed model for the analysis of the stability of 3D slopes with weak zones was verified on the basis of existing research results and analytical solutions. At the same time, the calculation results show that the above-mentioned characterization method combined with the finite element method is able to realize the automatic retrieval of the internal slip surface of the slope. By inverting the failure mode of the slope with weak zones, it can be found that the stabilities of slopes with weak zones re dominated by those weak zones.
Taking the Gangtou tunnel slope as the object of study, a stability analysis of a 3D slope with weak zones was carried out, and the influence of the inclination angle of the weak zones on the slope stability was analyzed using the 3D finite element model. The results show that the weak zones weaken the strength of the slope, leading to a decrease in safety factor. When the inclination angle of the weak zones was approximately parallel to the slope angle, the slope with the weak zones was more prone to instability. These

Conclusions
Aiming at the uneven spatial distribution of weak zones in slopes, this work investigated a mathematical theoretical model based on two-phase random media to characterize a 3D slope with weak zones. Based on the USDFILD platform in ABAQUS software, a Fortran program corresponding to the proposed model was developed, and the 3D finite element modeling of the slope with weak zones was performed. By setting the coordinate rotation angle of the anisotropic correlation structure, the inclination angle of the weak zones in the slope model was controlled. The validity of the proposed model for the analysis of the stability of 3D slopes with weak zones was verified on the basis of existing research results and analytical solutions. At the same time, the calculation results show that the above-mentioned characterization method combined with the finite element method is able to realize the automatic retrieval of the internal slip surface of the slope. By inverting the failure mode of the slope with weak zones, it can be found that the stabilities of slopes with weak zones re dominated by those weak zones.
Taking the Gangtou tunnel slope as the object of study, a stability analysis of a 3D slope with weak zones was carried out, and the influence of the inclination angle of the weak zones on the slope stability was analyzed using the 3D finite element model. The results show that the weak zones weaken the strength of the slope, leading to a decrease in safety factor. When the inclination angle of the weak zones was approximately parallel to the slope angle, the slope with the weak zones was more prone to instability. These research results provide an important reference for the further quantitative evaluation of landslide susceptibility in the engineering of the slope of the Gangtou tunnel and the selection of excavation slope angles.