# Seismic Stability of Loess Tunnel with Rainfall Seepage

Authors: Xuansheng Cheng; Wanlin Zhang; Jing Fan; Xiaoyan Zhang; Haibo Liu; Xinhai Zhou
Journal: Advances in Civil Engineering
 (2020)

---

## 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.


---

## Body

## 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]:(1)KBN=αNGR,CBN=ρcp,(2)KBT=αTGR,CBT=ρcS,where KBN is the normal spring stiffness, KBT is the tangential spring stiffness, CBN is the damping coefficient of normal dampers, CBT is the damping coefficient of tangential dampers, R is the distance from the wave beginning to the artificial boundary point, G is the medium shear modulus, cp is the P wave velocity of the medium, cs is the S wave velocity of the medium, αNis the normal viscoelastic boundary correction coefficient, and αT is the tangential viscoelastic boundary correction coefficient. In general, αN = 0.8–1.2, and αT = 0.35–0.65. In this paper, αN = 1.0, and αT = 0.5. cp and cs can be calculated according to the following formula:(3)cp=E1−μρ1+μ1−2μ,cs=Gρ,(4)G=E21+μ,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 as(5)fσ,τ=τ=c+σtanφ,where σ 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 tanφ 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,(6)C=αM+βK,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](7)ζi=α2ωi+βωi2i=1,2,…,n,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):(8)α2ω1+β⋅ω12=ζ;α2ω2+β⋅ω22=ζ.From equation (8),(9)α=2⋅ω1⋅ω2⋅ζω1+ω2;β=2⋅ζω1+ω2.(2) Dynamic Equation and Solution. The matrix differential equation of a structure under earthquake action can be expressed as [33,34](10)Mu¨t+Cu˙t+Kut=−Mu¨gt+pf,where u¨t, u˙t, and ut are the acceleration, velocity, and displacement vectors of a model node, respectively; u¨g is the earthquake acceleration; and pf is the vector of the surface load.The Newmark numerical integration is applied to equation (10), and we assume that(11)ut+Δt=ut+Δtu˙t+12−ςΔt2u¨t+ςΔt2u¨t+Δt,(12)u˙t+Δt=u˙t+1−τΔtu¨t+τΔtu¨t+Δt,where χ and δ are the constant values. At the time of t+Δt, the differential equation of motion is(13)Mu¨t+Δt+Cu˙t+Δt+Kut+Δt=−Mu¨gt+Δt+pf.Takingχ=1/2, δ=1/4, and Δt≤0.01Tmax (where Tmaxis 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(14)M+Δt2Cu¨t+Δt+Cu˙t+Δt2u¨t+Kut+Δt=−Mugt+Δt+pf.From equation (13),(15)u¨t+Δt=4Δt2ut+Δt−ut−4Δtu˙t−u¨t.Combining equation (15) with equation (14) yields(16)K+2ΔtC+4Δt2Mut+Δt=C2Δtut+u˙t+M4Δt2ut+4Δtu˙t+u˙t−Mu¨gt+Δt+pf.From equation (16), the value of ut+Δt 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:(17)c′=cη,(18)φ′=arctantanφη.Namely,(19)τ=c+σtanφ′=cη+σtanφ′,where c′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 as(20)F=J2+αI1−K=0,where I1 is the first invariant of the stress tensor, I1=σ1+σ2+σ3, and J2is the second invariant of the stress tensor, J2=1/6σ1−σ2+σ2−σ3+σ3−σ1.In the Mohr–Coulomb model, the parametersαand Kcan be expressed by using c (cohesion force) and φ (internal friction angle):(21)α=sinφ′33+sin2φ′,(22)K=3c′cosφ′33+sin2φ′.


