Stability of Saline Soil Subgrade Slopes under Surface Runoff and Seismic Action: A Case Study of the G0711 Urumqi–Yuli Expressway, China

1  Introduction

Subgrade slope stability refers to the ability of the subgrade slope to remain stable under external loads. In road engineering, maintaining the stability of subgrade slopes is crucial. Once slope stability issues arise, they can lead to subgrade collapse, road damage, and even severe traffic accidents. Saline soil is a general term for salt soil, alkali soil, and various degrees of salinized and alkalized soils, typically referring to soils with a soluble salt content greater than 0.3%, which exhibit undesirable properties such as subsidence, salt expansion, and corrosion [1]. Saline soil is mainly distributed in the arid regions of China, such as Xinjiang, Gansu, Qinghai, and the Qinghai-Tibet Plateau, where the climate is dry, precipitation is scarce, and the salt concentration in groundwater is high, leading to soil salt accumulation. Saline soil is a special type of soil containing soluble salts. As the soil moisture content decreases, the water in the soil crystallizes, and these salt crystals act as a skeleton between soil pores. When the moisture content increases, the soluble salts dissolve completely in the soil water, causing structural damage and reducing soil strength [2].

Currently, many scholars still primarily employ traditional analytical methods, such as limit equilibrium methods and strength reduction techniques, to evaluate the stability of subgrades and slopes under rainfall or seismic conditions. These methods have been widely used in engineering practice because of their clear physical meaning, high computational efficiency, and relatively simple parameter determination. Moreover, traditional approaches have been extensively adopted in design codes and engineering guidelines, providing strong practical applicability and operability. Alonso et al. [3] investigated the evolution of the slope safety factor under various rainfall histories, analyzing the effects of rainfall duration, rainfall amount, soil-water characteristic curve, and soil permeability on slope stability. A sensitivity analysis was also conducted, revealing that a higher soil permeability coefficient leads to more rapid changes in the safety factor. Rahardjo et al. [4] further simulated the variations in pore water pressure induced by various rainfall patterns and calculated the corresponding changes in the slope safety factor. The study found that initial rainfall has a significant impact on slope stability. Dou et al. [5] established a rainwater infiltration-redistribution model that accounts for the variability of saturated permeability coefficient, and proposed a probabilistic framework for rainfall-induced landslides, offering a new perspective for understanding the complexity of soil moisture movement. Gavin and Xue [6] introduced a simplified analytical method for rainfall infiltration in unsaturated soils. By developing a modified Green-Ampt model that incorporates the influence of matric suction, they achieved rapid prediction of wetting front progression. Sajjan et al. [7] proposed an improved coupled hydro-mechanical analysis method. Through the development of a distributed hydrological model that considers slope geometry and fill permeability characteristics, they systematically analyzed the runoff generation mechanism and stability evolution of embankment slopes under intense rainfall conditions. Song and Li [8] used ABAQUS’s pore fluid seepage function to calculate the impact of rainfall infiltration on subgrade slope stability, studying the effects of rainfall infiltration on the seepage field, stress field, and displacement of slopes. Guo and Zhao [9] combined experiments with numerical simulations to reveal the mechanical behavior and failure modes of an expansive soil subgrade under rainfall conditions, analyzing the effects of rainfall intensity, initial matric suction, and slope ratio. Fan et al. [10] found that rainfall infiltration leads to soil softening, reduced matric suction, and strength reduction, seriously threatening the safety of slopes and tunnels. Liu et al. [11] identified rainfall intensity, rainfall duration, and soil long-term strength as the main factors affecting slope stability under rainfall. Gao et al. [12] proposed that composite embankments can effectively reduce frost heave and uneven deformation compared to ordinary embankments under rainfall infiltration conditions, avoiding longitudinal cracks. These studies did not consider the minimal impact of low rainfall intensity on the moisture content of saline soil subgrade slopes in arid regions. Concentrated rainfall in arid regions can affect the moisture content of the shallow soil layers of saline soil subgrade slopes, leading to surface runoff and posing significant risks to subgrade slopes. Kanji and Das [13] demonstrated that the coupling between rainfall intensity and runoff processes governs the advancement depth of the wetting front and the spatial variability of soil moisture, thereby exerting a significant influence on slope stability and near-surface deformation. They further pointed out that, in arid and semi-arid regions, rainfall–runoff pathways play a critical role in regulating soil hydraulic responses and stress evolution induced by changes in soil moisture [14]. However, these methods usually involve simplified treatments of moisture distribution and strength degradation processes, making it difficult to accurately capture the coupled effects of seepage, strength evolution, and dynamic loading. This limitation highlights the need to develop analytical approaches that explicitly consider hydro-mechanical coupling mechanisms in order to more reliably assess subgrade slope stability.

As one of the most destructive natural disasters, earthquakes can severely affect slope stability in a short time, triggering secondary disasters such as landslides and mudflows. Rajabian and Vahedifard [15] proposed using ground anchors to effectively improve the seismic stability of concave slopes. Liao et al. [16] proposed a simulation-based evaluation method to analyze the impact of mainshocks and aftershocks on rock slope stability, finding that frequency domain parameters are particularly effective in assessing rock slope stability changes under high-magnitude earthquakes. Wu et al. [17] analyzed the effects of different factors on the safety factor of anchored slopes under tension-shear coupling and tension-only conditions using limit equilibrium theory and a modified pseudo-dynamic method. Wang et al. [18] found that under seismic conditions, the potential sliding surface with asynchronous strength parameter reduction is shallower than under static conditions but deeper than without asynchronous strength parameter reduction. Zhou et al. [19] used a three-dimensional dynamic finite difference method to establish a simplified rock slope model, studying the dynamic response mechanism of bedding slopes with alternately distributed soft and hard rock layers under different seismic wave excitation directions. Liu et al. [20] studied the stability of partially saturated soil slopes under the combined effects of seismic forces and pore water pressure using the kinematic method of limit analysis. Konovalov and Gensiorovsky [21] proposed a probabilistic random method for normalizing internal slope deformation to assess slope stability. Leshchinsky and San [22] investigated the stability of simple slopes under seismic loading using the pseudo-static method. By assuming a predetermined logarithmic spiral failure surface, they determined the normal stress distribution through numerical methods, and further derived the minimum safety factor satisfying all limit equilibrium equations. Hall [23] analyzed train vibration data under seismic action using a mathematical numerical model, finding relatively large shear stresses in soft soil layers. Massey et al. [24] conducted field investigations on rock slopes under strong seismic action to study their dynamic response characteristics, and performed numerical simulation analysis accordingly. Giri and Sengupta [25] carried out dynamic response tests on slopes and, in combination with the FLAC 3D finite element software, conducted computational studies on the dynamic response characteristics of slopes under earthquake conditions. Feng et al. [26] proposed a dynamic–numerical coupled approach capable of simultaneously simulating the seismic response of landslides. By coupling the PFC and FLAC programs, they achieved integrated numerical simulation. Cheng et al. [27] explored the impact of seismic action on the stability of saline soil subgrades under the influence of temperature on the strength of saline soil subgrade fillers. These studies considered the impact of earthquakes on slope stability but did not consider the significant impact of moisture content changes on the strength parameters of saline soil subgrade slopes due to their high salt content under surface runoff. Therefore, it is necessary to consider the impact of strength parameter changes on saline soil subgrade slopes when subjected to surface runoff and earthquakes.

