A Case Study of Using Numerical Analysis to Assess Slope Stability of National Highways in Northern Taiwan

1. Introduction

Many nations have constructed substantial road networks and residential communities in hilly and mountainous terrain as a result of global population increase and socioeconomic development. But occasionally, landslides happen, and they can be extremely dangerous for people and property, particularly in places with steep terrain and a lot of rainfall [1,2,3]. Thus, geotechnical researchers have focused increasingly on slope stability assessment [4,5,6,7,8,9,10,11,12,13,14,15,16].

In mountainous regions, landslides rank among the most dangerous natural calamities [17,18]. Taiwan's landscape is primarily made up of hills and mountains, which make up almost two-thirds of the island's overall area. As a result, safety in the event of landslides is a major concern for the appropriate government bodies. Since the construction of National Freeway 1, the demand for motorways has grown. The freeway segment between Keelung and Xizhi was constructed 46 years ago, and shotcrete has been used to cover the bedrock slope's surface. The shotcrete covering has been repaired and maintained, but aging has led to cracks and plant growth on the surface. It is worth noting that on National Freeway 5, several slopes are classified as high-risk according to the technical specifications of the Freeway Maintenance Manual of the Freeway Bureau of the Ministry of Transportation and Communications of the Republic of China in Taiwan.

The southbound exits at Wudu on National Freeways 1 and 5 with high-risk slopes were chosen as case studies in light of the aforementioned. Numerous landslides have occurred close to these locations. Slope displacement at normal and high groundwater levels was simulated using the PLAXIS 2D CE finite element analysis tool, and variations in acceleration, speed, and displacement in slope models were examined using pseudo-static analysis in seismic circumstances. In the future, the numerical simulation results might be helpful for preventing the collapse of roadway slopes when severe rains or other conditions impair their stability.

2. Literature review

2.1. Slope Stability

Under specific geological background circumstances, slope instability typically arises as a result of the excitation of initiating events like intense earthquakes and heavy rains [19]. A hillslope's stability can be affected by a variety of causes. Geological material, geological structure, topographic and environmental conditions, and engineering considerations are the four main causes of landslides, according to Hung [20].

• Geological material: The slope stratum is mainly composed of single or multiple geological materials. The cementation, particle size, and composition of the geological materials will directly affect the stability of the slope.
• Geological structure: In geological structures, the orientation of weak planes (such as strike and dip angle) and the type of slope (such as dip slope, escarpment, and oblique slope) are crucial characteristics. Slope stability is compromised by other geological features that split slopes into discontinuous or fragmented rocks, such as bedding surfaces, joints, folds, faults, and other discontinuities. This reduces the strength of the rock mass and increases weathering.
• Environmental and topographic factors: Topography includes the slope's general form and rough surface terrain. While slope terrain has its own height and gradient, surface terrain relates to the characteristics of the geological structure. Groundwater, earthquakes, and rainfall are examples of environmental elements that affect slope stability.
• Engineering aspects: Slope stability is affected directly or indirectly by human factors such as highways, tunnel excavation, blasting, overdevelopment of slopes, bad building or site selection, poor slope drainage systems, inadequate maintenance of protective structures, etc.

In Taiwan, accidents involving slopes are common in building projects. Slope patrol, slope monitoring, and ground anchor inspection are used to keep an eye on stability in order to prevent collapse due to the sudden failure of rock or soil mass on a slope. To increase the facilities' service life, slope safety evaluation and categorization are also done. To put it briefly, these facilities' preventive maintenance lowers the chance of failure and averts catastrophes. Slopes are divided into four classes—A, B, C, and D—per the Freeway Maintenance Manual. The following is a description of the classification criteria and remedies:

• Class A slope: Instability is obviously present. Increased patrols and observation are necessary, as is prompt intervention.
• Class B slope: Due to the possible indication of instability, additional patrols, monitoring, maintenance, reinforcing, and remediation are needed.
• Class C slope: It needs regular patrol or maintenance, with monitoring as necessary, even if there are no obvious signs of instability.
• Class D slope: Patrolling is still necessary even though it is steady.

2.2. Prediction of Slope Stability

Predicting slope stability has been the subject of numerous studies over the last few decades. These research' methodologies can be broadly categorized into three groups: finite element technique (FEM), limit equilibrium method (LEM), and limit analysis method (LAM) [21]. Due to its ease of use, the LEM is one of the most popular approaches among them [22]. However, several idealized assumptions about the LEM make its response imprecise. Because it puts strict restrictions on limit state solutions [23] and offers a straightforward calculation procedure, the LAM is commonly used in slope stability [24]. However, when dealing with complicated slope structures and variable soil deposition, it can be difficult to create theoretical formulations of upper or lower bounds [21].