#### 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 as(23)υ=−kJ,where υ 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 byit:(24)it=−Dθw∂θw∂z+kθw,where ∂θw/∂z is the water content gradient and Dθw 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 is(25)−∂∂xiρυi+ω=∂∂tρn,where ρ is the water density, υi is the Darcy velocity, n 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. θ=nSω, where Sω is the saturation and 0≤Sω≤1. 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(26)−∂∂xiρυi+ω=∂∂tρnSω.Assuming that the Darcy laws of saturated and unsaturated seepage are both available, Darcy’s law in unsaturated seepage can be expressed as(27)υi=−kijkrθ∂H∂xj.Combining equation (27) with the continuity equation (26) yields(28)−∂∂xiρkijkrθ∂H∂xj+ω=∂∂tρnSω,where kij is the saturated permeability tensor and krθ 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, 0≤krθ≤1. In the saturation zone, krθ=1.H=h+z, where z is the position head and h is the pressure head. Considering the functional relationship among θ, h, and kr, krθ=krθh, and(29)∂θ∂t=∂θ∂h∂h∂t,C=∂θ/∂h is defined as the water holding capacity. Considering that the compressibility of water is very low, ρ can be regarded as a constant. The saturation, Sω, can be replaced by the water content volume, θ. Differential equations of the saturated and unsaturated seepage can be obtained after calculating equation (29):(30)∂∂xikijkrh∂h∂xj+ki3krh+ω=Ch+βSs∂h∂t,where the water holding capacity, C, is equal to 0 in the saturation zone and Ss 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, Ss=0, 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 as(31)∂∂xku∂h∂x+∂∂yku∂h∂y+ω=Ch+βSs∂h∂t,where the unsaturated permeability coefficient is ku, ku=kskrh, and ks 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 Figure1.Figure 1
The energy of point A in they direction of an arbitrary datum.Considering the energy state of point A, point A has a gravitational potential energy,Eg:(32)Eg=mwgy,where mw is the water quality at point A and y is the height of point A relative to the datum plane.The energy,Ep, at point A induced by pressure can be expressed as(33)Ep=mw∫0uwυwmwduw=mw∫0uwduwρw=mwuwρw,where Ep is the pressure energy, uw is the pore water pressure at point A, υw is the water flow velocity at point A, and ρw is the water density at point A.Therefore, the kinetic energy of pointA can be expressed as(34)Eυ=mwυw22.The total potential energy of pointA can be expressed as(35)E=mwgy+mwuwρw+mwυw22.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:(36)hw=y+uwρwg+υw22g.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,hw=0. 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 i=1. The osmotic resistance can be expressed as(37)f=−iγw=−iρwg=−ρwg.The work of overcoming resistance can be expressed as(38)Wf=fy=−ρwgy.With no water infiltration, the value ofυw is 0. Based on the conservation of energy,(39)ymax=uwρwg⋅mwρw−mw.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 Figure1, under the condition of the plane state, water seepage follows Darcy’s law. The stable seepage of an anisotropic continuum equation is(40)∂∂xkx∂h∂x+∂∂yky∂h∂y=0,where h is the energy head, kx is the permeability coefficient in the x direction, and ky 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:(41)hx,y,t0=h0x,y.Head boundary condition:(42)hx,y,t|Γ1=h1x,y,t.Flow boundary condition:(43)−kijkrh∂h∂xj+ki3krhniΓ2=qn.Overflow boundary condition:(44)−kijkrh∂h∂xj+ki3krhniΓ3≥0,hΓ3=0,where qn and qθ are the normal vectors, with the direction of the vertical cross section in the positive direction. For rainfall special flow boundary conditions, qn = it; ni is the outer normal direction cosine of the boundary, t0 is the initial time, Γ1 is the energy head boundary, Γ2 is the flow boundary, and Γ3 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(45)H=∑i=1MNiHi.According to the variation principle and finite element theory, equation (40) can be expressed in matrix form as follows:(46)kh=F,where k is the seepage matrix, h is the node head matrix, and F 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 Table1. The material parameters of loess before and after rainfall over a certain range are compared with the loess material parameters before rainfall (Table 2).Table 1
Material parameters.
Material categoryElastic modulus (GPa)Poisson’s ratioGravity (kN/m3)Cohesive force (MPa)Internal friction angle (°)Loess0.10.3518.50.14533Initial lining300.225——Secondary lining300.225——Table 2
Material parameters.
CategoryElastic modulus (MPa)Poisson’s ratioGravity (kN/m3)Cohesive force (MPa)Internal friction angle (°)Before rainfall850.3518.50.125After rainfall300.3515.580.07519.8


### 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, whereHdis 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. KBN = 1.26×106N⋅S/m, CBN = 3.2×105N/m, KBT = 0.63×106N⋅S/m, and CBT = 6.58×105N/m. 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.Figure 2
Sketch map of the dynamic analysis model (unit: m).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 Figure3.Figure 3
Dynamic analysis model.When considering the train load for the model in Figure3, 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.Figure 4
Model of train load analysis.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.