Although previous studies have extensively investigated subgrade and slope stability under rainfall or seismic loading using traditional analytical methods and numerical techniques, several limitations remain. Most existing studies simplify the seepage process by assuming homogeneous permeability or averaged hydraulic parameters, making it difficult to capture the spatial heterogeneity of moisture migration in layered saline soils. Moreover, the coupling between rainfall-induced moisture redistribution and subsequent strength degradation is often treated implicitly or indirectly, rather than being explicitly incorporated into stability analysis frameworks. In addition, empirical relationships between moisture content and shear strength parameters have mainly been used as standalone laboratory results, while their integration with spatially distributed seepage fields in numerical stability analyses remains limited. As a result, the link between moisture migration, strength evolution, and slope stability is not fully established, particularly for saline soil subgrade slopes in arid regions subjected to combined rainfall and seismic effects.

To address these gaps, the present study proposes an integrated seepage–stability analysis framework that explicitly links unsaturated seepage simulation, moisture-dependent strength degradation, and numerical slope stability analysis. The key novelties of this work are as follows:

(1) A layered equivalence method is developed to systematically transform spatial moisture content fields obtained from COMSOL simulations into layer-specific strength parameters for FLAC3D modeling;

(2) Rainfall-induced strength degradation is explicitly considered by incorporating experimentally established moisture content–strength relationships into the numerical stability analysis; and

(3) The coupled effects of rainfall infiltration and seismic loading on saline soil subgrade slope stability are quantitatively evaluated using a representative expressway case in an arid region.

These contributions provide a transparent and physically consistent approach for assessing subgrade slope stability under coupled hydro-mechanical conditions and offer practical insights for engineering design in saline soil regions.

2  Engineering Background

The G0711 Urumqi–Yuli Expressway was selected as the study area due to its typical geological, climatic, and engineering characteristics, making it representative of saline soil subgrade slopes in arid regions of northwest China, as shown in Fig. 1. The expressway traverses areas dominated by sulfate and chloride-type saline soils and is exposed to a typical arid continental climate characterized by scarce annual precipitation, strong evaporation, and highly concentrated short-term rainfall events.

Figure 1: G0711 Urumqi–Yuli expressway.

The Korla section of the G0711 Expressway has experienced engineering problems such as slope surface erosion, shallow instability, and subgrade deformation induced by heavy rainfall. These issues are further exacerbated by the frequent seismic activity in the Xinjiang region, where moderate earthquakes occur with relatively high frequency. The combined effects of rainfall-induced moisture variation and seismic loading pose a significant threat to the long-term stability and service performance of saline soil subgrade slopes along this expressway. Therefore, the G0711 Urumqi–Yuli Expressway provides a practical and representative engineering case to investigate the coupled influence of surface runoff and seismic action on saline soil subgrade slope stability. The findings of this study are expected to provide technical reference and methodological support for the design, assessment, and protection of similar saline soil subgrade slopes in arid and seismically active regions.

The annual precipitation in Korla is less than 60 mm, while the annual evaporation is about 2800 mm. The monthly average precipitation and temperature from 1971 to 2024 were statistically analyzed, with the monthly average precipitation shown in Fig. 2 and the monthly average temperature in Fig. 3. Fig. 2 shows that precipitation is unevenly distributed throughout the year, with the highest average precipitation in July at 12.4 mm and the lowest in November at 0.6 mm. Fig. 3 shows that the highest monthly average temperature is in July at 26.4°C.

Figure 2: Monthly average rainfall in Korla.

Figure 3: Monthly average temperature in Korla.

Xinjiang is located in the northwestern frontier of China, at the junction of the Eurasian and Indian plates, making it one of the most earthquake-prone regions in China. Some earthquakes in the Xinjiang Uygur Autonomous Region are listed in Table 1, showing that although the magnitude of earthquakes in Xinjiang is not large, the frequency is high, making it necessary to design local roads for earthquake resistance. According to the China Seismic Ground Motion Parameter Zoning Map [28] (GB 18306-2015), the study area is classified as a Class II site, for which the basic seismic peak ground acceleration ranges from 0.1 to 0.2 g. Based on this zoning requirement, a peak ground acceleration of 0.2 g was adopted in this study to evaluate the seismic stability of the saline soil subgrade slope.

3  Shear Strength Test of Saline Soil

3.1 Physical Parameter Test of Soil Samples

The soil samples were taken from the saline soil along the G0711 Urumqi to Yuli Expressway, mainly sulfate or chloride saline soil. According to the “Standard for Soil Test Methods” (GB/T50123-1999), the selected soil samples were tested for relative density, liquid limit and plastic limit moisture content, maximum dry density and optimal moisture content, particle analysis, and soluble salt chemical tests. The physical characteristic parameters and salt content of the saline soil samples are shown in Tables 2 and 3 [29].

3.2 Shear Test of Soil Samples under Different Moisture Contents

To study the effect of moisture content on the strength of the subgrade soil in the test section, shear tests were conducted on the selected soil samples at 96% compaction under moisture content 19.9%, optimal moisture content and moisture content 13.9%. To ensure that the moisture content of the samples did not change, the samples prepared with ring knives were wrapped in plastic wrap. All samples were subjected to vertical stresses of 100, 200, 300, and 400 kPa. The stress-strain relationship curves of the soil samples under different moisture contents and vertical stresses are shown in Fig. 4 (where σ is the vertical stress) [29]. Under low vertical stress of 100 kPa, the soil samples exhibited significant strain softening behavior at different moisture contents.