Compared to the LEM and LAM, the FEM appears to be a more reliable and adaptable processing technique for slope performance modeling and evaluation [21,25,26,27]. Kafle et al. [7] constructed a finite element model to simulate the Bianjiazhai landslide near the Suofengying reservoir at the Wu River in Guizhou Province, China. Their results showed that sudden reservoir level fluctuations had a critical effect on the stability of the slope. In addition, a strength reduction analysis confirmed that the deformation behavior of the slope was mainly a shear-slip phenomenon, and the associated shear zones were identified. In order to accurately and efficiently simulate in situ stress prior to seismic loading, Kontoe et al. [28] used pseudo-static finite element analysis. They also investigated how sensitive the analysis findings were to the finite element mesh level.

The main metric used to assess slope stability is the safety factor [10], which is the ratio of the available shear strength to the equilibrium shear stress along a specific failure surface [4]. Assumptions on the form and location of the failure surface are not required when utilizing the FEM for a slope stability analysis. According to the FEM, failure occurs when the soil's shear strength is not strong enough to support the sliding force of the slope, and it is feasible to model how displacement and slope failure will develop. The FEM often employs the soil shear strength reduction technique to calculate the slope's safety factor [29]. During the slope failure process, soil cohesion c and internal friction angle φ continue to decrease. According to the literature, no matter how the slope is formed, when it reaches the failure state, the reduced cohesion tends to have the same value. Accordingly, Zienkiewicz et al. [30] reduced the strength parameters (c, tanφ) of a slope step-by-step through a reduction factor of 1/N, and the reductions for the two parameters were equal. When reduced to the failure state, the N value at this time can be defined as the safety factor of the slope.

Using the region close to the top of the neighboring slope toe as the chosen node, Giam and Donald [31] calculated the slope's overall safety factor by examining the link between the node displacement and the change in soil strength characteristics (cohesion and friction angle). Ugai [32] determined the system's overall safety factor iteratively by adopting a minimal system safety factor value for the material's strength (c and φ). To determine whether the system is unstable, they then performed finite element analysis using the converted material strength. Brinkgreve and Bakker [33] used the PLAXIS finite element formula they developed to calculate the slope's safety factor and the corresponding failure sliding surface position in addition to the shear strength reduction technique. The idea of shear strength reduction was used by Matsui and San [34] to assess slope safety factors and the emergence of possible slip surface failure. Additionally, the shear strain in the failure process that extends from the slope's base to its summit (i.e., total shear failure) was classified as slope failure. Ugai and Leshchinsky [35] extended the two-dimensional problem to three dimensions. They used 20-node isoparametric solid elements and reduced integration with 8 Gaussian integration points to conduct a three-dimensional vertical excavation slope simulation, taking into account the influence of pseudo-static seismic forces. The analysis results are compared with the results of the limit equilibrium method to verify their rationality. Griffiths and Lane [36] used the FEM to apply shear strength reduction to the analysis of slope stability, taking into account the effects of weak interlayers in the soil and groundwater. The occurrence of slope failure was determined by the non-convergence of the solution, because slope failure is often accompanied by a rapid increase in displacement. The position and shape of the critical slip surface were represented by a deformed mesh diagram and node displacement vector diagram. In addition, the analysis results of the safety factor were compared with those of the limit equilibrium method. Based on the incremental method of elastic–plastic mechanics and the bilinear projection operator, Sun et al. [29] combined the strength reduction method with the φ-ν inequality and proposed a virtual element strength reduction method for slope stability analysis. The deformation of homogeneous and heterogeneous slopes under different strength reduction coefficients was solved, and the grid dependence of this method was discussed. Furthermore, numerical examples were used to verify the correctness and effectiveness of the proposed method.

2.3. Pseudo-Static Analysis of Soil Inertial Force

One of the primary dynamic loads frequently observed on soil or rock slopes is seismic force [19]. Rock or soil slopes are subject to enormous dynamic stress, which has a significant impact on their stability. The pseudo-static analysis based on force equilibrium assumes that seismic force is a static inertial force operating laterally on the structure while assessing structural stability. Throughout the analysis, the seismic coefficient was progressively raised. To ascertain the seismic load that the structure can withstand as it approaches critical stability, a number of stability analyses were conducted.