## 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]:(1)KBN=αNGR,CBN=ρcp,(2)KBT=αTGR,CBT=ρcS,where KBN is the normal spring stiffness, KBT is the tangential spring stiffness, CBN is the damping coefficient of normal dampers, CBT is the damping coefficient of tangential dampers, R is the distance from the wave beginning to the artificial boundary point, G is the medium shear modulus, cp is the P wave velocity of the medium, cs is the S wave velocity of the medium, αNis the normal viscoelastic boundary correction coefficient, and αT is the tangential viscoelastic boundary correction coefficient. In general, αN = 0.8–1.2, and αT = 0.35–0.65. In this paper, αN = 1.0, and αT = 0.5. cp and cs can be calculated according to the following formula:(3)cp=E1−μρ1+μ1−2μ,cs=Gρ,(4)G=E21+μ,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 as(5)fσ,τ=τ=c+σtanφ,where σ 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 tanφ 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,(6)C=αM+βK,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](7)ζi=α2ωi+βωi2i=1,2,…,n,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):(8)α2ω1+β⋅ω12=ζ;α2ω2+β⋅ω22=ζ.From equation (8),(9)α=2⋅ω1⋅ω2⋅ζω1+ω2;β=2⋅ζω1+ω2.(2) Dynamic Equation and Solution. The matrix differential equation of a structure under earthquake action can be expressed as [33,34](10)Mu¨t+Cu˙t+Kut=−Mu¨gt+pf,where u¨t, u˙t, and ut are the acceleration, velocity, and displacement vectors of a model node, respectively; u¨g is the earthquake acceleration; and pf is the vector of the surface load.The Newmark numerical integration is applied to equation (10), and we assume that(11)ut+Δt=ut+Δtu˙t+12−ςΔt2u¨t+ςΔt2u¨t+Δt,(12)u˙t+Δt=u˙t+1−τΔtu¨t+τΔtu¨t+Δt,where χ and δ are the constant values. At the time of t+Δt, the differential equation of motion is(13)Mu¨t+Δt+Cu˙t+Δt+Kut+Δt=−Mu¨gt+Δt+pf.Takingχ=1/2, δ=1/4, and Δt≤0.01Tmax (where Tmaxis 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(14)M+Δt2Cu¨t+Δt+Cu˙t+Δt2u¨t+Kut+Δt=−Mugt+Δt+pf.From equation (13),(15)u¨t+Δt=4Δt2ut+Δt−ut−4Δtu˙t−u¨t.Combining equation (15) with equation (14) yields(16)K+2ΔtC+4Δt2Mut+Δt=C2Δtut+u˙t+M4Δt2ut+4Δtu˙t+u˙t−Mu¨gt+Δt+pf.From equation (16), the value of ut+Δt 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:(17)c′=cη,(18)φ′=arctantanφη.Namely,(19)τ=c+σtanφ′=cη+σtanφ′,where c′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 as(20)F=J2+αI1−K=0,where I1 is the first invariant of the stress tensor, I1=σ1+σ2+σ3, and J2is the second invariant of the stress tensor, J2=1/6σ1−σ2+σ2−σ3+σ3−σ1.In the Mohr–Coulomb model, the parametersαand Kcan be expressed by using c (cohesion force) and φ (internal friction angle):(21)α=sinφ′33+sin2φ′,(22)K=3c′cosφ′33+sin2φ′.


### 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.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 as(5)fσ,τ=τ=c+σtanφ,where σ 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 tanφ 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,(6)C=αM+βK,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](7)ζi=α2ωi+βωi2i=1,2,…,n,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):(8)α2ω1+β⋅ω12=ζ;α2ω2+β⋅ω22=ζ.From equation (8),(9)α=2⋅ω1⋅ω2⋅ζω1+ω2;β=2⋅ζω1+ω2.(2) Dynamic Equation and Solution. The matrix differential equation of a structure under earthquake action can be expressed as [33,34](10)Mu¨t+Cu˙t+Kut=−Mu¨gt+pf,where u¨t, u˙t, and ut are the acceleration, velocity, and displacement vectors of a model node, respectively; u¨g is the earthquake acceleration; and pf is the vector of the surface load.The Newmark numerical integration is applied to equation (10), and we assume that(11)ut+Δt=ut+Δtu˙t+12−ςΔt2u¨t+ςΔt2u¨t+Δt,(12)u˙t+Δt=u˙t+1−τΔtu¨t+τΔtu¨t+Δt,where χ and δ are the constant values. At the time of t+Δt, the differential equation of motion is(13)Mu¨t+Δt+Cu˙t+Δt+Kut+Δt=−Mu¨gt+Δt+pf.Takingχ=1/2, δ=1/4, and Δt≤0.01Tmax (where Tmaxis 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(14)M+Δt2Cu¨t+Δt+Cu˙t+Δt2u¨t+Kut+Δt=−Mugt+Δt+pf.From equation (13),(15)u¨t+Δt=4Δt2ut+Δt−ut−4Δtu˙t−u¨t.Combining equation (15) with equation (14) yields(16)K+2ΔtC+4Δt2Mut+Δt=C2Δtut+u˙t+M4Δt2ut+4Δtu˙t+u˙t−Mu¨gt+Δt+pf.From equation (16), the value of ut+Δt 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:(17)c′=cη,(18)φ′=arctantanφη.Namely,(19)τ=c+σtanφ′=cη+σtanφ′,where c′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 as(20)F=J2+αI1−K=0,where I1 is the first invariant of the stress tensor, I1=σ1+σ2+σ3, and J2is the second invariant of the stress tensor, J2=1/6σ1−σ2+σ2−σ3+σ3−σ1.In the Mohr–Coulomb model, the parametersαand Kcan be expressed by using c (cohesion force) and φ (internal friction angle):(21)α=sinφ′33+sin2φ′,(22)K=3c′cosφ′33+sin2φ′.