Figure 4: Different moisture content stress-strain curve.

3.3 Effect of Moisture Content on Shear Strength Parameters

Through the shear tests, the changes in cohesion and internal friction angle of the soil samples under different moisture contents are shown in Fig. 5.

Figure 5: The variation curves of cohesion and internal friction angle with water content.

When the moisture content of the soil samples decreased by 3% from the optimal moisture content, the cohesion was 66.0 kPa, an increase of 89.7% compared to the optimal moisture content; when the moisture content decreased by 1.5%, the cohesion was 54.3 kPa, an increase of 29.5%; when the moisture content increased by 1.5%, the cohesion was 21.1 kPa, a decrease of 39.4%; when the moisture content increased by 3%, the cohesion decreased to 16.0 kPa, a decrease of 54.0%. It can be seen that as the moisture content decreased from 3% below the optimal moisture content to 3% above, the cohesion of the soil samples continuously decreased, with a total decrease of 312.5%. As the moisture content of the samples increased, the internal friction angle gradually decreased, ranging from 34.67° to 34.06°, with little overall change.

The above tests indicate that variations in moisture content mainly affect the cohesion of saline soil, whereas the internal friction angle remains relatively insensitive. This behavior is primarily attributed to changes in matric suction, which is defined as the difference between pore air pressure and pore water pressure in unsaturated soils [30]. Matric suction provides an additional apparent cohesion by increasing interparticle tensile forces and enhancing soil shear resistance. As the moisture content increases, matric suction gradually decreases due to the reduction in negative pore water pressure, leading to a significant decrease in cohesion. This mechanism explains the pronounced sensitivity of saline soil cohesion to moisture variation.

The moisture content–strength relationships adopted in this study are expressed by Eqs. (1) and (2), which were originally proposed in Reference [29] based on laboratory shear tests on saline soils. These equations are directly adopted here to quantify the dependence of soil cohesion and internal friction angle on moisture content.