The limit equilibrium approach, also known as the Mononobe–Okabe (M-O) method, was created in Japan in the 1920s by Okabe [37] and Mononobe and Matsuo [38]. It serves as the foundation for the most widely used cantilever and gravity structure analysis. For pseudo-static analysis, the M-O approach is also frequently employed. It calculates the dynamic soil pressure coefficient based on force equilibrium after transforming the seismic load into the inertial force of soil mass. Dynamic soil pressure, or active soil pressure under dynamic conditions, is depicted in equation (1). Figure 1a displays the stresses produced under dynamic conditions, while Figure 1b displays the stress path produced by the soil mass's inertial force under static conditions.

$$P_{AE}={\displaystyle \frac{1}{2}}\left(1\pm K_V\right)K_{\mathit{AE}}\gamma H^2$$

(1)

In Equation (1), K_V is the vertical seismic coefficient, K_AE is the dynamic soil pressure coefficient, $\gamma$ is the unit weight of soil, and H is the wall height. Various modifications of the M-O method have been introduced. Seed and Whitman [39] separated dynamic soil pressure (P_AE) into static soil pressure (P_A) and seismically induced soil pressure (ΔP_dyn), as shown in Equation (2):

$$P_{AE}=P_A+\Delta P_{dyn}$$

(2)

3. Overview of the Case Studied

For this investigation, three locations were selected. Site 1 may be found at the Freeway 1 Wudu exit (Figure 2), while Site 2 is situated at Freeway 5 with a 7K+082 mileage (Figure 3), and Site 3 is situated at Freeway 5 with a 7K+410. The geological data and the overview of the slope protection project for each site are described below, and they act as a guide for the input values utilized in the subsequent numerical analysis.

According to the geological maps supplied by Taiwan's Central Geological Survey, the stratum at Site 1 is the Nangang Formation, which is made up of thick layers of dark shale and slate-gray tuffaceous sandstone. The slope is protected by RC curtain walls, free-styled grill beams, and grouted anchor bars. The grouted anchor bars are composed of three-meter-long #8 steel bars spaced three meters apart horizontally.

The Nangang Formation, the stratum at Site 2, is interlayered with sandstone and shale and is located 4.9 meters below the colluvium's ground surface. A total of 26 layers of 30 tons of prestressed ground anchors are used to protect the slope. The current ground anchors consist of a fourth-stage, 15-layer ground anchor; a free section of five meters; an anchor section of five meters; and a horizontal spacing of three meters. After that, eleven levels of ground anchors were placed, spaced three meters apart horizontally, with ten meters for the free section and five meters for the anchor section.

The Nangang Formation is the layer found at location 3. According to the results of the drilling test, there is a layer of colluvium four meters below ground level, with powdery sandstone underneath. Ten layers of thirty tons of prestressed ground anchors are used to protect the slope. The current ground anchors are seven layers in length, with a free portion of fourteen meters, an anchor section of eleven meters, and a horizontal spacing of two and a half meters. Three layers of ground anchors were then placed, with an anchor section of 11 meters, a free section measuring 14 meters, and a horizontal spacing of two and a half meters.

4. Numerical Analysis Software and Methods Used

The PLAXIS 2D CE finite element analysis program, developed by Plaxis B.V., was used for numerical analysis and calculations in this study. The program was first developed in 1987 by Delft University of Technology in the Netherlands. The company was founded in 1993, and the first Windows-based version, PLAXIS 2D, was released in 1998. After years of research and promotion, a three-dimensional (3D) version was developed in 2001, and then PLAXIS 2D AE was developed in 2014 and updated to PLAXIS 2D CE in 2019.

The Mohr–Coulomb model is a common basic rock mass material model for numerical analysis in geological engineering. Its stress–strain relationship covers linear elasticity and complete plasticity. The linear elasticity part follows Hooke’s law and the required stiffness parameters are Young’s modulus E, and Poisson’s ratio ν, whereas the plasticity part follows the Mohr–Coulomb failure criterion, of which the maximum principal stress (${\sigma }_1^{\text{'}}$) vs. the minimum principal stress (${\sigma }_3^{\text{'}}$) is expressed in Equation (3), and observes the principle of non-associated plasticity. The strength parameters needed are cohesion c, angle of internal friction

, and angle of dilatancy ψ. The parameters required for rock masses in the Mohr–Coulomb model are listed in Table 1.

$${\sigma }_1^{\text{'}}={\displaystyle \frac{2c\mathrm{c}\mathrm{o}\mathrm{s}}{1-\mathrm{s}\mathrm{i}\mathrm{n}}}+{\displaystyle \frac{1+\mathrm{s}\mathrm{i}\mathrm{n}}{1-\mathrm{s}\mathrm{i}\mathrm{n}}}{\sigma }_3^{\text{'}}$$

(3)

5. Numerical Analysis Results and Discussion

