Introduction

A ground surface that is inclined at an angle to the horizontal is called a slope. Slopes are formed by inclining soil masses and can occur naturally or can be manmade. Slopes are commonly required in the construction of highways, railway embankments, earth dams, levees and canals.

When the ground surface is inclined, the gravity component acting parallel to the slope tends to move the soil mass downward. To ensure stability and prevent failure, it is therefore essential to analyze the slope behavior.

Various analytical and numerical techniques exist for slope stability assessment. In this example, a realistic slope is analyzed using the Finite Element Method, accounting for the elasto-plastic behavior of the soil. Specifically, a cohesive, clay-like soil is subjected to its self-weight, body load due to gravity and an additional vertical load applied near the slope face. The objective is to investigate the soil’s response and evaluate the safety of the slope in terms of potential sliding or instability [Chowdhury, 2009].

Cohesive soils are materials whose particles adhere to each other due to electrochemical attractions and water–clay interactions. For cohesive soils, a common yield model is the Drucker–Prager or Mohr–Coulomb criterion, which relates shear strength to the soil’s cohesion and internal friction angle . The yield condition can be expressed as a function of the mean stress and the deviatoric stress, capturing how the soil resists both shear and volumetric deformation.

The following application first presents the general elasto-plastic response of the loaded slope. Then, the position of the vertical load is varied to investigate different types of failure mechanisms. In the end, the slope is assessed in terms of its factor of safety, which quantifies how far the current loading condition is from the critical failure state [Yang, 2025].

For an introduction into the subject of plasticity, please consult the Plasticity monograph.

Geometry and Load Scheme

The model geometry represents an embankment with a sloping face subjected to a uniform surface load applied near the crest. Owing to the slope’s uniformity along the horizontal direction, the problem can be idealized as a plane strain case, accounting for its depth. The geometry is first defined as a 2D polygon and subsequently meshed with refinement in regions expected to develop plastic strains earlier than others.

3D representation of the soil with a sloped surface and a load applied near the slope crest.

The base of the model is assumed to be fixed, representing the contact with the underlying ground not depicted in the figure. The left boundary is positioned far enough from the slope to remain unaffected by its deformation, effectively modeling an extension of the surrounding terrain. To limit computational cost, only the portion shown is simulated, with a vertical sliding boundary condition applied to the left face, allowing vertical movement while restricting horizontal displacement. Although such a condition can influence the slope response if the domain is too narrow, this is not the case here.

Similarly, the bottom-right region is constrained to move only vertically under the load transmitted by the sloped section, while horizontal motion is restricted due to the confining effect of adjacent ground. The soil mass is substantial, and given the clay-like properties of the material, self-weight effects are also considered as an additional mechanical load acting on the slope.

A two-dimensional view of the soil. The bottom boundary is fixed, and the lateral boundaries are constrained to move only in the vertical direction. The load is applied near the crest.

In the following, the geometrical features are defined, and the domain is meshed, with particular attention to areas expected to experience high stress and thus be involved in the plastic analysis.

It is worth noting that, given the dimensions of the top and bottom bases and the overall height, the geometry could also be represented as a 2D polygon and meshed with a uniform boundary grid. However, the loading configuration indicates that the vertical load and the sloped face are the primary sources of high stress concentrations and potential plastic zones. For this reason, a nonuniform mesh is adopted.

Each polygon edge is first discretized into a sequence of points with a specified minimum spatial resolution. This can be efficiently achieved using the ToGradedMesh function, which maps each segment into a graded distribution of points. The polygon is thus represented by a denser set of boundary nodes in critical regions. The final mesh is then generated by imposing an element size distribution that follows the graded boundary spacing, ensuring finer resolution where higher stress gradients are expected.

As mentioned earlier, the top and bottom sides require higher mesh resolution near the right boundary, while the sloped face demands finer discretization toward its ends to better capture stress gradients. In contrast, the left boundary does not require special refinement; its resolution is automatically determined during the meshing process once the full polygonal domain is generated.

It should be noted that the primary purpose of the graded mesh is twofold: to ensure high resolution in critical regions while keeping the total number of elements limited, thereby reducing computational cost. This aspect is particularly important for elasto-plastic analyses, which require multiple nonlinear iterations and significantly longer solving times compared to purely elastic simulations.

Soil Behavior

The material considered is a clay-like soil, which is prone to the formation of high shear zones that can lead to plasticization and potential sliding. However, its cohesive properties contribute to resisting deformation and delaying failure by increasing shear strength. The Drucker–Prager and Mohr–Coulomb yield criteria have both shown good capability in representing this type of soil behavior. In the following analysis, the Drucker–Prager model is adopted. It is worth noting, and will be discussed later, that under specific assumptions, these two criteria can be considered equivalent and compared by properly scaling their parameters.

