Abstract
Rainfall seepage changes the mechanical properties of loess masses. Considering fluid-solid coupling, the calculation model of a loess tunnel is established according to the finite element method (FEM). Based on porous media seepage theory and rainfall infiltration depth theory and considering the infiltrated depth of the loess surface for different rainfall intensities over a certain period, the stability of a loess tunnel under different rainfall amounts and loess cover thicknesses is studied using the dynamic finite element static strength reduction method. The results show that under the same rainfall intensity, the safety factor increases with the depth of the tunnel; the safety factor of the loess tunnel with the same loess cover thickness decreases with increasing infiltration depth. The plastic strain is mainly distributed on both sides of the vault and the arch feet. The stability of the loess tunnel is directly related to the loess cover thickness and rainfall seepage.
1. Introduction
Most loess is distributed in arid and semiarid mountainous areas. The groundwater depth is relatively deep, and the loess overburden layer is often in an unsaturated state. Thus, the seepage action is mainly affected by rainfall. During the rainy season, with heavy rainfall, especially under earthquake action, a variety of secondary geological disasters can easily be triggered when structures are influenced by rain seepage. Li et al. [1] considered a series of engineering problems resulting from the complex process of permeation-collapsibility-permeation during construction in collapsible loess areas and conducted permeability tests based on loess collapsibility. Yang and Li [2] studied the stability of the roof of a shallow tunnel and ground subsidence using reliability theory and the limit analysis method. Mei et al. [3] conducted a series of experimental studies on the deformation and shear strength properties of compacted loess. Qian et al. [4] carried out field uplift load tests on 18 straight-sided and 15 belled shafts at three collapsible loess sites under an arid environment on the Loess Plateau in Northwest China. In addition, for tunnel engineering, a large number of engineering accident cases show that tunnel collapses, landslides, etc., are induced by surface rainfall; in particular, in the case of rainfall, the pore ratio, saturation, and permeability coefficient of loess are changed by rain infiltration. In this way, the loess layer of a certain depth rapidly reaches a saturated state, and the mechanical properties of the loess mass are changed. Thus, it is necessary to study the stress state of a loess mass in the case of rainwater seepage.
Due to the unfavorable characteristics of the loess, topography, hydrogeological conditions, climate, and other factors in the western region, the loess tunnel is located in a complex vehicle stress and infiltration field, and the coupling effect on the structure of the loess tunnel is a complex state that changes constantly with time and space. Earthquakes are a highly random natural disaster, and the uncertain and complex loess sedimentation process is affected by tectonic movement and human activities. Considering the action of ground motion, the soil around the lining of a tunnel generates a certain acceleration, and the deformation and movement of the soil react to the infiltrated rainwater around the tunnel section, which drives the movement of the soil and intensifies the action, resulting in the instantaneous action of additional excitation. This process promotes interactions between the pore water of the soil mass and the tunnel lining and the surrounding soil mass, which requires a transition from the single method of unsafe factor analysis to the method of multifactor coupling analysis and from a single field analysis to a study of the complex state of multifield coupling. In engineering calculations, the study of seepage starts from the soil. The seepage theory for a porous continuous medium has been used to analyze soil seepage [5]. Analysis of the seepage of the fractured medium starts from the single fracture and extends to the complex network system model. Yoo [6] used a 3D finite element model (FEM) of the spatial characteristics of the weak zone to study the process of tunnel excavation and conduct an in-depth analysis of the spatial characteristics of the tunnel vault, the side wall 3D displacement of the side wall of the tunnel and the weak zone, and the initial stress state. Yang and Yan [7] used nonlinear criteria based on the upper bound theorem and a novel failure mechanism to study the influence of the seepage force and discussed the possible collapse form of deep buried tunnels in layered soil. Liu [8] proposed the seepage flow law of deformable fractured media under a load according to the closed deformation rule of the structural plane. Anagnostou and Kovári [9] developed a hypothetical premise of a groundwater tunnel with the excavation surface in a stable condition. Broere [10] summarized the limit of pore water pressure for a groundwater tunnel with a stable excavation surface based on experimental testing. Wang et al. [11] simulated metro traffic loading for the Xi’an Metro Line 2 by establishing an analytical model of a subgrade train-track system and a material mechanical model of a loess soil mass under the effects of long-term cyclic loading. Soubra [12] applied the upper bound method of limit analysis theory to calculate the active and passive limit pressures before a pressure shield. Lee et al. [13] used the limit equilibrium theory to consider the stability of a soil body under the action of water. Wang et al. [14] used FLAC3D to carry out a complex three-dimensional numerical simulation and a safety evaluation of the surrounding rock of the spillway tunnel of the Yangqu Hydropower Station and obtained 3D slide arcs with a good shape. Using the Xiamen Subsea Tunnel as an example, Li et al. [15] analyzed the stability of the tunnel under fluid-solid coupling by using the FEM. Using the ADINA simulation structure and the fluid field, Cheng et al. [16] established a calculation model of a curved wall loess tunnel with a superlarge section (LTSLS). Wang et al. [17] tested the relationship between the critical groundwater seepage velocity and the subsurface erosion using self-made instruments. The results of tests and numerical simulations indicate that the pipelines not only change the original seepage field of a loess foundation but also accelerate the underground erosion of the foundation. Lee and Nam [18,19] studied the seepage distribution range of groundwater during the excavation of rocks surrounding a shallow buried circular tunnel. Saito et al. [20] proposed that decreases in the permeability of surrounding rocks lead to changes in the pore water pressure of the rock in surrounding tunnels. Farhadian et al. [21] optimized the analytic equation of groundwater seepage in tunnels. Li et al. [22] studied the deformation law of the rock surrounding a double-arch tunnel under the action of groundwater seepage. Huang and Yang [23] studied and compared variations in the safety factor of a tunnel soil mass and the distribution range of the elastoplastic zone under fluid-solid coupling action with different external environments. Yin et al. [24] analyzed the strain and stress distribution of a rock mass by calculating the safety factor of a cross sea tunnel. Pan et al. [25] considered the action of groundwater seepage and surface runoff and analyzed and discussed the failure mechanism of loess landslides. Cheng et al. [26] studied the influence of safety and stability on a subsea tunnel under fluid-solid coupling and dynamic action using the FEM. Farhadian et al. [21] used practical engineering as a reference to optimize the analytical equation of a tunnel with groundwater seepage. Wu et al. [27] conducted a series of physical tests to simulate rain-induced slope failures, and the results confirmed the hypotheses that the pore water pressure and water content in a loose soil slope change rapidly and that water infiltration into cracks in the slope greatly impacts landslide development. Cheng et al. [28] studied the stability of the seismic and dynamic response of a subsea tunnel under different temperatures and seepage amounts, which is beneficial to tunnel engineering applications. Zheng et al. [29] analyzed and compared the settlement law of the rock/soil mass of a tunnel with different permeability coefficients, and the calculated results were consistent with observations.
Although previous studies have considered the static and dynamic qualitative and quantitative mechanical properties of loess tunnel structures during construction, most studies have only considered a single factor. The number of studies that considers loess tunnel structures under the action of various dynamic response factors or that quantitatively and qualitatively analyzes the stress, strain, displacement, acceleration, and pore water pressure variation rule for loess tunnel structures is relatively small. Based on the above research, the seismic stability of loess tunnels with rainwater seepage is studied using the dynamic finite element strength reduction method, which provides a theoretical reference for practical engineering applications, safety analyses, and design of loess tunnels.
2. Calculation Parameters and Analysis Model
2.1. Boundary Conditions
Many vibration and wave problems [30] related to elastic infinite foundations exist in civil engineering. To improve the computational efficiency, a computational model is established over a certain range, and an artificial boundary is placed at the boundary of the computational model, which allows the infinite distance condition to be converted to the artificial boundary condition, thereby reasonably eliminating the influence of the infinite domain. In recent years, the practicality of this approach has become a concern. The simple boundary, transmission boundary, viscous boundary, uniform boundary, Smith superposition boundary, and viscoelastic boundary have been widely used, and the normal and tangential spring stiffness and damping coefficient of the viscoelastic boundary can be calculated according to the following formula [31]:where is the normal spring stiffness, is the tangential spring stiffness, is the damping coefficient of normal dampers, is the damping coefficient of tangential dampers, is the distance from the wave beginning to the artificial boundary point, is the medium shear modulus, is the P wave velocity of the medium, is the S wave velocity of the medium, is the normal viscoelastic boundary correction coefficient, and is the tangential viscoelastic boundary correction coefficient. In general, = 0.8–1.2, and = 0.35–0.65. In this paper, = 1.0, and = 0.5. and can be calculated according to the following formula:where is the Poisson ratio, refers to the water content, and E is Young’s modulus.
2.2. Stability Analysis Theory
2.2.1. Soil Strength Theory
The Mohr–Coulomb strength theory is obtained by expanding the Coulomb strength formula to a three-dimensional state. In general, in engineering calculations, the linear strength criterion of soil can be expressed by the Mohr–Coulomb yield criterion. The equation can be expressed aswhere is the normal stress on the shear sliding surface, is the shear strength of soil, c is the cohesive strength of soil, is the internal friction angle of soil, and is the internal friction coefficient of soil.
2.2.2. Strength Reduction Method
(1) Modal Analysis. In general, in the dynamic finite element analysis of various projects, the Rayleigh damping matrix model is often used. Assuming that the damping model is related to the loess mass and stiffness matrix of the structure,where C is the damping matrix, M is the loess mass matrix, and K is the stiffness matrix.
According to the orthogonality of the matrix, the relationship between , , and can be expressed as [32]with fixed, the values of and only relate to the natural frequency of the structure. The first-order and second-order frequencies of the structure can be obtained by modal analysis. The frequencies are introduced into equation (7). and can be obtained by combining equations (7) and (8):
From equation (8),
(2) Dynamic Equation and Solution. The matrix differential equation of a structure under earthquake action can be expressed as [33,34]where , , and are the acceleration, velocity, and displacement vectors of a model node, respectively; is the earthquake acceleration; and is the vector of the surface load.
The Newmark numerical integration is applied to equation (10), and we assume thatwhere and are the constant values. At the time of , the differential equation of motion is
Taking , , and (where is the maximum natural vibration period of the model), the Newmark numerical integration is unconditionally stable, and the result is satisfactory.
Inserting equations (11) and (12) into equation (13) yields
From equation (13),
Combining equation (15) with equation (14) yields
From equation (16), the value of can be obtained along with the horizontal and vertical displacement vectors. For the plane problem, the stress calculation expression can be obtained by the displacement method.
(3) Strength Reduction Method. In general, the instability of loess is mainly caused by decreases in the loess strength under the influence of external factors. The strength reduction method [33–36] is used to simulate the structural instability by reducing the shear strength parameters of the loess until the structure reaches the critical failure state. The core concept of the method is that the shear strength index of loess (cohesion force c and internal friction angle ) is divided by a reduction factor . The results are shown in equations (17) and (18). Then, we can obtain the new shear strength index (cohesion force c′ and internal friction angle ) through reduction, which can be used to replace the original larger shear strength index. The calculation and analysis process is repeated until the shear strength parameter reaches the value that makes the structure to reach the critical failure state. Thus, the reduction coefficient corresponding to the limit equilibrium state is the maximum safety factor of the tunnel:
Namely,where is the cohesion force of loess after reduction and is the internal friction angle after reduction.
2.2.3. Yield Criteria
The yield criteria are the conditions used to judge if a point in a material has entered a plastic state from an elastic state and to determine the sign of a point in a material that has entered the plastic stress stage with the continuous development of this deformation state. In the definition of the geotechnical yield condition, the yield of loess is related to the hydrostatic pressure. The Mohr–Coulomb yield criterion is often used for loess, and the equation can be expressed aswhere is the first invariant of the stress tensor, , and is the second invariant of the stress tensor, .
In the Mohr–Coulomb model, the parameters and can be expressed by using c (cohesion force) and (internal friction angle):
2.2.4. Calculation Process of the Finite Element Strength Reduction Method
In an actual tunnel project, structural failure usually depends on the shear strength of the loess. Therefore, with the dynamic stability analysis of the tunnel loess mass, the shear strength parameters of the loess can be continuously reduced. In this paper, the dynamic finite element strength reduction method, which can be used to conduct stability analyses for the tunnel, is proposed. Based on the dynamic analysis model, the data of the model analysis are input into the static analysis model, and the shear strength parameters (cohesion force and internal friction angle) of the loess are reduced. By decreasing the shear strength parameters of the loess mass, the static instability failure of the loess mass is shown until the model does not converge, and the obtained reduction parameter is the safety reserve of the loess mass. The specific implementation steps are as follows: first, a structure analysis model is established, and modal analysis is selected as the analysis type to obtain the natural frequency of the structure. Second, the seismic wave and Rayleigh damping coefficient are input into the dynamic analysis model to obtain the moment of maximum horizontal displacement of the model vertices, along with the horizontal displacement values of each node at the left and right sides of the model. Third, in the static analysis model, static strength reduction is used to obtain the safety factor of the loess mass by taking the horizontal displacement as the initial displacement, applying the horizontal displacement to each node at the left and right sides of the model, and reducing the shear strength parameters of the internal friction angle and cohesion force of the loess mass until the calculation result does not converge due to excessive plastic deformation.
2.3. Percolation Theory
2.3.1. Porous Media Seepage Theory
Seepage mechanics is a branch of fluid mechanics that mainly focuses on the movement forms and laws of fluid in porous or fractured media. Because loess is a heterogeneous, anisotropic porous medium, the flow and flow direction of water and the changes of the water and seepage fields for the loess caused by seepage are closely related to actual engineering problems.
For the seepage problem of loess material, the relationship between seepage energy loss and seepage velocity was first defined by the French hydraulic engineer Darcy in 1956 based on a large number of experimental studies, which revealed that the flow through a certain cross-sectional area for a unit time in porous media is inversely proportional to the seepage distance and proportional to the cross-sectional area of water and the head loss of water pressure. The results can be expressed aswhere is the fluid velocity, k is the permeability coefficient, and J is the hydraulic gradient.
2.3.2. Seepage Theory of Rainfall
In the loess region, due to the deep groundwater level and extensive distribution of unsaturated loess, the water content near the surface is low. Consequently, according to changes in the loess moisture content, the water content distribution during water infiltration under ponding conditions can be divided into the saturation zone; the transition zone, which exhibits an obvious change in the moisture content; the conduction zone, which shows little change in water content; and the wetting zone, in which the moisture content decreases rapidly to the initial value. The distribution range of these areas is directly related to the physical parameters of the loess, and the infiltration rate can be used to express the flow of rainwater into the loess per unit time. Usually, the infiltration rate at any time can be expressed by :where is the water content gradient and is the diffusion rate of unsaturated loess.
Based on the above data, the absolute value of the moisture gradient for any location on earth’s surface is very large, and the corresponding permeability value is also large under the condition of a stable water source in the initial stage of seepage. With the passage of time, the moisture content gradient and the permeability continuously decrease, and the infiltration rate also tends to stabilize. Under stable rainfall seepage, the loess is regarded as an equivalent continuum, and the continuity equation of saturated seepage iswhere is the water density, is the Darcy velocity, is the porosity, and is the source sink term.
For unsaturated seepage flow, the porosity, n in equation (25), can be regarded as the moisture content of the unsaturated loess medium. , where is the saturation and . The saturated Darcy flow rate can be transformed into the Darcy velocity in the unsaturated seepage field. The continuity equation of the unsaturated seepage problem is
Assuming that the Darcy laws of saturated and unsaturated seepage are both available, Darcy’s law in unsaturated seepage can be expressed as
Combining equation (27) with the continuity equation (26) yieldswhere is the saturated permeability tensor and is the ratio of the unsaturated permeability coefficient to the saturation permeability coefficient, which is a function of the saturation or pressure head. In the unsaturated zone, . In the saturation zone, .
, where is the position head and is the pressure head. Considering the functional relationship among , , and , , and is defined as the water holding capacity. Considering that the compressibility of water is very low, can be regarded as a constant. The saturation, , can be replaced by the water content volume, . Differential equations of the saturated and unsaturated seepage can be obtained after calculating equation (29):where the water holding capacity, , is equal to 0 in the saturation zone and is the unit storage capacity. For an unsaturated rock loess mass, the value is 0; for saturated loess, the value is a constant. In many cases, , and the value of is 1 in the saturation zone. In the unsaturated zone, the value of is 0.
In the case of homogeneous isotropy, equation (30) can be expressed aswhere the unsaturated permeability coefficient is , , and is the isotropic permeability coefficient in the saturated state.
2.3.3. Rainfall Infiltration Depth
Due to the collapsible property of loess, the water body rates continue to improve under the action of rain seepage, and the shear strength of the loess declines rapidly, which threatens the safety of the tunnel loess mass, especially under the conditions of continuous rainfall. The changes in the water content in the deeper loess layer with rainfall and rainwater infiltration time show a continuous, gradual increase. In recent years, some achievements have been made in related research. With high-intensity rainfall, continuous long-term rainfall, and a large amount of rainfall infiltration, the slope stability worsens. Rainfall infiltration increases loess saturation, which increases the bulk density of overlying loess. This results in an increased sliding force of the loess mass. The increase in the loess side slope moisture caused by rainwater infiltration and the shear strength parameters of loess and the matrix suction of the loess mass are decreased by the sharp decline in loess matric suction, which decreases the slope stability.
The infiltration depth refers to the vertical distance between the wet and dry boundary layers and the ground after a continuous rainfall process and in the absence of water on the surface. Due to some differences in the relative altitude between two existing points, rainfall can infiltrate in loess, which drives water flow downward.
A point in a general water phase relative to a datum consists of three kinds of energy: gravitational potential energy, pressure energy, and kinetic energy. The selection of the reference surface during infiltration is shown in Figure 1.
Considering the energy state of point A, point A has a gravitational potential energy, :where is the water quality at point A and is the height of point A relative to the datum plane.
The energy, , at point A induced by pressure can be expressed aswhere is the pressure energy, is the pore water pressure at point A, is the water flow velocity at point A, and is the water density at point A.
Therefore, the kinetic energy of point A can be expressed as
The total potential energy of point A can be expressed as
The total energy of point A is called the energy head. The energy head of the point can be obtained by dividing equation (33) by the weight of water at point A:
In general, the infiltration of rainwater can only reach a certain depth. When reaching a certain depth, the potential energy between the base level and the maximum penetration depth is 0. That is, . Because some resistance is involved in the rainwater infiltration process, the energy consumed by the loess resistance during rainwater infiltration needs to be considered. Assuming that the infiltration path is the height of vertical infiltration, the value of the hydraulic slope is . The osmotic resistance can be expressed as
The work of overcoming resistance can be expressed as
With no water infiltration, the value of is 0. Based on the conservation of energy,
Equation (39) is the theoretical calculation formula of the maximum depth of rainwater infiltration. Due to the limit of the infiltration depth, for continuous rainfall, surface flow can form. Therefore, the maximum infiltration depth determines the maximum infiltration capacity.
2.3.4. Differential Equation of Seepage
As shown in Figure 1, under the condition of the plane state, water seepage follows Darcy’s law. The stable seepage of an anisotropic continuum equation iswhere h is the energy head, is the permeability coefficient in the x direction, and is the permeability coefficient in the y direction.
(1) Initial Condition. For saturated and unsaturated seepage fields, the entire seepage region can be divided into saturated and unsaturated regions. The boundary conditions include the head boundary, the known flow boundary, and the overflow boundary. The seepage research area is formed by uniting the saturated and unsaturated zones without considering the free surface flow supply boundary. The free surface flow supply boundary can be ignored.
Initial condition: the energy head is a function of the following coordinates:
Head boundary condition:
Flow boundary condition:
Overflow boundary condition:where and are the normal vectors, with the direction of the vertical cross section in the positive direction. For rainfall special flow boundary conditions, = ; is the outer normal direction cosine of the boundary, is the initial time, is the energy head boundary, is the flow boundary, and is the saturated overflow boundary.
(2) Finite Element Seepage Solution. According to the finite element method, assuming that the entire seepage field is continuously composed of finite elements, the energy head distribution of the seepage field depends on the energy head value of each node. The energy head function of any point in the unit can be expressed as
According to the variation principle and finite element theory, equation (40) can be expressed in matrix form as follows:where k is the seepage matrix, is the node head matrix, and is the node flow matrix.
According to the above three kinds of boundary conditions, the energy head value of each node, h, can be obtained. Based on the above equations, the solution of the entire seepage field can be obtained.
2.4. Calculation Parameters
The lining structure of the tunnel is defined as a composite concrete lining. The constitutive loess mass in the tunnel is defined as a Mohr–Coulomb material. The double-layer concrete lining around the tunnel is defined as an elastic-plastic material. The specific material parameters of the concrete lining and loess are shown in Table 1. The material parameters of loess before and after rainfall over a certain range are compared with the loess material parameters before rainfall (Table 2).
2.5. Analysis Model
The section type of the tunnel is the curved wall type. The span of the curved wall section is 15 m, where Hdis the overlying thickness and the height is 11.25 m. Considering the speed and precision of the numerical simulation process, the length of 5 times the tunnel span is taken on the left and right sides of the tunnel. The length of 5 times the tunnel height is taken at the bottom of the tunnel. When defining the boundary conditions, the artificial viscoelastic boundary is set along the left and right sides of the tunnel. When the property of the viscoelastic boundary is set, based on equations (3) and (4), the normal and tangential spring stiffness and damping coefficient are obtained, respectively. = , = , = , and = . The bottom of the tunnel is provided with a fixed hinge constraint, and the tunnel top is free. The dynamic analysis model is shown in Figure 2. If the viscoelastic boundary between the left and right sides of Figure 2 is removed and changed to a horizontal constraint, the static analysis model of the structure is obtained.
Finite element analysis software was used to establish the structure analysis model of the loess tunnel, and the viscoelastic boundary of the model was defined as a spring element when the corresponding element was defined and the mesh was divided. The surrounding rock and lining structure are divided by a mapping network, and the grid division of the model is shown in Figure 3.
When considering the train load for the model in Figure 3, the tunnel bottom padding, concrete foundation, and three two-dimensional graphic track plate units are defined in turn from the bottom of the tunnel to the top. The fixed points are set at an interval distance as the position of the moving train point load. The network partitioning and a detailed analysis of the model are shown in Figure 4.
Considering earthquake action, the safety of a loess tunnel is the most disadvantaged one under the action of a near-field impulse earthquake. In this paper, only near-field pulse seismic action is considered, and the seismic wave is input along the horizontal direction of the model. When defining the element of the structural model, the material definition and element setting of the loess mass before and after rainfall are needed, respectively, and the loess around the rainfall is uniformly set as a linear elastic Mohr–Coulomb material. The seepage effect of the tunnel is considered for overlying loess thicknesses of 30 m, 60 m, and 80 m and under moderate rain, heavy rain, and hard rain action.
In this paper, it is assumed that the initial moisture content of the loess is less than 15%. The structure of the loess tunnel encountered moderate rain, heavy rain, and hard rain erosion, corresponding to three different rainfall intensities of 30 mm/d, 65 mm/d, and 75 mm/d, respectively. Considering the rainwater infiltration depth, based on the theoretical solution [37] of the reference rainfall infiltration depth, the rainwater infiltration depth was 1 m, 1.5 m, and 2 m, respectively. Considering the seepage situation of a loess tunnel with different overlying thicknesses, the infiltration depth remains unchanged.
3. Seismic Stability
3.1. Seismic Stability under Moderate Rain Action
The rainfall intensity is set as 30 mm/d under moderate rain seepage action. With the rainfall penetrating to 1 m below the surface, the thickness of the overlying loess layer is taken as 30 m, 60 m, or 80 m. The near-field seismic pulse is input along the horizontal direction of the structural analysis model, and the horizontal displacement-time history curves of the model’s right vertices of 818, 596, and 4097 can be obtained from the postprocessing module. The results are shown in Figure 5.
(a)
(b)
(c)
The maximum horizontal displacement values of the tunnel model with different loess cover thicknesses can be obtained at 10.67 s, 11.91 s, and 12.10 s. Then, taking the displacement values of each point in the analysis model with different loess cover thicknesses into each node at the left and right boundaries of the static model, the shear strength parameters of the surrounding loess mass are continuously reduced until the computational results show no convergence. The plastic strain nephograms of the loess tunnel before critical failure under moderate rain seepage action are shown in Figure 6.
(a)
(b)
(c)
From Figure 6, for a loess tunnel with different loess cover thicknesses, the maximum safety factor of a superdeep tunnel is 2.673, and the safety factor is 2.278 in the case of a burial depth of 60 m. However, the minimum safety factor for a shallow burial depth is 2.267. The critical plastic strain is mainly distributed on both sides of the arch and toe of the tunnel and expands with the tunnel depth around the surrounding mass. The maximum strain is located at the arch-toe of the tunnel.
3.2. Seismic Stability under Heavy Rain Action
Considering heavy rain seepage and earthquake action with the near-field pulse, the rainfall intensity is set as 65 mm/d. The rainwater infiltration depth is 1.5 m. The overburden thickness is taken as 30 m, 60 m, or 80 m. The horizontal displacement-time history curves of the model’s top-right corner vertices can be obtained for the different loess overlying thicknesses. The results are shown in Figure 7.
(a)
(b)
(c)
Figure 7 shows the maximum horizontal displacement values of the tunnel model with different loess cover thicknesses at 10.67 s, 11.91 s, and 12.10 s, which were obtained by the horizontal displacement-time history curves. The horizontal displacement of each node at the left and right edges of the model can also be derived at these moments. Then, taking the horizontal displacement with different thicknesses of the loess cover at the time into the corresponding node in the static analysis model, the plastic strain nephograms and safety factors of the loess tunnel with different thicknesses of loess in the critical failure state are shown in Figure 8.
(a)
(b)
(c)
From Figure 8, when the tunnel with different overburden thicknesses is in the critical failure state, the safety factor increases with increasing burial depth. The safety factor of the shallow buried tunnel is 2.265, and the safety factor of the deep buried tunnel is 2.266. The safety factor of the superdeep buried tunnel is 2.669. The plastic strain is also symmetrically distributed on both sides of the vault and arch-toe, and the maximum value occurs at the arch-toe part.
3.3. Seismic Stability under Hard Rain Action
Due to the abundant rainfall in the rainy season, it is necessary to study the seismic stability of the loess mass under hard rain action. The rainfall intensity of hard rain is set to 75 mm/d. The rainwater infiltration depth is 2.0 m. The overburden thickness is taken as 30 m, 60 m, or 80 m. The horizontal displacement-time history curves of the model’s top-right corner vertices of 456, 261, and 239 can be obtained for different loess overlying thicknesses. The results are shown in Figure 9.
(a)
(b)
(c)
From Figure 9, the maximum horizontal displacement values of the tunnel model with different loess cover thicknesses at 10.67 s, 11.91 s, and 12.10 s can be obtained by the horizontal displacement-time history curves. The horizontal displacement of each node at the left and right edges of the model can be derived for each moment. Taking the horizontal displacement with different loess cover thicknesses at each moment into the corresponding node in the static analysis model, the plastic strain nephograms of the loess tunnel in the critical failure state are shown in Figure 10.
(a)
(b)
(c)
From Figure 10, when the tunnel with different overburden thicknesses under the action of hard rain and seepage is in the critical failure state, the safety factor decreases with decreasing burial depth. The safety factor of the shallow buried tunnel is 2.261, and the safety factor of the deep buried tunnel is 2.266. The safety factor of the superdeep buried tunnel is 2.664. The partial plastic strain is distributed on both sides of the vault, and the main damage area occurs on both sides of the arch-toe before the critical damage state occurs.
The contrast results of the safety factor under different rainfall conditions can be obtained by contrast analysis of the factors that affect the tunnel stability. The results are shown in Figure 11.
From Figure 11, the following observations can be made.
The safety factor of the loess mass increases with increasing tunnel depth under the critical failure state and the same rainfall intensity as the loess tunnel with different cover thicknesses. The safety factor of the tunnel increased from 2.267 to 2.278 and finally reached 2.673 under moderate rain seepage action. The safety factor of the tunnel increased from 2.265 to 2.266 and finally reached 2.669 under heavy rain seepage action. The safety factor of the tunnel increased from 2.258 to 2.261 and finally reached 2.664 under hard rain seepage action. The plastic strain also increases with increasing depth, and the distribution range of the plastic strain increases. The plastic strain is mainly distributed on both sides of the arch and toe, and the maximum value occurs at both sides of the feet.
For the loess mass of the tunnel with the same overburden thickness, considering the influence of rainfall intensity on the tunnel stability, the safety factor of the loess tunnel decreases with increasing infiltration depth. The safety factor of the tunnel decreased from 2.267 to 2.265, finally reaching 2.258 in the case of the shallow buried tunnel. The safety factor of the tunnel decreased from 2.278 to 2.266, finally reaching 2.261 in the case of the deep buried tunnel. The safety factor of the tunnel decreased from 2.673 to 2.669, finally reaching 2.664 in the case of the superdeep buried tunnel. The two points mentioned above show that the stability of the tunnel is directly related to the overburden thickness and the load action of the tunnel.
4. Conclusions
Under the common action of rainfall seepage and near-field pulse ground motion, the safety factor of a tunnel with different burial depths increases with increasing overburden thickness. Under different rainfall intensities, the safety factor of a loess tunnel structure decreases with increasing rainfall. With different rainfall intensities and loess cover thicknesses, the distribution range of the critical strain nephogram of the loess tunnel structure is located on both sides of the arch-toe and vault. The maximum strain value also occurs at the toe of the arch. Under the common action of rainfall seepage and near-field pulse ground motion for the same rainfall intensity and loess thickness, the safety factor of the surrounding tunnel mass structure is reduced to some extent. However, the distribution range and maximum plastic strain are increased to some extent. Under the common action of rainfall seepage and near-field pulse ground motion, at the same rainfall intensity, when the loess cover thickness is less than 60 m, the change in the safety factor is not obvious. After the loess cover thickness exceeds 60 m, the safety factor increases linearly with increasing loess thickness. However, the change in the safety coefficient with different rainfall intensities is very small under the same cover thickness.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This paper was supported by the National Natural Science Foundation of China (Grant numbers: 51478212 and 51968045) and the Education Ministry Doctoral Tutor Foundation of China (Grant number: 20136201110003).