c={−10.4Δw+34.25w∈[wopt−3% ,wopt]−6.93Δw+34.25w∈[wopt,wopt+3% ](1)

φ={0.06Δw+34.25w∈[wopt−3% ,wopt]−0.09Δw+34.25w∈[wopt,wopt+3% ](2)

where, c is the cohesion of the soil sample (unite: kPa); φ is the internal friction angle of the soil sample (unite: °); Δw is the increase or decrease in moisture content relative to the optimum moisture content (unite: %); w is the moisture content (unite: %); wopt is the optimum moisture content (unite: %).

Currently, many scholars use traditional methods to treat the subgrade as a whole and assign a single shear strength parameter value when studying subgrade slope stability. However, traditional slope stability analysis methods, such as the limit equilibrium method and conventional strength reduction techniques, do not fully consider the spatial variability of soil permeability coefficients. In these approaches, the infiltration process is commonly simplified by assuming homogeneous or averaged permeability parameters within the slope, or by indirectly accounting for rainfall effects through pore water pressure increments. As a result, the heterogeneity of permeability among different soil layers and its influence on moisture migration cannot be explicitly represented. This simplification may lead to an inaccurate description of the seepage process and, consequently, an incomplete assessment of rainfall-induced strength degradation and slope stability. Therefore, it is necessary to adopt analysis frameworks that explicitly incorporate spatial variations in permeability and their coupling with mechanical responses. Therefore, based on the segmented variation law between moisture content and shear strength parameters, a subgrade moisture content-strength parameter equivalent layered calculation model was introduced to better reflect the impact of moisture content changes on saline soil subgrade slopes.

4  Unsaturated Soil Seepage Strength Theory

Unsaturated soil seepage theory studies the seepage problem in soils where both gas and water phases coexist. In unsaturated soils, the pores contain both water and gas. The seepage law of unsaturated soils is different from that of saturated soils.

4.1 Unsaturated Soil Seepage Theory

In the analysis of the seepage field of subgrade soils, considering that rainwater infiltration is instantaneous, the compressibility of water and seepage consolidation are usually not considered. The partial differential equation for unsaturated soil seepage is as follows [29]:

∂∂x(kx∂h∂x)+∂∂z(ky∂h∂z)+W=mwρwg∂h∂t(3)

where, kx, ky are the permeability coefficients in x and y directions respectively; h is the total water head of rock fracture; W is source term; mw is specific water capacity; ρw is density of water.

It should be noted that Eq. (3), which governs unsaturated soil seepage in this study, is subject to several inherent limitations. First, the equation assumes isotropic and homogeneous hydraulic properties within each soil layer, and thus does not explicitly account for small-scale spatial variability in permeability. Second, hysteresis effects in the soil–water characteristic curve during wetting and drying cycles are neglected, which may influence moisture redistribution under repeated rainfall events. In addition, the formulation does not explicitly consider preferential flow paths or soil structural anisotropy that may develop in highly fissured or heterogeneous soils. Despite these limitations, Eq. (3) has been widely adopted in unsaturated seepage analyses due to its clear physical basis and computational efficiency. In the present study, its application is justified because the analysis focuses on rainfall-induced infiltration under short-term loading conditions, and the soil properties are represented using layer-averaged parameters derived from laboratory tests. Moreover, the primary objective of this study is to capture the macroscopic redistribution of soil moisture and its influence on strength degradation, rather than detailed pore-scale flow behavior.

4.2 Strength Reduction Method

The principle of the strength reduction method [31] is to keep the gravitational acceleration of the slope soil constant. Before calculation, a small value is selected as the initial reduction factor, and this value is used to continuously reduce the cohesion and internal friction angle of the soil. If convergence is not achieved or the plastic zone is not penetrated within a given number of iterations, the reduction factor is slowly increased, and the new cohesion and internal friction angle are divided by the new reduction factor for calculation until the plastic zone is penetrated or non-convergence occurs, indicating that the soil has reached a critical sliding state. The calculation formula can be expressed as [32]:

c′=cFr(4)

φ′=arctan⁡(1Frtan⁡φ)(5)

where, c and c′ are the cohesion of rock and soil before and after reduction respectively; φ and φ′ are the internal friction angles of rock and soil before and after reduction respectively; Fr is strength reduction factor.

According to the basic theory of elastoplastic mechanics, the initial deformation stage of rock and soil belongs to elastic deformation. Therefore, in the calculation process, a sufficiently small reduction factor Fr (this paper takes Fr0 = 0.8) is first selected for calculation, and then the value of Fr is gradually increased. By continuously reducing the shear strength parameters of the rock and soil, the calculation is performed until the plastic deformation zone of the slope is penetrated or the calculation result does not converge. It should be emphasized that, in the strength reduction analysis, the initial strength reduction factor Fr0 serves only as the starting value for the iterative procedure and is introduced to ensure numerical stability and convergence efficiency at the early stage of computation. Its selection does not affect the determination of the final critical strength reduction factor, as the slope safety factor is solely defined by the value of Fr corresponding to the onset of instability in the numerical model. Given that the slope remains clearly stable and no plastic zone develops when Fr = 1.0, initiating the reduction process from a value slightly lower than 1.0 facilitates the gradual establishment of a reasonable stress field and a smooth transition toward the critical failure state.

Currently, the judgment methods for the critical state of slopes mainly include the following three aspects: (1) Check whether the finite element program calculation converges; (2) Observe whether the displacement curve of the characteristic point shows a sudden change; (3) Confirm whether the plastic zone of the slope is penetrated.

5  Impact of Surface Runoff on the Stability of Saline Soil Subgrade Slopes

5.1 Introduction to Numerical Simulation Software

This section uses the numerical simulation software COMSOL Multi-physics, employing solid mechanics and unsaturated soil seepage knowledge to analyze the stability of saline soil subgrade slopes. The changes in the moisture field of the saline soil subgrade slope under surface runoff were analyzed, and the strength reduction method was used to calculate the stability of the saline soil subgrade slope under different moisture content layers.

5.2 Model Establishment

Taking the subgrade slope of the G0711 Urumqi to Yuli Expressway as the research object, a subgrade model was established. The road width is 27 m, and the slope is 1:1.5 when the subgrade height is less than or equal to 8 m; the slope is 1:1.75 when the subgrade height is between 8 and 20 m. The subgrade fill is saline soil with a compaction degree of 96%, and the foundation soil is silty clay. The complete mesh of the slope model consists of 698 domain units and 112 boundary units, as shown in Fig. 6. The initial moisture content of the subgrade slope is shown in Fig. 7.

Figure 6: Subgrade slope grid division model.

Figure 7: Initial water content of subgrade slope.

At the start of the simulation, the entire subgrade slope was assumed to be in an initial equilibrium state. The distribution of initial pore water pressure was established based on the region’s average annual climatic conditions and the typical state of the subgrade after construction. Specifically, the initial moisture content was distributed non-uniformly across the slope model, with the shallow soil moisture content close to the optimum water content (~16.9%), varying slightly with depth toward the foundation. This simulates the initial moisture condition of the subgrade under long-term natural conditions (the revised manuscript’s Fig. 7 shows a cloud diagram of the initial moisture content distribution).

Rainfall Infiltration Boundary: During the rainfall period (4 h), the slope top and both slope surfaces were defined as flux boundaries. Based on the set rainfall intensities (1, 2, 3, 4, 5 mm/h), the rainfall rate was converted into an equivalent water flux applied to these boundaries. We used the lesser value between the actual rainfall intensity and the soil’s saturated permeability coefficient as the infiltration flux. This aligns with the basic concept of the Green-Ampt model: when rainfall intensity is less than the soil’s infiltration capacity, all rainfall infiltrates; otherwise, the excess forms runoff. This directly models the conditions for surface runoff generation. Before and after the rainfall, these boundaries were set as zero-flux boundaries (i.e., no moisture exchange), simulating moisture migration driven only by gravitational potential and matric suction gradients under rain-free conditions.

Model Bottom Boundary: Assumed to be a relatively impermeable layer or where the groundwater table is extremely deep, thus it was set as a zero-flux boundary with no vertical moisture exchange. This is consistent with the arid geological background and deep groundwater table of the study area.

Model Lateral Boundaries: Considering the symmetry of the subgrade slope model and the dominance of vertical infiltration and slope surface runoff on moisture movement, the lateral boundaries were set as zero-flux boundaries, neglecting moisture exchange with the far-field sides.

5.3 Impact of Rainfall on the Surface Moisture Content of Saline Soil

Clearly states that the model is based on a typical cross-section of the subgrade from the G0711 Urumqi–Yuli Expressway, Road surface width: 27 m. Two standard ratios based on subgrade height—1:1.5 when height ≤ 8 m; 1:1.75 when height is between 8 and 20 m. Includes not only the embankment fill but also the underlying foundation soil (silty clay) to ensure reasonable boundary conditions and account for the interaction between the foundation and the subgrade.

Subgrade Fill: Specified as saline soil with a compaction degree of 96%. Its key hydraulic parameters (e.g., soil-water characteristic curve, permeability function for the Richards equation) and mechanical parameters (initial cohesion, internal friction angle, elastic modulus, Poisson’s ratio, etc.) are all derived from the laboratory test results described in Section 2 and reference [29].

Foundation Soil: Specified as silty clay, with its relevant geotechnical parameters selected as typical values based on regional geological survey reports and engineering experience.

Core Modeling Strategy—“Moisture Content-Strength Parameter Equivalent Layered Model”: This section explicitly states that the core innovation of the model lies in the fact that the strength parameters are not constants. In the simulation, we first calculated the spatial moisture content distribution after rainfall using COMSOL. Then, based on the c-ω and φ-ω relationships established in Fig. 5, Eqs. (1) and (2), the computational domain was divided into horizontal layers according to moisture content contours. Finally, the average moisture content of each layer was converted into strength parameters and assigned to the corresponding zones in the FLAC3D model. This fully describes the key steps in model establishment from the “moisture field” to the “strength field.” Mesh Information: Specifies the type and scale of mesh used, for example, stating that the final model consists of 698 domain elements and 112 boundary elements.

The change in moisture content of the saline soil subgrade is mainly affected by groundwater and rainfall. Rainfall increases the moisture content of the subgrade soil, thereby affecting the shear strength of the soil. Under different rainfall conditions, due to differences in soil permeability, the moisture content changes in different parts of the subgrade may be uneven, affecting the overall stability of the subgrade. To study the impact of different rainfall intensities on the moisture content of the saline soil subgrade slope, five rainfall intensities of 1, 2, 3, 4, and 5 mm/h were set, each lasting 4 h, to analyze their impact on the slope.

Fig. 8 illustrates the spatial distribution of moisture content in the saline soil subgrade slope under different rainfall intensities. Specifically, Fig. 8a–e correspond to rainfall intensities of 1, 2, 3, 4, and 5 mm/h, respectively.

Figure 8: Water content distribution cloud map of subgrade slope under different rainfall intensity.

As shown in Fig. 8a, when the rainfall intensity is 1 mm/h, the increase in moisture content is mainly confined to the shallow surface layer of the slope, with negligible variation in the deeper soil. In Fig. 8b,c, corresponding to rainfall intensities of 2 and 3 mm/h, the wetted zone gradually extends downslope and slightly deepens, indicating enhanced infiltration along the slope surface.

When the rainfall intensity increases to 4 and 5 mm/h, as shown in Fig. 8d,e, the moisture content in the shallow soil layer increases significantly, and localized surface runoff becomes evident. Although the overall moisture variation within the subgrade remains limited due to the arid climatic conditions, higher rainfall intensities clearly intensify moisture accumulation near the slope surface, which may adversely affect slope stability. Fig. 9 shows that as the rainfall intensity increases, the moisture content of the shallow soil layers of the subgrade slope increases from 16.7% to 19.7%.

Figure 9: The curve of shallow water content of subgrade slope changing with rainfall intensity.

5.4 Impact of Surface Runoff on Slope Stability

This paper considers the impact of surface runoff on the stability of the saline soil subgrade. The strength parameters of the subgrade soil at specific locations are determined based on moisture content, and the subgrade is divided into horizontal layers to calculate the safety factor of the subgrade slope and analyze its stability. Using the functional relationship between soil strength parameters and moisture content determined in the laboratory, a subgrade moisture content-strength parameter equivalent layered calculation model was adopted. Fig. 10 shows the schematic diagram of the subgrade moisture content layering.

Figure 10: Subgrade layering schematic model.

Equivalent plastic strain is a physical quantity that reflects the location of the yield strain of the slope material, representing the potential sliding surface at the critical failure of the slope [33]. Fig. 11 presents the equivalent plastic strain distributions of the saline soil subgrade slope after strength reduction under different rainfall intensities. Fig. 11a shows the equivalent plastic strain distribution corresponding to a rainfall intensity of 1 mm/h, while Fig. 11b represents the case with a rainfall intensity of 5 mm/h.

Figure 11: Equivalent plastic strain nephogram of subgrade slope under different rainfall intensities.

As illustrated in Fig. 11a, under low rainfall intensity, the plastic strain is mainly concentrated near the slope toe, and the plastic zone is relatively limited, indicating that the slope remains in a stable state. In contrast, Fig. 11b demonstrates that when the rainfall intensity increases to 5 mm/h, the equivalent plastic strain zone expands upward along the slope surface and develops into a continuous arc-shaped band. This arc-shaped plastic strain concentration zone can be identified as the potential sliding surface at the critical failure state of the slope. The comparison between Fig. 11a,b clearly indicates that increased rainfall intensity enhances moisture-induced strength degradation, thereby promoting the development and penetration of the plastic zone within the slope.

Fig. 12 shows the slope safety factor under different rainfall intensities. It can be seen that when the rainfall intensity is 1 mm/h, the slope safety factor is 2.41. As the rainfall intensity increases, the shallow water content of the slope not only increases but also leads to surface runoff, and the slope safety factor decreases. When the rainfall intensity reaches 5 mm/h, the slope safety factor drops to 1.98, a decrease of 17.8%.

Figure 12: Curves of slope safety factor and shallow water content of subgrade changing with rainfall intensity.

6  Influence of Surface Runoff and Seismic Action on the Stability of Saline Soil Subgrade Slopes

6.1 Introduction to Seismic Action Numerical Simulation Software

FLAC3D is an extension of the two-dimensional finite difference program FLAC2D, capable of simulating the mechanical properties and plastic flow analysis of three-dimensional structures in soil, rock, and other materials. It adjusts polyhedral units in a three-dimensional grid to fit actual structures. FLAC3D employs an explicit Lagrangian algorithm and mixed-discrete partitioning technology, allowing for accurate simulation of material plastic failure and flow. Since it does not require the formation of a stiffness matrix, it can solve large-scale three-dimensional problems with relatively small memory space, making it a widely referenced numerical simulation software in the geotechnical field.

6.2 Initial Geostress Analysis

To ensure physical consistency between the seepage analysis and the subsequent stability analysis, a layered equivalence approach was adopted to systematically convert the moisture content fields obtained from COMSOL into strength parameter inputs for the FLAC3D model. First, COMSOL Multiphysics was employed to simulate the unsaturated seepage process of the saline soil subgrade slope under different rainfall intensities, yielding the spatial distribution of moisture content within the slope after rainfall infiltration. Subsequently, based on the moisture content contour characteristics, the subgrade was divided into several horizontal layers with approximately uniform moisture content, and representative average moisture content values for each layer were extracted from the COMSOL results. This layered discretization effectively captures the spatial heterogeneity of moisture distribution within the subgrade induced by rainfall infiltration.

On this basis, the experimentally established relationships between moisture content and shear strength parameters were applied by substituting the average moisture content of each layer into the fitted functions to calculate the corresponding cohesion and internal friction angle. The laboratory results indicate that variations in moisture content predominantly affect soil cohesion, whereas the internal friction angle is only marginally influenced. Therefore, in the layered model, cohesion varies significantly among different layers as a function of moisture content, while the internal friction angle remains nearly constant. Finally, the calculated layer-specific strength parameters were assigned to the corresponding zones in the FLAC3D model, enabling the numerical simulation to explicitly account for rainfall-induced strength degradation and establishing a clear and traceable data transfer process between the COMSOL seepage analysis and the FLAC3D stability analysis.

Boundary Conditions (Static and Dynamic Stages):

Static Stage: Fixed in all displacement directions (fix x y z).

Dynamic Stage: To simulate seismic wave propagation from the bedrock upward, the bottom boundary was changed to a viscous boundary to absorb reflected waves and reduce boundary effects. The seismic acceleration time history (baseline-corrected and filtered Kobe wave, peak 0.2 g) was applied as a stress on this boundary as the input seismic load.

Model Lateral Boundaries:

Static Stage: Normal displacement was constrained, while tangential sliding was allowed (simulating “roller support”) to approximate the lateral constraint conditions of a semi-infinite space.

Dynamic Stage: Similarly set as free-field boundaries or viscous boundaries to simulate seismic wave propagation in an infinite domain, allowing the model to undergo realistic deformation under seismic action while absorbing waves propagating outward from inside the model.

Model Top Boundary: Always a free surface, unconstrained, allowing free deformation.

The soil was modeled using the Mohr-Coulomb elastoplastic constitutive model. Rayleigh damping was applied in the dynamic analysis to simulate the energy dissipation characteristics of the soil, with the damping ratio set based on empirical values for the relevant soil type.

Fig. 13 presents the distribution of vertical stress within the saline soil subgrade model under initial geostatic equilibrium conditions prior to the application of seismic loading. The stress field was generated through static equilibrium analysis in FLAC3D, considering only gravitational forces. This figure establishes the initial stress baseline of the subgrade slope, which is essential for subsequent seismic response analysis. The absence of external loads (e.g., traffic, earthquake) ensures that the stress state is solely due to soil self-weight, providing a reference for evaluating stress changes induced by seismic action. Fig. 14 illustrates the total displacement field of the subgrade model resulting from the initial geostatic equilibrium process. The displacements are induced by soil self-weight during the static balancing phase and represent elastic settlement before dynamic loading. This figure provides the initial displacement baseline for the subgrade. In subsequent seismic simulations, all dynamic displacements are computed relative to this initial state. The displacement field demonstrates that the model is properly initialized and ready for transient dynamic analysis.

Figure 13: Stress cloud diagram under initial state.

Figure 14: Displacement cloud diagram under initial state.

6.3 Stability Analysis of Subgrade Slopes under Runoff and Seismic Influence

Seismic load, also known as seismic force, is the combined effect of inertial force, earth pressure, and water pressure on a structure due to an earthquake. Seismic waves propagate through the soil layer in the form of longitudinal and transverse waves. Longitudinal waves cause vertical vibrations in structures, while transverse waves cause horizontal oscillations. Extensive seismic damage phenomena indicate that horizontal seismic action is the main cause of slope failure. Therefore, this paper only considers the influence of transverse waves on slope stability. The selected seismic wave is the Kobe wave, and the selected time period is a representative segment from 4 to 24 s, suitable for seismic time history analysis. Since the basic seismic peak acceleration of the study area is 0.1 to 0.2 g, the peak acceleration is adjusted to 0.2 g. The acceleration time history curve is shown in Fig. 15, with the peak acceleration of −0.2 g occurring at 8.52 s.

Figure 15: Peak adjusted Han Shenbo.

Since the input is an acceleration time history curve of the seismic wave, the final displacement obtained by integrating the acceleration in the software is not zero. If the seismic wave acceleration is directly input into the calculation, the model will exhibit continued velocity and residual displacement at the bottom at the end of the dynamic calculation. Therefore, baseline correction is required before dynamic analysis. The principle of baseline correction is to add a low-frequency polynomial or periodic function waveform (a low-frequency sine wave in the software) to the original acceleration time history to ensure that the final velocity and displacement are zero, as shown in Fig. 16.

Figure 16: Baseline correction.

The purpose of filtering is to remove the high-frequency components of the seismic wave, as the maximum frequency of the seismic wave significantly affects the grid size. The higher the maximum frequency of the seismic wave, the smaller the grid size required to meet accuracy conditions. Therefore, to avoid excessively high grid size accuracy requirements, the selected seismic wave needs to be filtered to save computation time. The relationship between the maximum grid size and the shortest wavelength of the input seismic wave is expressed as Eqs. (6) and (7) [34]. The filtering process is shown in Fig. 17.

Δl≤(18:110)λ(6)

where, Δl is the maximum size for meshing; λ is the shortest wavelength of the input seismic wave.

Figure 17: Filtering processing.

According to elastic wave propagation theory, the propagation speed of transverse waves in a medium can be expressed by the following formula [35]:

Cs=Gρ(7)

where, G is shear modulus of medium; ρ is medium density.

Analysis of Subgrade Surface Settlement and Displacement

Numerical simulation analysis was conducted on the subgrade under seismic action at different rainfall intensities. The settlement cloud diagrams of the subgrade slope at rainfall intensities of 1 and 5 mm/h are shown in Figs. 18 and 19. It can be observed that as the rainfall intensity increases, the settlement displacement of the subgrade slope gradually increases, with the maximum settlement displacement concentrated at the subgrade slope.

Figure 18: Displacement cloud diagram of subgrade slope under seismic action with rainfall intensity of 1 mm/h.

Figure 19: Displacement cloud diagram of subgrade slope under seismic action with rainfall intensity of 5 mm/h.

From the Figs. 18 and 19, it can be seen that the maximum displacement of the subgrade slope changes consistently with the rainfall intensity. Under seismic action, the subgrade slope has the smallest settlement at a rainfall intensity of 1 mm/h, with a settlement value of 18.16 mm. At a rainfall intensity of 5 mm/h, the settlement reaches its peak at 25.60 mm, a decrease of 40.97%. Since this paper only considers the reduction of subgrade strength due to surface runoff, the seismic action merely promotes the destruction of the subgrade, exacerbating the degree of subgrade settlement. This makes the pattern of subgrade settlement, which would normally require long-term cumulative damage to become apparent, manifest clearly in a short period. However, it does not significantly affect the pattern of subgrade settlement deformation. Therefore, the following recommendations are proposed for subgrade protection: To ensure the driving comfort and safety of the road, special attention should be paid to seismic warnings in July when rainfall is abundant. Long-term monitoring of subgrade settlement should be conducted, with a focus on observing slope settlement changes in July, and necessary reinforcement and prevention measures should be taken when needed.

7  Discussion and Limitations

The present study relies on several secondary equations to describe key physical processes, including unsaturated seepage behavior, moisture-dependent strength degradation, and dynamic response of the subgrade slope. While these formulations provide an efficient and practical framework for numerical analysis, their inherent assumptions and limitations should be carefully acknowledged.

In the study, several secondary governing equations were employed to describe unsaturated seepage, moisture-dependent strength degradation, and the dynamic response of the subgrade slope. Although these formulations are widely used in engineering practice due to their efficiency and practicality, their underlying assumptions and applicability require clarification. First, the unsaturated seepage analysis is based on the Richards equation, which assumes continuous pore-water connectivity and does not explicitly account for hysteresis effects in the soil–water characteristic curve. For saline soils with complex structures and compositions, localized preferential flow paths and salt crystallization processes may influence moisture migration, which are not fully captured by this formulation. In addition, the dynamic analysis adopts equivalent elastic soil parameters and simplified wave propagation assumptions. While such simplifications are commonly used to achieve a balance between computational efficiency and physical realism, they may underestimate the nonlinear soil behavior under strong seismic excitation, thereby limiting the accuracy of the predicted dynamic response under extreme loading conditions. Therefore, the results of this study should be interpreted primarily as indicative of relative trends and comparative effects rather than precise predictions of real-field behavior. Future research may incorporate hysteretic seepage processes, fully coupled hydro-mechanical models, and nonlinear dynamic constitutive relationships to further enhance the reliability and resolution of stability assessments for saline soil subgrade slopes.

8  Conclusion

This study established a layered moisture content–strength parameter equivalence model to investigate the stability of saline soil subgrade slopes under the combined effects of surface runoff and seismic action. The transient moisture fields induced by rainfall infiltration were first obtained using COMSOL Multiphysics, and the corresponding layer-specific strength parameters were subsequently incorporated into FLAC3D to evaluate the seismic response of the subgrade slope. Based on the numerical results, the following conclusions can be drawn.

(1) Variations in moisture content have a pronounced influence on the shear strength of saline soils, with cohesion being significantly more sensitive to moisture changes than the internal friction angle. An increase in moisture content leads to substantial cohesion degradation, which plays a dominant role in reducing slope stability under rainfall and seismic loading conditions.

(2) Rainfall intensity primarily affects the moisture content of the shallow soil layers of the subgrade slope. As rainfall intensity increases, shallow moisture accumulation and surface runoff become more pronounced, resulting in a continuous decrease in the slope safety factor. Under seismic action, this moisture-induced strength degradation is further amplified, leading to increased surface settlement and deformation of the subgrade slope.

(3) From an engineering perspective, the results indicate that saline soil subgrade slopes in arid and semi-arid regions are particularly vulnerable during periods of concentrated rainfall combined with seismic activity. Engineering design should therefore emphasize effective surface drainage systems to reduce shallow moisture accumulation and consider the moisture-dependent degradation of soil strength in seismic stability assessments. In regions with seasonal rainfall concentration, enhanced monitoring and reinforcement measures are recommended during high-risk periods.

(4) Several limitations of this study should be acknowledged. The rainfall process was simplified as a single short-duration infiltration event, and the cumulative effects of multiple rainfall cycles and long-term moisture migration were not considered. In addition, the seismic loading was represented by a single input motion, and the potential influence of different seismic wave characteristics was not examined.

(5) Future research should focus on incorporating long-term rainfall–seepage processes and repeated seismic loading into the coupled numerical framework. Further laboratory and field investigations are also required to better quantify the coupled effects of moisture, salinity, and cyclic loading on the mechanical behavior of saline soils, thereby improving the reliability of stability evaluations and engineering design methods for saline soil subgrades.

Acknowledgement: The authors wish to express their sincere appreciation to the relevant research institutions and laboratories for providing the technical support and resources that made this study possible.

Funding Statement: The research was supported the Gansu Joint Research Fund Project (25YFFA023 and 25JRRA1119).

Author Contributions: Bingbing Han contributed to the conceptualization of the study, methodology development, numerical modeling, and preparation of the original manuscript draft. Ruidong Li participated in laboratory testing, data curation, parameter determination, and assisted in numerical simulation and result analysis. Yu Zhang was responsible for the overall study design, supervision of the research, interpretation of results, and critical revision of the manuscript. Hui Zhang contributed to software implementation using COMSOL Multiphysics and FLAC3D, simulation analysis, and visualization of results. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: All data, models, and codes that support the findings of this study are available from the corresponding author upon reasonable request.

Ethics Approval: This study does not involve any human participants, human data, animals, or animal-derived materials. All analyses were conducted using numerical simulation methods and laboratory-derived parameters reported in the literature. Therefore, ethical approval was not required for this research.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. GB 50021-2001. Code for investigation of geotechnical engineering. Beijing, China: China Standards Press; 2001. [Google Scholar]

2. Zhang Y, Tian RZ, Wang TL. Study on salt expansion mechanism of subgrade-culvert transition section in saline soils and cold regions. Cold Reg Sci Technol. 2023;205(6):103701. doi:10.1016/j.coldregions.2022.103701. [Google Scholar] [CrossRef]

3. Alonso E, Gens A, Lloret A, Delahaye C. Effect of rain infiltration on the stability of slopes. In: Proceedings of the First International Conference on Unsaturated Soils; 1995 Sep 6–8; Paris, France. p. 241–9. [Google Scholar]

4. Rahardjo H, Li XW, Toll DG, Leong EC. The effect of antecedent rainfall on slope stability. Geotech Geol Eng. 2001;19(3):371–99. doi:10.1023/A:1013129725263. [Google Scholar] [CrossRef]

5. Dou HQ, Han TC, Gong XN, Zhang J. Probabilistic slope stability analysis considering the variability of hydraulic conductivity under rainfall infiltration–redistribution conditions. Eng Geol. 2014;183:1–13. doi:10.1016/j.enggeo.2014.09.005. [Google Scholar] [CrossRef]

6. Gavin K, Xue J. A simple method to analyze infiltration into unsaturated soil slopes. Comput Geotech. 2008;35(2):223–30. doi:10.1016/j.compgeo.2007.04.002. [Google Scholar] [CrossRef]

7. Sajjan AK, Gyasi-Agyei Y, Sharma RH. Rainfall-runoff modelling of railway embankment steep slopes. Hydrol Sci J. 2013;58(5):1162–76. doi:10.1080/02626667.2013.802323. [Google Scholar] [CrossRef]

8. Song X, Li H. Analysis of rainfall infiltration stability of high liquid limit red clay roadbed in Africa. Comput Intell Neurosci. 2022;2022:8222075. doi:10.1155/2022/8222075. [Google Scholar] [PubMed] [CrossRef]

9. Guo Z, Zhao Z. Numerical analysis of an expansive subgrade slope subjected to rainfall infiltration. Bull Eng Geol Environ. 2021;80(7):5481–91. doi:10.1007/s10064-021-02274-7. [Google Scholar] [CrossRef]

10. Fan SY, Song ZP, Zhang YW, Liu NF. Case study of the effect of rainfall infiltration on a tunnel underlying the roadbed slope with weak inter-layer. KSCE J Civ Eng. 2020;24(5):1607–19. doi:10.1007/s12205-020-1165-0. [Google Scholar] [CrossRef]

11. Liu J, Yang C, Gan J, Liu Y, Wei L, Xie Q. Stability analysis of road embankment slope subjected to rainfall considering runoff-unsaturated seepage and unsaturated fluid-solid coupling. Int J Civ Eng. 2017;15(6):865–76. doi:10.1007/s40999-017-0194-7. [Google Scholar] [CrossRef]

12. Gao J, Jiang D, Du W. Laboratory study on the water-heat behavior and freeze-thaw deformation of a composite embankment considering rainfall in seasonally frozen regions. Math Probl Eng. 2022;2022:9706867. doi:10.1155/2022/9706867. [Google Scholar] [CrossRef]

13. Kanji S, Das S. Assessing predictability of spatio-temporal soil-moisture heterogeneity with satellite-data using GeoDetector and ensemble machine learning models considering environmental and climatic factors. J Hydrol. 2025;663:134312. doi:10.1016/j.jhydrol.2025.134312. [Google Scholar] [CrossRef]

14. Kanji S, Das S. Exploring the morpho-tectonic nature, hydrological and physical characteristics of a watershed and prioritizing sub-watersheds surface runoff potentialities by integrating MCDM and ensemble machine learning models. J Environ Manag. 2025;386:125772. doi:10.1016/j.jenvman.2025.125772. [Google Scholar] [PubMed] [CrossRef]

15. Rajabian A, Vahedifard F. Seismic stability analysis of anchored concave slopes. Transp Geotech. 2025;51:101505. doi:10.1016/j.trgeo.2025.101505. [Google Scholar] [CrossRef]

16. Liao Z, Jiang Z, Ma K, Gao Z, Ke H, Wei A. Investigation on cumulative response evolution and stability assessment of rock slope under mainshock and aftershocks. Comput Geotech. 2025;180:107092. doi:10.1016/j.compgeo.2025.107092. [Google Scholar] [CrossRef]

17. Wu X, Cui Z, Yang B, Lian B. Seismic stability of anchored slopes considering tension-shear coupling effect by modified pseudo-dynamic approach. Soil Dyn Earthq Eng. 2025;190:109206. doi:10.1016/j.soildyn.2024.109206. [Google Scholar] [CrossRef]

18. Wang R, Yan H, Yao J, Li Z. Stability analysis of slopes under seismic action with asynchronous discounting of strength parameters. Appl Sci. 2025;15(1):169. doi:10.3390/app15010169. [Google Scholar] [CrossRef]

19. Zhou Y, Zhao F, Shi Z. Dynamic response mechanism of bedding slopes with alternatively distributed soft and hard rock layers under different seismic excitation directions: insights from numerical simulations. Materials. 2024;17(23):5939. doi:10.3390/ma17235939. [Google Scholar] [PubMed] [CrossRef]

20. Liu C, Li Y, Wang L. Stability of slopes in partially saturated soils: incorporating the combined effects of seismic forces and pore water pressure. Soil Dyn Earthq Eng. 2024;187:108996. doi:10.1016/j.soildyn.2024.108996. [Google Scholar] [CrossRef]

21. Konovalov AV, Gensiorovsky YV. Assessment of slope stability taking into account the uncertainty of standard seismic impacts and material parameters of the site. Russ J Pac Geol. 2024;18(5):590–600. doi:10.1134/S1819714024700349. [Google Scholar] [CrossRef]

22. Leshchinsky D, San KC. Pseudostatic seismic stability of slopes: design Charts. J Geotech Engrg. 1994;120(9):1514–32. doi:10.1061/(asce)0733-9410(1994)120:9(1514). [Google Scholar] [CrossRef]

23. Hall L. Simulations and analyses of train-induced ground vibrations in finite element models. Soil Dyn Earthq Eng. 2003;23(5):403–13. doi:10.1016/s0267-7261(02)00209-9. [Google Scholar] [CrossRef]

24. Massey C, Della Pasqua F, Holden C, Kaiser A, Richards L, Wartman J, et al. Rock slope response to strong earthquake shaking. Landslides. 2017;14(1):249–68. doi:10.1007/s10346-016-0684-8. [Google Scholar] [CrossRef]

25. Giri D, Sengupta A. Dynamic behavior of small scale nailed soil slopes. Geotech Geol Eng. 2009;27(6):687–98. doi:10.1007/s10706-009-9268-x. [Google Scholar] [CrossRef]

26. Feng ZY, Lu YR, Shen ZR. A numerical simulation of seismic signals of coseismic landslides. Eng Geol. 2021;289:106191. doi:10.1016/j.enggeo.2021.106191. [Google Scholar] [CrossRef]

27. Cheng J, Zhang Y, Ma Y, Chen X, An N. The effect of seismic action on stability of saline soil subgrade in cold region based on isothermal stratification method. Earthq Res Adv. 2024;4(3):100271. doi:10.1016/j.eqrea.2023.100271. [Google Scholar] [CrossRef]

28. GB 18306-2015. China seismic ground motion parameter zoning map. Beijing, China: China Standards Press; 2015. [Google Scholar]

29. Zhang Y, Luo Y, Xu AH, Fang J. Stability of high-salinity saline soil subgrade slope under the influence of moisture content and shear strength relationship. J Chang Univ Nat Sci Ed. 2020;40(3):22–32. doi:10.1088/1755-1315/384/1/012072. [Google Scholar] [CrossRef]

30. Li ZM. Research on coupling mechanism of hydro-thermo-mechanical field in frozen soil and engineering application [dissertation]. Harbin, China: Harbin Institute of Technology; 2017. [Google Scholar]

31. Xu W, Huang X, Huang J, Yang Z. Structural analysis of backfill highway subgrade on the lower bearing capacity foundation using the finite element method. Adv Civ Eng. 2021;2021(1):1690168. doi:10.1155/2021/1690168. [Google Scholar] [CrossRef]

32. Zhang Y, Yang Z, Liu J, Fang J. Impact of cooling on shear strength of high salinity soils. Cold Reg Sci Technol. 2017;141(4):122–30. doi:10.1016/j.coldregions.2017.06.005. [Google Scholar] [CrossRef]

33. Liu HY. Evaluation of shear strength and slope stability of unsaturated soil in geological disaster-prone areas with different soil types [dissertation]. Lanzhou, China: Lanzhou University; 2023. [Google Scholar]

34. Kuhlemeyer RL, Lysmer J. Finite element method accuracy for wave propagation problems. J Soil Mech Found Div. 1973;99(5):421–7. doi:10.1061/jsfeaq.0001885. [Google Scholar] [CrossRef]

35. Chen YM, Xu DP. FLAC/FLAC3D foundation and engineering examples. Beijing, China: China Water & Power Press; 2013. [Google Scholar]