To solve the elasto-plastic problem, an initial purely elastic analysis is first performed using the defined material parameters and boundary conditions. The resulting elastic stress and strain fields serve as the initial approximation for the subsequent elasto-plastic solution, which is obtained by iteratively solving the same governing equations while incorporating the plasticity model.

The external load is applied on the top surface, near the sloped face. Since the governing equations will be solved multiple times, an auxiliary function is defined to handle not only the material parameters and problem variables but also the geometric domain and the loaded region. This modular approach allows the automatic generation of the corresponding PDEs while accommodating changes in domain properties or load configuration.

The deformation results show a settlement of the top surface caused by both the self-weight of the soil and the applied external load. Before proceeding with the elasto-plastic analysis, the yield criterion can be evaluated based on the current stress state.

The evaluation of the yield criterion requires computing the first stress invariant and the square root of the second invariant of the deviatoric stress, . These quantities provide insight into the proximity of the material to yielding, as defined by the Drucker–Prager criterion. Further details can be found in the Drucker-Prager section in the Plasticity monograph.

For stiff clay-like soils, the Drucker–Prager (DP) parameters indicate a cohesion parameter greater than zero and a noticeable level of pressure sensitivity [Rafiei, 2020].

As shown in the stress plot, the region near the toe of the slope experiences elevated stress levels that exceed the yield limit. Therefore, an elasto-plastic analysis is required to accurately predict the final deformation and stress redistribution within the soil mass.

When simulations run for a long time, it may make sense to store the obtained results to hard disk such that they are readily available later without rerunning the simulation.

The elasto-plastic solution shows a pronounced settlement in the loaded area. This behavior results from plasticization, which redistributes part of the load from the vertical direction into horizontal shear, causing a slight displacement of the sloped face toward the right. However, this horizontal shift remains limited and does not yet indicate a critical risk of slope failure. The potential for failure will be analyzed later; first, by examining the principal stress directions, it is possible to identify the likely orientation of the critical plane that could initiate and propagate a failure mechanism within the slope.

Before proceeding with the discussion, it is useful to clarify an important notation convention. In solid mechanics, tensile stresses are typically considered positive and compressive stresses negative, as shown in the Solid Mechanics monograph. However, in soil mechanics, the convention is reversed: compressive stresses are often taken as positive since they are the most relevant.

Here the standard convention is used, where the positive value represents a tensile stress. That is why in the following compression plot the scale is negative.

Furthermore, principal stresses are ordered from the largest, , to the smallest, . Consequently, in soil mechanics applications, greater attention is often given to the third principal stress, since a higher compressive magnitude, i.e. a larger positive value under the soil convention, represents a critical stress condition for the material.

The Drucker–Prager and Mohr–Coulomb criteria describe soil yielding as governed by a combination of shear stress and mean normal stress. These models account for the influence of confinement and friction between soil particles. Failure occurs when the shear stress acting on a plane reaches a critical value determined by the material’s cohesion and internal friction angle.

The potential failure planes correspond to orientations where the shear stress to normal stress ratio reaches its maximum, typically at angles related to the material’s friction angle. Along these planes, plastic yielding initiates, marking the onset of irreversible deformation and the development of shear bands in the soil mass.

The Drucker–Prager and Mohr–Coulomb parameters can be converted one to the other, as equivalent forms. To facilitate the analysis of the failure planes, it is more convenient to refer to the angle of internal friction of the Mohr–Coulomb model, since it directly governs their orientation. For additional details, refer to the Plasticity monograph.

The failure lines are modulated by the internal friction angle, which shifts their orientation relative to the maximum shear direction. While the maximum shear direction rise is 45° with respect to the line connecting the first and third principal stress directions, in this case the angle became .

The failure criterion predicts two symmetric directions oriented at equal angles about the principal stress directions. The positive direction and the negative direction are computed using a positive and negative sign for the angle , respectively. The streamlines of these directions define the lines of potential failure.

Failure happens along the curved red lines in the high shear stress areas.

The failure lines are shown in blue and red and their opacity is scaled with respect to maximum shear stress . The red lines represent the classical sliding pattern typically observed in slope failures. The blue lines show the mechanically equivalent conjugate failure line that completes the pair of possible failure orientations predicted by the Drucker–Prager and Mohr–Coulomb models. Although both orientations satisfy the failure criterion mathematically, only one of them, the red one, is physically activated. This line cannot be determined a priori, but is selected by the physics of the problem, depending on the slope geometry, loading direction and boundary conditions.

The plot also highlights regions of high shear stress beneath the loaded region with whiter areas. Together with the curved red failure lines, these illustrate the characteristic rotational movement that develops along these paths during a slope sliding failure.