## 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 as(23)υ=−kJ,where υ 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 byit:(24)it=−Dθw∂θw∂z+kθw,where ∂θw/∂z is the water content gradient and Dθw 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 is(25)−∂∂xiρυi+ω=∂∂tρn,where ρ is the water density, υi is the Darcy velocity, n 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. θ=nSω, where Sω is the saturation and 0≤Sω≤1. 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(26)−∂∂xiρυi+ω=∂∂tρnSω.Assuming that the Darcy laws of saturated and unsaturated seepage are both available, Darcy’s law in unsaturated seepage can be expressed as(27)υi=−kijkrθ∂H∂xj.Combining equation (27) with the continuity equation (26) yields(28)−∂∂xiρkijkrθ∂H∂xj+ω=∂∂tρnSω,where kij is the saturated permeability tensor and krθ 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, 0≤krθ≤1. In the saturation zone, krθ=1.H=h+z, where z is the position head and h is the pressure head. Considering the functional relationship among θ, h, and kr, krθ=krθh, and(29)∂θ∂t=∂θ∂h∂h∂t,C=∂θ/∂h is defined as the water holding capacity. Considering that the compressibility of water is very low, ρ can be regarded as a constant. The saturation, Sω, can be replaced by the water content volume, θ. Differential equations of the saturated and unsaturated seepage can be obtained after calculating equation (29):(30)∂∂xikijkrh∂h∂xj+ki3krh+ω=Ch+βSs∂h∂t,where the water holding capacity, C, is equal to 0 in the saturation zone and Ss 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, Ss=0, 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 as(31)∂∂xku∂h∂x+∂∂yku∂h∂y+ω=Ch+βSs∂h∂t,where the unsaturated permeability coefficient is ku, ku=kskrh, and ks 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 Figure1.Figure 1
The energy of point A in they direction of an arbitrary datum.Considering the energy state of point A, point A has a gravitational potential energy,Eg:(32)Eg=mwgy,where mw is the water quality at point A and y is the height of point A relative to the datum plane.The energy,Ep, at point A induced by pressure can be expressed as(33)Ep=mw∫0uwυwmwduw=mw∫0uwduwρw=mwuwρw,where Ep is the pressure energy, uw is the pore water pressure at point A, υw is the water flow velocity at point A, and ρw is the water density at point A.Therefore, the kinetic energy of pointA can be expressed as(34)Eυ=mwυw22.The total potential energy of pointA can be expressed as(35)E=mwgy+mwuwρw+mwυw22.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:(36)hw=y+uwρwg+υw22g.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,hw=0. 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 i=1. The osmotic resistance can be expressed as(37)f=−iγw=−iρwg=−ρwg.The work of overcoming resistance can be expressed as(38)Wf=fy=−ρwgy.With no water infiltration, the value ofυw is 0. Based on the conservation of energy,(39)ymax=uwρwg⋅mwρw−mw.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 Figure1, under the condition of the plane state, water seepage follows Darcy’s law. The stable seepage of an anisotropic continuum equation is(40)∂∂xkx∂h∂x+∂∂yky∂h∂y=0,where h is the energy head, kx is the permeability coefficient in the x direction, and ky 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:(41)hx,y,t0=h0x,y.Head boundary condition:(42)hx,y,t|Γ1=h1x,y,t.Flow boundary condition:(43)−kijkrh∂h∂xj+ki3krhniΓ2=qn.Overflow boundary condition:(44)−kijkrh∂h∂xj+ki3krhniΓ3≥0,hΓ3=0,where qn and qθ are the normal vectors, with the direction of the vertical cross section in the positive direction. For rainfall special flow boundary conditions, qn = it; ni is the outer normal direction cosine of the boundary, t0 is the initial time, Γ1 is the energy head boundary, Γ2 is the flow boundary, and Γ3 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(45)H=∑i=1MNiHi.According to the variation principle and finite element theory, equation (40) can be expressed in matrix form as follows:(46)kh=F,where k is the seepage matrix, h is the node head matrix, and F 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.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 as(23)υ=−kJ,where υ 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 byit:(24)it=−Dθw∂θw∂z+kθw,where ∂θw/∂z is the water content gradient and Dθw 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 is(25)−∂∂xiρυi+ω=∂∂tρn,where ρ is the water density, υi is the Darcy velocity, n 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. θ=nSω, where Sω is the saturation and 0≤Sω≤1. 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(26)−∂∂xiρυi+ω=∂∂tρnSω.Assuming that the Darcy laws of saturated and unsaturated seepage are both available, Darcy’s law in unsaturated seepage can be expressed as(27)υi=−kijkrθ∂H∂xj.Combining equation (27) with the continuity equation (26) yields(28)−∂∂xiρkijkrθ∂H∂xj+ω=∂∂tρnSω,where kij is the saturated permeability tensor and krθ 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, 0≤krθ≤1. In the saturation zone, krθ=1.H=h+z, where z is the position head and h is the pressure head. Considering the functional relationship among θ, h, and kr, krθ=krθh, and(29)∂θ∂t=∂θ∂h∂h∂t,C=∂θ/∂h is defined as the water holding capacity. Considering that the compressibility of water is very low, ρ can be regarded as a constant. The saturation, Sω, can be replaced by the water content volume, θ. Differential equations of the saturated and unsaturated seepage can be obtained after calculating equation (29):(30)∂∂xikijkrh∂h∂xj+ki3krh+ω=Ch+βSs∂h∂t,where the water holding capacity, C, is equal to 0 in the saturation zone and Ss 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, Ss=0, 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 as(31)∂∂xku∂h∂x+∂∂yku∂h∂y+ω=Ch+βSs∂h∂t,where the unsaturated permeability coefficient is ku, ku=kskrh, and ks 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 Figure1.Figure 1
The energy of point A in they direction of an arbitrary datum.Considering the energy state of point A, point A has a gravitational potential energy,Eg:(32)Eg=mwgy,where mw is the water quality at point A and y is the height of point A relative to the datum plane.The energy,Ep, at point A induced by pressure can be expressed as(33)Ep=mw∫0uwυwmwduw=mw∫0uwduwρw=mwuwρw,where Ep is the pressure energy, uw is the pore water pressure at point A, υw is the water flow velocity at point A, and ρw is the water density at point A.Therefore, the kinetic energy of pointA can be expressed as(34)Eυ=mwυw22.The total potential energy of pointA can be expressed as(35)E=mwgy+mwuwρw+mwυw22.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:(36)hw=y+uwρwg+υw22g.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,hw=0. 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 i=1. The osmotic resistance can be expressed as(37)f=−iγw=−iρwg=−ρwg.The work of overcoming resistance can be expressed as(38)Wf=fy=−ρwgy.With no water infiltration, the value ofυw is 0. Based on the conservation of energy,(39)ymax=uwρwg⋅mwρw−mw.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 Figure1, under the condition of the plane state, water seepage follows Darcy’s law. The stable seepage of an anisotropic continuum equation is(40)∂∂xkx∂h∂x+∂∂yky∂h∂y=0,where h is the energy head, kx is the permeability coefficient in the x direction, and ky 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:(41)hx,y,t0=h0x,y.Head boundary condition:(42)hx,y,t|Γ1=h1x,y,t.Flow boundary condition:(43)−kijkrh∂h∂xj+ki3krhniΓ2=qn.Overflow boundary condition:(44)−kijkrh∂h∂xj+ki3krhniΓ3≥0,hΓ3=0,where qn and qθ are the normal vectors, with the direction of the vertical cross section in the positive direction. For rainfall special flow boundary conditions, qn = it; ni is the outer normal direction cosine of the boundary, t0 is the initial time, Γ1 is the energy head boundary, Γ2 is the flow boundary, and Γ3 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(45)H=∑i=1MNiHi.According to the variation principle and finite element theory, equation (40) can be expressed in matrix form as follows:(46)kh=F,where k is the seepage matrix, h is the node head matrix, and F 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 Table1. The material parameters of loess before and after rainfall over a certain range are compared with the loess material parameters before rainfall (Table 2).Table 1
Material parameters.
Material categoryElastic modulus (GPa)Poisson’s ratioGravity (kN/m3)Cohesive force (MPa)Internal friction angle (°)Loess0.10.3518.50.14533Initial lining300.225——Secondary lining300.225——Table 2
Material parameters.
CategoryElastic modulus (MPa)Poisson’s ratioGravity (kN/m3)Cohesive force (MPa)Internal friction angle (°)Before rainfall850.3518.50.125After rainfall300.3515.580.07519.8