This study focused on two aspects: the displacement of the slope at normal and high groundwater levels and the displacement of the slope model while it is subjected to seismic loading during pseudo-static analysis.

5.1. Normal Groundwater Level and High Groundwater Level

The "normal" groundwater level employed in the slope stability analysis was the average groundwater level as established by survey monitoring wells. The high groundwater level was estimated to be two-thirds of the total level.

5.2. Seismic Loading in the Pseudo-Static Analysis

The seismic loading utilized in the pseudo-static analysis was taken from the Seismic Design Specifications for Highway Bridges published by the Ministry of Transportation and Communications of the Republic of China in Taiwan. The case was located on a type 1 rock bed in the Shiding District of New Taipei City; hence the seismic coefficient used in the pseudo-static analysis was as follows:

$$k_h=0.5 \text{×} 0.4S_{DS}=0.5 \text{×} 0.4 \text{×} 0.7 g = 0.14 g$$

(4)

$$k_v=0.5\left(k_h\right)=0.07 g$$

(5)

where $S_{DS}$ is the designed short-period horizontal seismic acceleration of the site, $k_h$ is horizontal seismic acceleration, and $k_v$ is vertical seismic acceleration.

5.3. PLAXIS Analysis

PLAXIS 2D was used to build slope simulation models, with the Mohr–Coulomb model used for soil layers. The soil parameters used in the model were unsaturated unit weight ${\gamma }_{\mathrm{u}\mathrm{n}\mathrm{s}\mathrm{a}\mathrm{t}}$, saturated unit weight ${\gamma }_{\mathrm{s}\mathrm{a}\mathrm{t}}$, Young’s modulus E, Poisson’s ratio ν, cohesion c, internal friction angle φ, and dilatancy angle ψ. The model consisted of grouted anchor bars, free-styled grill beams, RC curtain walls, free and fixed lengths of ground anchor, and a typical retaining wall.

In this study, the corresponding simulation elements were set according to the site structure. The slope protection structure included grouted anchor bars, lattice beams, RC curtain walls, stiffened retaining walls, and ground anchors. Among them, the grouted anchor bars, the stiffened grid in the stiffened retaining wall, and the anchor section of the ground anchor were simulated by embedded beam elements; the lattice beam, RC curtain wall, and RC panel of the stiffened retaining wall were simulated; the road surface was simulated by the plate element; and the free segment of the ground anchor was simulated by the node-to-node anchor element.

5.3.1. Case Simulation of Site 1

Table 1 displays the soil properties for Site 1. Table 2 and Table 3 display the material parameters of Site 1's slope protection structures. Figure 4 displays Site 1's analysis model. The blue line represents simulated free-form lattice beams, the yellow line represents simulated RC curtain walls, the blue arrows represent the pavement load, and the pink parts represent simulated grouted anchor bars. From top to bottom, the soil layer was separated into two layers: sandstone and worn rock.

Table 1. Soil parameters of Site 1.

Soil layer	Parameters
	${\mathit{\gamma }}_{\mathbf{u}\mathbf{n}\mathbf{s}\mathbf{a}\mathbf{t}}$ (kN/m3)	${\mathit{\gamma }}_{\mathbf{s}\mathbf{a}\mathbf{t}}$ (kN/m3)	E (kN/m2)	ν	c (kN/m2)	φ ( ° )	ψ ( ° )
Weathered rock layer	22.0	22.5	2.5 × 104	0.25	10	28	0
Sandstone layer	25.2	25.7	4 × 106	0.25	150	26	0

Table 1. Soil parameters of Site 1.

Soil layer	Parameters
	${\mathit{\gamma }}_{\mathbf{u}\mathbf{n}\mathbf{s}\mathbf{a}\mathbf{t}}$ (kN/m3)	${\mathit{\gamma }}_{\mathbf{s}\mathbf{a}\mathbf{t}}$ (kN/m3)	E (kN/m2)	ν	c (kN/m2)	φ ( ° )	ψ ( ° )
Weathered rock layer	22.0	22.5	2.5 × 104	0.25	10	28	0
Sandstone layer	25.2	25.7	4 × 106	0.25	150	26	0

Table 2. Material parameters of grouted anchor bars of Site 1.

Simulation elements	Young’s modulus (kN/m2)	Unit weight (kN/m3)	Diameter (cm)	Horizontal spacing (m)	Front end side friction resistance (kN/m)	Rear end side friction resistance (kN/m)
Embedded beam row element	2.5 × 107	10	25	3	509	509

Table 2. Material parameters of grouted anchor bars of Site 1.