Slope Failure Zones

In slope and foundation stability analysis, different failure mechanisms can develop depending on geometry, loading and soil strength. Three slope failure modes will be investigated next. As will be seen, the previously computed solution is typical for a face failure. The naming will become apparent later. Furthermore the toe and Prandtl failure modes will be introduced and compared. Together, these mechanisms describe the range of possible failure modes in soils subjected to shear and compressive stresses.

In the following, the load position is varied to highlight different failure mechanisms. Before proceeding, all the components required to solve the elasto-plastic behavior of the soil are encapsulated in an auxiliary function.

To illustrate the different failure mechanisms, the load position is varied. Before proceeding, all the components required to solve the elasto-plastic behavior of the soil are encapsulated in an auxiliary function.

In order to investigate a different failure mechanism, the load position is changed while keeping the geometry and material properties unchanged.

The plot highlights two distinct failure mechanisms. The first shows a more localized displacement concentrated near the upper part of the slope, gradually decreasing toward the base. In contrast, the second mechanism exhibits a broader zone of deformation, with a lower peak displacement but a wider affected area, particularly extending from the slope base.

To better highlight the difference between the two failures, the horizontal displacement is extracted along the sloped edge.

The plot shows that the horizontal displacement of the slope, directed toward the sliding side, is greater in the face failure, where the load is applied close to the exposed slope face. In this mechanism, the face itself tends to collapse outward, hence the name face failure. By contrast, toe failure is characterized by movement concentrated in the lower part of the slope, where the failure surface emerges at the bottom section, known as the toe of the slope. In the plot, the effect of face failure is most pronounced near the crest, at point A, while at the base, point B, the displacement values are similar because they are constrained by the external boundary.

Both face and toe failures are types of sliding mechanisms associated with slope instability. Another form of soil failure is the Prandtl failure, which is instead related to the classical settlement behavior of the soil under loading. To observe this effect, the load must be positioned sufficiently far from the slope. This can be achieved by extending the right side of the domain while maintaining the same slope angle.

Face failure occurs when the slip surface originates near the slope face and extends outward, leading to a shallow, progressive sliding of the surface material. Toe failure develops when the failure surface passes through or just beneath the slope toe, typically involving deeper soil layers and resulting in larger displacements. Prandtl failure, on the other hand, represents a classical plastic mechanism observed beneath loaded areas such as footings, characterized by a symmetrical fan-shaped pattern of shear zones connecting to the surface, not shown here.

The failure mechanism can also be interpreted from the failure lines. For this geometry, the negative failure direction , obtained using , is sufficient to represent the expected failure orientation. Furthermore, when scaled by the maximum shear stress, these lines directly illustrate the predicted failure pattern.

The failure lines are represented in black, and their opacity is scaled with respect to the maximum shear stress value . In all cases, the failure lines originate from the loaded region. Face and toe failure mechanisms exhibit similar patterns, but in the toe failure the mechanism is shifted away from the slope face, producing an indirect effect by pushing into the lower part of the slope. In contrast, the Prandtl failure is largely confined to the loaded area. Although some failure lines remain visible near the slope face due to the self-weight of the soil, the maximum shear intensity there is very low, indicating that collapse in that region is unlikely.

As previously introduced, there is a clear distinction between the different soil failure mechanisms, most of which arise from the geometrical configuration of the slope. Due to the combined effects of external loading, self-weight and domain geometry, the slope may fail through face failure, where displacement is concentrated near the slope crest, or through toe failure, where a large portion of the soil mass begins to move from the slope base, dragging the entire slope with it.

The Prandtl-type failure, on the other hand, is characterized by a localized plasticization of the soil, leading to increased displacement compared to a purely elastic model. It primarily represents a settlement-type mechanism rather than a slope instability.

It is also worth noting that the general magnitude of displacement remains comparable across all configurations, which is mainly associated with the vertical component. While the externally loaded zone undergoes movement, the effect of the soil’s self-weight remains constant across the three geometries, producing a compressive state that can contribute to slope collapse under critical conditions.

Shear Strength Reduction

To compute a factor of safety for slope stability, the shear strength reduction method is commonly employed [Rafiei, 2020]. In this approach, the soil’s plasticity parameters are gradually reduced by a strength reduction factor, , until the model reaches failure. The value of at which numerical instability occurs represents the safety factor of the slope.

To provide a clearer physical interpretation of the parameters involved, it is preferable to discuss them within the Mohr–Coulomb framework, since the corresponding conversion from the Drucker–Prager formulation has already been introduced. The two key parameters are the cohesion parameter , and the angle of internal friction , which together define the shear strength of the soil. They determine the maximum admissible shear stress that the soil can sustain for a given normal stress acting on a potential failure plane, expressed as .