## 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, whereHdis 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. KBN = 1.26×106N⋅S/m, CBN = 3.2×105N/m, KBT = 0.63×106N⋅S/m, and CBT = 6.58×105N/m. 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.Figure 2
Sketch map of the dynamic analysis model (unit: m).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 Figure3.Figure 3
Dynamic analysis model.When considering the train load for the model in Figure3, 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.Figure 4
Model of train load analysis.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 Figure5.Figure 5
Horizontal displacement-time history curves for (a) the 30  m layer at node 818, (b) the 60 m layer at node 596, and (c) the 80 m layer at node 4097.
(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 Figure6.Figure 6
Plastic strain nephograms of the rock surrounding the tunnel with different loess cover thicknesses: (a) 30 m,η = 2.251; (b) 60 m, η = 2.269; (c) 80 m; η = 2.278.
(a)(b)(c)From Figure6, 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 Figure7.Figure 7
Horizontal displacement-time history curves for (a) the 30 m layer at node 398, (b) the 60 m layer at node 261, and (c) the 80 m layer at node 239.
(a)(b)(c)Figure7 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.Figure 8
Plastic strain nephograms of the loess mass surrounding the tunnel for different overburden thicknesses: (a) 30 m,η = 2.265; (b) 60 m, η = 2.266; (c) 80 m, η = 2.669.
(a)(b)(c)From Figure8, 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 Figure9.Figure 9
Horizontal displacement-time history curves for (a) the 30 m layer at node 456, (b) the 60 m layer at node 261, and (c) the 80 m layer at node 239.
(a)(b)(c)From Figure9, 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.Figure 10
Plastic strain nephograms of the rock surrounding the tunnel with different overburden thicknesses: (a) 30 m,η = 2.258; (b) 60 m, η = 2.261; (c) 80 m, η = 2.664
(a)(b)(c)From Figure10, 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 Figure11.Figure 11
Safety factor of the loess tunnel with different rainfall intensities.From Figure11, 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.


## 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 Figure5.Figure 5
Horizontal displacement-time history curves for (a) the 30  m layer at node 818, (b) the 60 m layer at node 596, and (c) the 80 m layer at node 4097.
(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 Figure6.Figure 6
Plastic strain nephograms of the rock surrounding the tunnel with different loess cover thicknesses: (a) 30 m,η = 2.251; (b) 60 m, η = 2.269; (c) 80 m; η = 2.278.
(a)(b)(c)From Figure6, 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 Figure7.Figure 7
Horizontal displacement-time history curves for (a) the 30 m layer at node 398, (b) the 60 m layer at node 261, and (c) the 80 m layer at node 239.
(a)(b)(c)Figure7 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.Figure 8
Plastic strain nephograms of the loess mass surrounding the tunnel for different overburden thicknesses: (a) 30 m,η = 2.265; (b) 60 m, η = 2.266; (c) 80 m, η = 2.669.
(a)(b)(c)From Figure8, 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 Figure9.Figure 9
Horizontal displacement-time history curves for (a) the 30 m layer at node 456, (b) the 60 m layer at node 261, and (c) the 80 m layer at node 239.
(a)(b)(c)From Figure9, 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.Figure 10
Plastic strain nephograms of the rock surrounding the tunnel with different overburden thicknesses: (a) 30 m,η = 2.258; (b) 60 m, η = 2.261; (c) 80 m, η = 2.664
(a)(b)(c)From Figure10, 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 Figure11.Figure 11
Safety factor of the loess tunnel with different rainfall intensities.From Figure11, 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.



---

*Source: 8147950-2020-09-10.xml*