Simulation elements	Young’s modulus (kN/m2)	Unit weight (kN/m3)	Diameter (cm)	Horizontal spacing (m)	Front end side friction resistance (kN/m)	Rear end side friction resistance (kN/m)
Embedded beam row element	2.5 × 107	10	25	3	509	509

Table 3. Material parameters of free-form lattice beams and RC curtain walls of Site 1.

Simulation elements	Material type	Axial stiffness (kN/m3)	Flexural stiffness (kNm2/m)	Weight (kN/m/m)	ν
Plate element	Elasticity	7.53 × 106	$\mathrm{56,475}$	7.06	0.17

Table 3. Material parameters of free-form lattice beams and RC curtain walls of Site 1.

Simulation elements	Material type	Axial stiffness (kN/m3)	Flexural stiffness (kNm2/m)	Weight (kN/m/m)	ν
Plate element	Elasticity	7.53 × 106	$\mathrm{56,475}$	7.06	0.17

Figure 4. PLAXIS 2D model of Site 1.

Figure 4. PLAXIS 2D model of Site 1.

Shading diagrams of slope displacement in various conditions are displayed in Figure 5, Figure 6 and Figure 7. The force acting on the slope was centered on the slope's sliding surface because Site 1 had a normal groundwater level, as seen in Figure 5. In contrast, Figure 6 illustrates how the force acting on the slope went up to the top of the slope at high groundwater levels. The pseudo-static analysis of Site 1 indicated that the force acting on the slope extended from the top to the sliding surface when seismic loading was added, as seen in Figure 7. Table 4 displays the displacement of Site 1 in each scenario.

5.3.2. Case Simulation of Site 2

Table 5 displays the soil properties for Site 2. Table 6 and Table 7 display the material parameters of Site 2's slope protection structures. Figure 8 displays Site 2's analysis model. The road load is indicated by the blue arrow, the free segment of the ground anchor is indicated by the black line, and the anchor segment is indicated by the red line. From top to bottom, the soil layer is separated into worn sandstone and colluvium.

The shading diagrams of slope displacement under various conditions are displayed in Figure 9, Figure 10 and Figure 11. Figure 9 shows that the forces operating on the slope were concentrated on the sliding surface and a portion of the slope top because Site 2 was at the normal groundwater level. With the high groundwater level, Figure 10 shows t that the force acting on the slope moved up to the top of the slope. When seismic conditions were added, Figure 11 shows t that the shading diagrams of slope displacement looked like those with a normal groundwater level. The displacement of Site 2 in each scenario is shown in Table 8.

5.3.3. Case simulation of Site 3

Table 9 displays the soil properties for Site 3. Table 10, Table 11 and Table 12 display the material specifications of the slope protection structures at Site 3. Figure 12 displays Site 3’s analysis model. The retaining wall is shown by the blue line, the ground anchor's free section is shown by the black line, and the anchor segment is shown by the red line. From top to bottom, the soil layer is separated into a colluvial layer and sandstone and shale that are interbedded.

Slope displacement shading distribution diagrams under various conditions are displayed in Figure 13, Figure 14 and Figure 15. Figure 13 indicates that the force acting on the slope was focused below the ground anchor since Site 3 was at normal groundwater level. However, as Figure 14 illustrates, the force acting on the slope was focused at the top in the high groundwater level scenario. As seen in Figure 15, the pseudo-static analysis result revealed that the force acting on the slope extended from the slope's top to its toe when seismic conditions were included. The force acting on the slope was found to be centered beneath the ground anchors. The force acting on the slope, however, was focused near the top of the slope in the high groundwater level analysis. The force acting on the slope extended from the slope's top to its toe, according to the pseudo-static analysis result shown in Figure 15. Table 13 displays Site 3's displacement in each scenario.

6. Conclusions

Due to Taiwan's high-risk slopes on National Freeways 1 and 5, the PLAXIS 2D CE program was utilized in this study to simulate the slope stability of many existing slope protection schemes. The following conclusions were drawn from the numerical simulation results:

• Because the gradients at Sites 1 and 2 were similar, the simulation results show that, with a normal groundwater level, the force acting on the slope was mostly located at the position of the sliding surface. In the high groundwater level analysis, the force acting on the slope extended to the top.
• Due to the superior strength of the soil, Site 1 saw comparatively little displacement under both normal and high groundwater levels during the analysis.
• The displacement at the top of the slope increased as the groundwater level rose, regardless of whether the slope analysis was conducted with a normal or high groundwater level.
• The results of the simulation research demonstrate that the slope protection measures in place at a number of high-risk slopes on Taiwan's national freeways were secure under all conditions.