The shear strength reduction method progressively decreases the shear strength parameters by a reduction factor in order to determine how close the material is to its failure limit.

The values for and have been computed above.

In the following section, the shear strength reduction method is implemented. To simplify its application under different loading conditions and mesh configurations, all necessary procedures are wrapped in an auxiliary function. This function performs the following steps:

• Reduction of the Mohr–Coulomb parameters according to the specified reduction factor , and their conversion to the equivalent Drucker–Prager parameters.

• Initialization of the elastic solution as the starting state.

• Execution of the elasto-plastic iterative loop to update stresses and plastic strains.

• A convergence check to determine whether the elasto-plastic solution has reached failure.

Regarding the last step, in this case, the elasto-plastic solution is intentionally advanced toward numerical failure as the shear strength reduction method requires. However, due to the formulation of the elasto-plastic problem, convergence of the plastic flow may still be achieved even when the material is plastically yielding, resulting in a continuously increasing plastic strain. Although such a case does not represent numerical divergence, it leads to progressively larger displacements and an eventual loss of physical stability and meaning.

To account for this behavior, a tolerance threshold is introduced, based on the displacement error computed at the last iteration. The iterative displacement error, which can be extracted from the updated solution parameters, determines if the elasto-plastic computation is valid or not. While increasing the factor , the last valid iteration determines the factor of safety. To speed up the computation, a lower number of iterations is used.

The loop terminates automatically once the solution no longer achieves convergence. The last step therefore corresponds to the maximum value of the reduction factor , which represents the factor of safety.

A factor of safety below 1 is not permissible, while in practice, higher factors of safety are required.

The distribution of plastic strain shows that plasticization develops primarily along the slope, aligning with the characteristic failure surfaces associated with slope instability.

The plot illustrates the deformation pattern at the final iteration of the shear strength reduction analysis, in red, showing a clear movement of the slope toward failure. The displacement at the crest is significantly higher than that observed with , obtained with the original material parameters, in blue.

Slope Angle

As a final consideration, the inclination of the slope plays a fundamental role in stability analysis. The slope angle is one of the most significant geometric factors governing failure, and instability may occur solely due to the soil’s self-weight when the angle exceeds the material’s shear resistance.

In the following, the slope geometry is modified from its original inclination to four different slope angles, and for each configuration, the factor of safety is evaluated.

While all results indicate that the factor of safety is above 1, the last results are just barely above it.

The factor of safety (FoS) decreases sharply as the slope angle increases. In particular, the last mesh, characterized by the steepest face, shows a very low value close to one, indicating almost no margin of stability. This is also confirmed by the deformed mesh, which clearly illustrates the progressive collapse of the slope. The failure is caused by the combined effect of several mechanisms: the self-weight of the soil compresses the material and induces a slight downslope movement; the external load further compresses the crest area, increasing the stress concentration in the underlying region; and finally, the accumulation of plastic strain needed to accommodate the extreme load leads to large displacements that drive the slope toward sliding and eventual failure.

3D Visualization

As a final step, since the example originates from a 3D slope geometry but is analyzed under plane strain conditions, it is possible to map the 2D deformation field onto the 3D model by assuming that the deformation remains constant along the embankment width.

The displacement field can be mapped to three dimensions by assigning the planar coordinate to the 3D direction, and the planar coordinate to the 3D vertical direction, while the 3D direction remains constant with zero displacement.

To highlight the slope failure, several points along the sloped face are extracted to plot the displacement vectors. To achieve this, it is necessary to determine the equation of the plane representing the slope surface and then extract the points lying on that plane.

To achieve the desired effect of simulating the slope failure, the displacement field can be progressively scaled over several steps, generating a sequence of images that can be combined into an animation.

The animation illustrates how the slope progressively reaches failure through sliding. Under the effect of the external load, the soil mass is compressed against the lower section, resulting in a horizontal displacement of the slope that evolves into a sliding failure. It should be noted that the final frames of the animation appear overly deformed due to the displacement scale factor applied to enhance visualization of where and how failure develops in the case with the steepest slope.

References

1.  Chowdhury, R. (2009) Geotechnical Slope Analysis. ISBN: 978-1-135-19210-5. https://doi.org/10.1201/9780203864203.

2.  Yang, G. (2025) "Ultimate Bearing Capacity of Embedded Strip Footings on Rock Slopes under Inclined Loading Conditions Using Adaptive Finite Element Limit Analysis." Scientific Reports, 15-1. https://doi.org/10.1038%2Fs41598-025-91234-2.

3.  Rafiei, R. (2020) "Factor of Safety of Strain-Softening Slopes." Journal of Rock Mechanics and Geotechnical Engineering, 12-3. https://doi.org/10.1016/j.jrmge.2019.11.004.