Deformation and failure analysis of heterogeneous slope using nonlinear spatial probabilistic finite element method

1.   Introduction

Slope failures refer to the sudden or gradual movement of soil, rock, or debris down a slope. These failures can range from small landslides to catastrophic events, posing substantial dangers to infrastructure, natural ecosystem, and people's lives ^[1,2]. Slope failures can result from natural causes like excess precipitation, seismic movement, environmental weaknesses, and erosion ^[3]. Slope instability can also be caused by human-induced processes, including mining, building, deforestation, and inappropriate land use. By comprehending the origins, methods, and consequences of slope failures, researchers can develop effective strategies to mitigate risk and increase strength in failure-prone regions ^[4,5,6].

The mechanisms of slope failures vary depending on factors like slope geometry, geological conditions, and environmental forces ^[7]. Major failure mechanisms include sliding along weak planes or surfaces, stepper slopes subjected to gravitational forces and external loads, etc. The safety factors or failure characteristics for a slope are estimated in two major ways: Traditional limit equilibrium methods (LEM) and numerical analysis approaches. Many researchers offer comprehensive examinations of limit equilibrium methods of slope stability exploration ^[8]. Many LEM-based methods, such as the ordinary method of slices, force equilibrium method, Bishop's modified method, Jambu's method, Morgenstern & Price's method, and Spencer's method, were employed to determine the failure surface and factor of safety values ^[9,10]. The soil or rock mass must be separated into slices to utilize LEM methods. It is assumed that the forces acting on each slice of the slope are oriented in certain ways. The ability to differentiate among different limit equilibrium approaches depends critically on the assumption ^[11]. The Methods of limit equilibrium necessitate a continuous surface that traverses the soil bulk. The surface is crucial for determining the minimum factor of safety (FOS) against sliding or shear failure. Limit equilibrium methods require the artificial construction of equations of equilibrium, including side forces and their directions, before computation of slope stability. Moreover, LEM methods do not provide information about comprehensive deformations. The variability in slope geometry and material nonlinearities introduced more complexity in assumptions and further analysis using limit equilibrium methods ^[12,13].

In other parts, numerical methods are a better approach to defining, understanding, and analyzing slope stability problems in general. Finite element methods (FEM), finite difference methods (FDM), and discrete element methods (DEM) are prominent numerical methods used to analyze slopes ^[14,15]. FEM enables the analysis of more complex geometry, material heterogeneity, and boundary conditions compared to limit equilibrium methods. Slope geometry, material properties, and loading conditions are discretized into finite elements in finite element analysis. Many numerical approaches-based commercial software are available to analyze slopes with predefined procedures for designing and analyzing the slope in a constraint software environment. The numerical methods involve large numbers of equations that need to be solved simultaneously, which arises from complexity, high computation times, and large numbers of system variables. MATLAB is one of the software programs that uses matrix datatypes to efficiently handle large numbers of system variables and provides a suitable platform to deal with many simultaneous equations and looping operations ^[16,17].

In plane strain analysis, it is assumed that the slope mass under consideration is very long in one direction (along the plane), and there is no change in dimension along that direction. This enables the simplification of the three-dimensional problem into a two-dimensional one. Plane strain analysis greatly reduces the system complexity with certain assumptions ^[18,19]. The assumption includes that the deformation in the direction perpendicular to the plane of analysis is negligible, which means the strain change in the out-of-plane direction is zero. Additionally, the material properties are assumed to be uniform in the out-of-plane direction. Many commercial tools have employed the plane strain method to provide a platform for two-dimensional deterministic slope stability analysis ^{[20,21,22,23]}.

Typically, many researchers employ analytical and numerical methods to evaluate spatial and temporal variability on slopes to determine safety factors. In fact, slope failure can occasionally result from a deterministic strategy that overlooks spatial variability, even in cases where the "deterministic" factor of safety (FoS) is greater than one. Spatial inconsistency is essential as the behavior of the landscape is influenced by its geo-mechanical characteristics based on deformation, and failure typically starts and spreads from the most susceptible region. Probabilistic methods provide more systematic information, expressed in terms of reliability based on deformation and failure probability, allowing for a more robust evaluation ^{[24,25,26,27]}. The implementation of probabilistic methods to analyze slope stability problems is not new and has been in use since the 1970s ^{[28,29,30,31]}. These methods involved combining limit equilibrium techniques with various probabilistic approaches ^[32]. The probabilistic analysis offers several advantages. Most of the methods have incorporated random finite element methods (RFEM) and the Monte Carlo Simulation (MCS) approach ^{[32,33,34,35,36]}. The RFEM and MCS techniques compute the likelihood of slope failure by utilizing the uncertainties in material properties ^[32]. The MCS method considers the randomness of variables while also incorporating temporal variation. It allows for a straightforward determination of standard errors in the results. Furthermore, the size of the uncertainties does not affect the approach's convergence speed. Furthermore, the simulation procedure does not influence the sophistication of the execution function. One drawback of the MCS and probabilistic methods is that when the likelihood of failure decreases, the number of simulations and time required for calculation increase significantly. Numerical computations of complex systems also utilize probabilistic methods, particularly MCS, to modify the topography (plastic and/or elastic) characteristics of rock joints ^[37]. Similarly, the finite difference method (FDM), discrete element method (DEM), and finite volume method (FVM) have been embraced in recent years ^{[38,39,40,41]}. In the RS-DEM method (Random Set Discrete Element Method), the focus is on considering just the gaps of the input parameters rather than the actual values of the parameters themselves ^{[42,43,44,45,46,47,48,49]}. A variety of factors, including the composition of different materials, their granulometry, the construction method, the interaction with runoff, and consolidation processes, cause heterogeneity, also known as spatial and temporal variability. Instead of assuming a simple vertical gradient, as is common for uniform soil layers, the presence of such a wide range of differences in attributes supports the idea of examining their spatial distribution ^{[50,51,52,53,54]}.

Major factors contributing to spatial uncertainties in slope materials include measurement methods, sensor placements, topography, vegetation, weathering, evapotranspiration, etc. For instance, utilizing various testing and empirical methods to measure slope material parameters always raises questions about the precision of the measured values. In general, instruments measure quantities with a 5–10% tolerance, and these are widely acceptable in measurement and instrumentation. Repeated testing of the rock or soil samples reveals variations in the measured values of material properties. It is also not possible to measure the material parameters for each space point of the slope. Therefore, it is crucial to consider the random variation in material parameters when creating a slope model. The fusion of probabilistic methods and finite element modeling provides a systematic framework to incorporate and analyze the uncertainty in slope parameters. This probabilistic finite element method maps the spatial uncertainties in material properties using a probability distribution function (PDF) and coefficient of variation (CoV).

In the absence of comprehensive data, it is reasonable to use a uniform probability distribution for material attributes in preliminary slope stability analyses. This offers a useful method for modeling variability and uncertainty, making initial evaluations easier to understand and analysis. This method prevents bias, and guarantees a wide range of potential values, and it is simple to incorporate into computer models. Uniform distributions provide a useful foundation for understanding slope behavior and guiding preliminary design decisions, even though they are an oversimplification and can improved and compared with more precise data and other probability distributions in the future study. In general, horizontal correlation scales are often greater in sedimentary deposits, but various factors, such as complex depositional environments, post-depositional disturbances, tectonic activity, natural discontinuities, and varying sediment sources, can result in comparable or greater variability perpendicular to the layers. These factors suggest that the assumption of a greater correlation scale along the layer is not always valid and requires site-specific consideration.

In this study, the benchmark two-dimensional slope problem is taken from the literature authored by Arai and Tagyo in 1985 ^[55]. The problem resides in multi-layered cohesive soil layers that have certain slope inclinations and heights. This benchmark problem has been also used by many researchers for instability analysis ^{[56,57,58,59,60,61]}. The benchmark slope geometry and material parameters are presented in Table 1, where the additional probability distributions in material properties are incorporated. The material properties are assumed spatially random and distributed around a deterministic value to incorporate more realistic mapping and measuring uncertainties in this study. So, the benchmark slope model is additionally incorporated by uncertainty of uniform pdf and 10% CoV in all six material properties. The research elucidates the deformation and failure characteristics of heterogeneous slopes, highlighting the influence of material spatially uncertainties on the distribution of displacements, stresses, and strains, and failure region. The nonlinear probabilistic finite element method (PFEM) allows a detailed investigation of detailed deformations, stresses and strains patterns, and the investigation of localized failure mechanisms in the slope. Focusing on plane strain conditions allows for a simplified yet realistic representation of slope behavior. The study captures important aspects of slope response while making computations a lot simpler by looking at two-dimensional deformations in a plane perpendicular to the direction of loading.

The research work is to develop and analyze a MATLAB code-based mathematically simulated model of the testbench heterogeneous slope. The work originated with the modeling of a two-dimensional benchmark heterogenous slope consisting of three different layers. The slope is also incorporating the spatially distributed material properties (cohesion, friction angle, elastic modulus, dilation angle, and unit weight) with different degrees of variations in terms of %CoV. In this way, nonlinearities and randomness are introduced additionally to the slope model in terms of geometry, PDF, and CoV using the probabilistic finite element method approach. Further, Mohr-Coulomb's soil model and failure conditions are utilized for elastoplastic analysis of the slope and produce detailed deformations and slope failure characteristics under self-weight and step loading scenarios. The failure characteristics of the slope include the determination of the pre-failure zones, shear stress patterns, different strain patterns, and the probability of a nodal failure parameter. Also, new quantitative indicators, 'nodal failure' and 'probability of nodal failure' are presented to add more depth to failure investigations. Moreover, the deformations and failure characteristics of the testbench slope are examined at various levels of uncertainty by altering %CoV.

2.   Methodology

The systematic process flow chart of the research work is shown in Figure 1. The flowchart illustrates the four key stages of the proposed slope stability study: Initialization, pre-processing, elastoplastic analysis, and post-processing.

The initialization stage involves defining slope geometry, slope material parameters and their non-correlated random distribution, finite element parameters, a generalized slope constitutive model, and elastoplastic analysis parameters. In the preprocessing stage, the finite element method is initially used to identify elements and nodes for the given slope geometry. The material properties are demarcated for each node, including the spatial randomness of the space. The boundary conditions are managed to include the restricted motion in the X and Y directions. The global parameter matrices are calculated using the spatially distributed material parameters for the slope. Further, the initial stress profiles are also determined using these global parameters and the self-weight load. In the next stage, the elastoplastic analysis is employed for the slope model using the generalized Mohr-Coulomb's slope model. In this process, the slope model is subjected to uniform step-loading at the top edge. In the final stage, the impact of step-loading is observed in the form of detailed profiles of deformation, stresses, and strains. Moreover, information on failure nodes is used to identify pre-failure zones and nodal failures, and the probability of nodal failures is assessed specifically. Furthermore, the constructed model and processing techniques are simulated for various values of %CoV to examine the implications of uncertainties on slope physical characteristics and slope behaviors using a limit analysis scheme.

2.1.   Finite element & probabilistic formulation

The finite element method discretizes the two-dimensional slope geometry. In finite element analysis (FEA), nodes and meshes are fundamental components used to discretize and analyze complex geometries and solve the governing equations ^[62]. Meshing enables the analysis of complex structures by dissecting them into simpler sections, known as elements. These elements collectively cover the entire slope's domain. Each element in the mesh is defined by a set of nodes known as its connectivity. Figure 2(a) illustrates the definition of node connectivity, which is defined as an anti-clockwise sequence for each element of the slope.

This connectivity also determines how the shape functions interpolate the field variables within the element. Nodes are also defined as discrete points within the geometry being analyzed. They serve as locations where various quantities, such as displacement, stress, and strain, are calculated or applied. Each node typically has associated degrees of freedom (DoF), representing the directions in which the node can move or experience deformation or restraint.

The choice of mesh density, element type, and nodal distribution significantly impacts the accuracy and efficiency of the analysis ^[62,63]. Finite element analysis employs a variety of shape-based elements, including triangles, quadrilaterals, tetrahedral, and hexahedral, based on the geometry and the desired level of accuracy ^[64]. The present analysis discretizes the slope using quadrilaterals, bilinear, and four-node elements, and represents it as a 'Q4' element in two-dimensional analysis.

This quadrilateral-shaped Q4 element is intended to section the slope with minimal loss of coverage. It provides the optimum coverage with less complexity and loss of generality compared to other shaped elements. More elements or higher-order elements can be used, but the number of governing equations, complexity, and simulation time are increased effectively, but no significant changes in output profiles are observed.

The behavior of the field variables within each element is interpolated using shape functions defined at the element nodes. The shape function is defined in terms of the natural coordinate system, as shown in Figure 2(b). The natural coordinate system enables the specification of a point within an element through a spectrum of dimensionless values, ensuring its magnitude never surpasses unity. The two natural coordinates within the element are represented by ξ and η in the X direction and Y direction, respectively, where -1 ≤ ξ ≤ 1 and -1 ≤ η ≤ 1. Equation (1) provides the shape function for a Q4 element.

These shape functions can be derived based on the bilinear interpolation of the coordinates of four nodes of the element. Each shape function f (ξ, η) represents the contribution of the corresponding node to the overall behavior within the element. They are linear combinations of the natural coordinates (ξ, η) and satisfy the Kronecker delta property, which means evaluating them to 1 at the node they correspond to and 0 at all other nodes for an element ^[65].

These shape functions are then used to interpolate the field variables within the element based on the nodal values. They are also used to calculate the derivatives of field variables for numerical integration and other calculations within the element. The bilinear shape functions for the Q4 element provide a simple and efficient way to interpolate within the element and can capture complex deformations efficiently.

The slope on MATLAB is developed using a versatile constitutive Mohr-Coulomb's soil model, illustrated by six material properties ^[66,67]. The six material properties or parameters of slope are Young's modulus (E), Poisson's ratio ($\upsilon$), unit weight ($\gamma$), cohesion (c), friction angle ($\mathrm{\varnothing }$), and dilation angle ($\psi$). The finite element method incorporates slope heterogeneity by defining different value sets of these parameters for different segments of slope and fused probabilistic behavior by randomizing these parameters for stochastic analysis. After that, the global parameter matrices of these material parameters are established to incorporate features of nonlinear, probabilistic finite element-based heterogeneous slope.

One of the global parameters is the global elastic constitutive matrix (C), which is used by system-governing equations to establish relationships among field variables. Equation (2) presents the general formula for a node's elastic constitutive matrix (Cn). The global elastic constitutive matrix (C) comprises the collections of elastic constitutive matrices evaluated for all nodes of the benchmark slope.

The primary field variable of the finite element analysis is the displacement vector (U), which is defined in two spatial coordinates (u, v) of the cartesian coordinate system. The mathematical formulation of displacements is given in Eqs (3.1) and (3.2) in terms of shape functions. Here, u and v are the horizontal and vertical displacements, respectively, for a point within the Q4 element. Also, u_i, v_i, and N_i are horizontal, vertical, and shape functions of i^th node of the element. The nodal displacement is represented by the U vector, as given in Eqs (3.3). In the finite element analysis for the benchmark slope, the field variable horizontal and vertical displacements are determined at each node of the slope mesh.

The partial derivative of shape functions with regards to the cartesian coordinates is used to calculate the strain-displacement matrix (B) using the mathematical relationship mentioned in Eq (4).

The strain displacement matrix (B) and global linear constitutive matrix (C) are further used to calculate strains and stress tensors, as represented in Eqs (5.1) and (5.2).

These stress and strain states are defined for each gaussian point of all elements in the slope domain. The stress and strain states have a total of four components in different planes, as shown in Eqs (5.1) and (5.2). A strain state can be computed in four directions with a cartesian coordinate arrangement (xx, yy, xy, zz). The zz-components of the strain state are in a long spatial direction and assumed to be zero considering the plane strain condition.

Equation (6) uses this strain-displacement matrix and the linear constitutive matrix to determine the global stiffness matrix.

Here, the weak form of the governing partial differential equations leads to numerical integration using the finite element method. The present analysis employs the gaussian point integration method with two integration orders. In the present work, FEM-based computational methods can efficiently and accurately compute these integrals, leading to accurate solutions to governing equations. Equation (7) illustrates the further use of the global stiffness matrix to determine nodal displacements with the application of nodal force (F_n).

The fusion of the boundary conditions, which include the nodal degree of freedom in calculation and analysis, also simplified the finite element analysis. These conditions are incorporated in the corresponding displacement vector and their degree of freedom in such a way that deformations are only calculated at unrestrained nodes in the slope.

2.2.   Yield function & elasto-plasticity

A stress state is enumerated with four components with a cartesian coordinate arrangement (σxx, σyy, τxy, σzz) for each gaussian point available on the slope surface. The global stress matrix comprises four components of the stress value of all four gaussian points for each of the elements available on the slope surface. The stress state can also be represented in the form of a principal stress coordinate system using (σ₁, σ₂, σ₃) or in the form of stress invariants (p, q, $\theta$). These stress invariants (p, q, $\theta$), along with cohesion and frictional angle, form the failure condition or yield function (F), as mentioned in Eqs (8.1) and (8.2).

This yield function or failure function (F) is used to determine the elastic or plastic state for Mohr-Coulomb's model of elastoplastic slope and distinguish between elastic and plastic states. The yield function is determined and assessed for each node of the slope separately.

An elastoplastic material acts as either an elastic solid or a plastic solid, depending on the present stress state. If the value of the yield function is less than zero, then the point exhibits elastic behavior; otherwise, it is plastic. The yield or failure criteria create a surface in the three-dimensional principal stress space and describe the change from elasticity to plasticity as presented in Figure 3. This Mohr-Coulomb failure region has a hexagonal cross-section. Moreover, the yield surface in a two-dimensional stress space is represented in Figure 3. Stress states are classified as plastic when they are on the yield surface and as elastic when they are inside the surface. Stress states that are outside of the failure surface are unacceptable and need to be reallocated iteratively because when the material distorts plastically, stresses must stay on the failure surface.

In general, two methods have been used to solve iterative finite element-based governing equations. These methods are the tangent stiffness method and the constant stiffness method. The present work utilized constant stiffness techniques to adjust the loads on the system iteratively in order to attain convergence through repeated elastic solutions.

The displacement variation (${\mathit{\Delta}} u$) for each load increment (${\mathit{\Delta}} f$) must be solved, and the strain-displacement matrix is then utilized to obtain the total strain increments as presented in Eqs (9.1) and (9.2).

Changes in strain states respond to the stresses during elastoplastic analysis. The failure or yield function (F) is determined by evaluating the current stress state. If the failure function is equal to zero, it means stress states at the failure surface, and less than zero represents the elastic straining that occurs. If this yield function is greater than zero, material starts deforming plastically, and present stress state should not cross the failure surface, so the strain states are distributed between plastic strain (${\mathit{\Delta}} {\varepsilon }^{p}$) and elastic strain [${\mathit{\Delta}} {\varepsilon }^{e}$] as shown in Eq (10). The yield function is determined for each gaussian point of all elements separately, which provides precise point information about the slope in the plastic region.

The change in elastic strain is further used to determine the stress component ($\Delta {\rm{ \mathsf{ σ}}}$) using the elastic constitutive matrix (C), as depicted in Eq (11). This final stress state is achieved by accumulating this change in stress state ($\Delta {\rm{ \mathsf{ σ}}}$) for all iterations.

The generated stresses ($\Delta {\rm{ \mathsf{ σ}}}$) are added to preexisting stresses, and the resultant stress is again utilized to determine the failure condition as given in Eq (12).

If yield occurs (F > 0), then slope points undergo plastic deformation. Plastic straining happens when the stress condition cannot be maintained for an extended period. The yield function's value, which expresses how much the current stress state surpasses the failure condition, determines the size of the plastic strain rate, or visco-plastic strain rate. The plastic strain rate is expressed in terms of pseudo-time step (${\mathit{\Delta}} t$), yield function, and plastic potential derivative ($\partial G/\partial \sigma$) using Eq (13), considering associated flow conditions.

The visco-plastic strains are further utilized to calculate the residual nodal force for each gaussian point. The strain state represents the accumulated strain changes from all iterations and is mathematically expressed in Eq (14).

Now, the plastic strains are further utilized to determine the residual forces, which represent the internal forces remaining in a structure or material after the loads have been removed, which means incorporating the plasticity. The residual forces arise due to the redistribution of stresses during plastic deformation and can significantly affect the behavior and stability of the structure. The collective residue forces of an iteration are again added to the previous step load, and the process repeats until convergence occurs. The convergence conditions integrate the incremental displacements, small values of yield function, and strain values, which can be modulated by tolerance limit conditions.

In elastoplastic analysis, a slope undergoes a number of external loading sequences. We utilized the yield function and Mohr-Coulomb's failure criteria to pinpoint the precise nodes where the slope yields. These nodes are referred to as failure points. This work introduces two novel quantitative failure parameters: 'nodal failure' and 'probability of nodal failure'. The nodal failure is formulated as the ratio of the failure node to the total number of nodes for each external loading sequence. The information about nodal failure can be used to investigate the different loading types and different magnitudes of loading increments on a slope. Also, the probability of nodal failure is defined as the average of nodal failures. The probability of nodal failure is characterized by certain loading types, load increments, and loading sequences. These parameters can be used to investigate slope failures and comparatively analyze the different types and levels of uncertainties introduced in probabilistic finite element analysis. The present spatial probabilistic finite element analysis investigates these parameters for the heterogenous nonlinear benchmark slope.

3.   Simulation and result analysis

3.1.   Probabilistic finite element modelling of the testbench slope

The benchmark two-dimensional slope geometry is 32 m high and 66 m wide. MATLAB code models the benchmark slope and presents it in Figure 4. Finite element analysis primarily divides this slope into finite numbers of small pieces, known as 'elements'. This is like substituting a structure with a finite number of degrees of freedom for one with a limitless number of degrees of freedom. This process is termed meshing, and the sectionized slope is known as the meshed slope. Figure 5 displays the results of simulating the slope for meshing. The slope is converted into finite numbers of Q4 elements in the X and Y directions. The element used is of quadrilateral shape and is referred to as a 'Q4' element. It is speculated that each Q4 element has four corner points, or that elements are connected at specific places known as nodes. Two spatial coordinates, X and Y, can be used to characterize the geometry, material characteristics, and field variables of the system in a plain strain analysis.

Each element and node of the meshed slope is numbered to use their geometry, material characteristics, and field variables for calculation and analysis. The first element and first node are the left-most and bottom-most element and node, respectively. Elements are incrementally numbered in a left-to-right and bottom-to-top approach. In the same way, nodes are also incrementally numbered from left to right and from bottom to top. The total number of elements in the X and Y directions is 33 and 16, respectively. The total number of nodes in the X and Y directions is 34 and 17, respectively. Figure 5 represents the finite element mesh slope using Q4-shaped elements.

The slope model can map the slope's natural support at the bottom and side boundaries using appropriate boundary conditions. The boundary condition is a basic way to restrict the movement of boundary nodes. Natural boundary conditions restrict the movement of side-boundary nodes only in the Y direction. The bottom node of the slope is restrained in both the X and Y directions. This can be achieved by defining in X and Y the directional degree of freedom matrix for the respective nodes. Figure 5 of the meshed slope depicts the restricted nodes in both directions as circles and the X-direction restrained nodes as squares. These boundary conditions are used not only to map natural conditions, but also to reduce the model's complexity. Now, deformations are calculated for unrestrained nodes of the benchmark slope model only.

The developed slope model is defined using six material properties. These are Young's modulus (E), Poisson's ratio ($\upsilon$), unit weight ($\gamma$), cohesion (c), friction angle ($\mathrm{\varnothing }$), and dilation angle ($\psi$). The benchmark slope has both material nonlinearity and geometric nonlinearity. The benchmark slope has three distinct material layers, comprising different material parameters as listed in Table 1.

We encompass the spatial probabilistic uniform distribution for material parameters, which is listed in Table 1 and is presented over the slope space in Figure 6. The random distribution for each material property has a center value, which is the deterministic value of the material property. The distribution is characterized by the percentage coefficient of variation (%CoV), which is a measure of dispersion around the deterministic value of material properties. The uniform probability distribution function with a 10% CoV in spatial values of material properties is labeled as the first randomized material properties of the benchmark slope. The frequency distributions of all material parameters for the first randomized set are represented in Figure 7. In this Figure 7, three material parameters (cohesion, friction angle, and dilation angle) are randomized around three different respective deterministic values of three different layers, as mentioned in Table 1. These histograms also confirmed that the first set of material parameters is randomized with a CoV value of 10%.

The detailed distributions of material properties are further used to determine the global material parameters of the slope. These global material parameters are marked for each node of the benchmark heterogenous meshed slope. These global material parameters are further used to determine the global elastic constitutive matrix (C), global pseudo-time stepping (dt), and global stiffness matrix (K). These global parameter matrices are also distinctly defined for each node of the slope to effectively map spatial probability distributions of material properties over the slope domain. These parameters are used to prepare the Mohr-Coulomb elastoplastic slope model and are further used to calculate the deformations, strains, and stresses globally. This requires solving a large number of variables and simultaneous equations in loops. A matrix-based approach to the MATLAB numerical tool can enable and manage complex calculations, providing a suitable platform for better user-control analysis and representation of the slope system.

3.2.   Detailed deformations and failure analysis

In this probabilistic finite element elastoplastic analysis, the testbench heterogenous slope is considered with the inclusion of the first set of randomizations in material properties (uniform pdf and 10% CoV), as described in previous Section 3.1. The initial state of global deformations, stresses, and strains is determined considering the elastic nature of the slope. After that, the slope model undergoes a step load increment for a fixed number of 100 steps, and corresponding deformation changes, stress changes, and strain changes for each node of the slope space are assessed. Moreover, the stress states are used to check the Mohr-Coulomb failure criterion for each gaussian point of the element. In this case, the slope demonstrates elastic characteristics, leading to the calculation of corresponding deformations, stress changes, and strain changes.

If the failure condition is met, then excess stresses and plastic strains are calculated when considering associative flow conditions. When the failure state materializes, the slope begins to exhibit plastic behavior, leading to the calculation of plastic strains. Then, in the process of overestimating stress distribution, plastic strains determine the excess load, and the next loading sequence utilizes this surplus load. When the stress state only temporarily violates yield, plastic straining occurs. The value of the failure function, which expresses how much the current stress state surpasses the yield condition, determines the visco-plastic strain rate. As a result of increasing visco-plastic stresses over time or through iterations, the material relaxes, and the yield function and visco-plastic strain rate decrease. When the visco-plastic strain rate moved towards a very low value or convergence, the incremental plastic strain and stress change, respectively, were equivalent to the accumulated visco-plastic strains and the corresponding stress change. Moreover, the nodal forces accumulated due to plastic strains are determined and used in the next iteration until convergence.

The slope model is simulated for 100 steps with a step load increment of -0.1 kN/m, ranging from 0 kN/m to -10 kN/m. The maximum of 200 iterations can be used to adjust the excess force generated due to stress distribution while full yielding occurs. The top layer of the slope's nodes receives the applied load. The Mohr-Coulomb elastoplastic model is employed as a constitutive slope model, and deformations for each node are calculated and stored in two arrays of size 578*1, one for X-direction deformations and the other for Y-direction deformations.

Figure 8 represents the deformations of all nodal points over the slope space, obtained during the numerical simulation of the testbench slope with 10% CoV. The maximum deformation in the +X and -Y directions is 0.2129 m and 1.5634 m, respectively. Also, the minimum deformation in the +X and -Y directions is -0.097 m and 0 m, respectively. The frequency distribution plots of both deformations were also obtained during analysis and depicted in Figure 9, which gives exact information on the number of nodes displaced in the X-direction and Y-direction of the XY-plane, where the slope is situated.

A three-dimensional matrix for the stress tensor is acquired from the simulation of a MATLAB based benchmark slope model. The size of the stress tensor matrix is three dimensional and 4*4*528. The stress matrix stores all four stress components (sigma_xx, sigma_yy, sigma_xy, and sigma_zz) for each gaussian point in a row. Each gaussian point in an element is associated with four rows. This 4*4 arrangement sub-matrix is spread in the Z-direction, representing the stress information of each individual element. Figure 10 presents the systematic representations of different axial and shear stress profiles. Figure 11 also depicts the frequency distribution plots for these stress components, determined for critical slope analysis. Similarly, a MATLAB-based slope model simulation determines the strain tensor, and Figure 12 presents axial, shear, and volumetric strain profiles.

The simulation of the two-dimensional slope model under plane strain conditions results in a relatively long slope in the spatial Z-direction, leading to the assumption of zero axial strain in the zz-component. The strain profiles comprise strains in xx-component, yy-component, xy-component, and volumetric strain. Figure 13 also depicts the frequency distribution plots for these strain components. In this frequency distribution, strain in the XX-plane having both positive and negative values means nodes are exhibiting tensile and compressive horizontal deformations, while the strain in the YY-plane is negative, meaning nodes are facing a major compressive effect.

The simulation results also provide information on failure nodes and where the slope yields. The distribution of these nodes, where a yield function fails or nodes enter the plastic region, is presented. Figure 14 plots the spatial distribution of failure nodes, or potential yielding regions, over the benchmark slope and presents these nodes with red asterisks. The recorded data correspond to a heterogeneous testbench slope with homogeneous material uncertainties of 10% CoV when the slope was subjected to a step load of -10 kN/m (100th step load sequence) at the top.

The benchmark slope yields a total of 231 to 235 failure nodes, with corresponding nodal failure percentages of 39.96% and 40.65% under step-load values of 2 kN and 200 kN, respectively. Figure 15 presents the variation of the nodal failure percentage with the step-load sequence, taking into account a 10% CoV. When a total of 100 step-loading sequences with a step increment of -0.1 kN/m are applied to the top edge of the benchmark slope, the percentage probability of nodal failure is 40.646%.

In this elastoplastic analysis, each failure node enters the plastic region a number of times. The occurrence of a node in the plastic zones depends on the access nodal stress generated due to self-weight load, external load increment, and residual force exerted against visco-plastic stress distributions. The visco-plastic excess stress is distributed using the constant stiffness method. In this iterative analysis, the frequency of yielding of each of the failure nodes is recorded and plotted as depicted in Figure 16.

The variation in frequent appearance in yielding can be utilized to define the severity level of the failure zone. The high value of the failure frequency of a zone represents the more sensitive failure region. Based on the normalized failure frequency values, Table 2 categorizes failure nodes into four groups and divides the overall yielding region into four zones. Table 2 also shows and lists the number of failure nodes in each zone.

The spatial distributions of zone-wise failure nodes are represented on the testbench slope, as shown in Figure 17. The failure nodes of zones 1, 2, 3, and 4 are represented by blue- dots, black- crosses, green asterisks, and red triangles, respectively, in Figure 17.

3.3.   Effect of different coefficient of variation (CoV) values

The developed slope model is simulated for different values of %CoV using MATLAB. Here, six randomized sets of material properties are generated, characterized by %CoV values of 0%, 5%, 10%, 20%, 30%, 40%, and 50%, respectively. The top edge of the testbench slope is subjected to stepwise external loading of 100 steps with a uniform load increment value of -0.1 kN/m. The changes in %CoV of uniformly randomized material properties are used to analyze the corresponding variations in slope deformations, stress limits, strain limits, and failure characteristics during the statistical analysis of the testbench slope. All simulation results are investigated using upper and lower bound values under finite element limit analysis for different values of CoV.

The boxplot of horizontal and vertical deformations for all six values of %CoV is presented in Figure 18. It is observed that the maximum and minimum values of horizontal deformation are shifted upwards (+X direction) without significantly changing the median value, and the lower bound value of Y-direction deformation is reduced as CoV increments.

Figures 19 and 20 also record and represent the optimal points of stress and strain components, with varying values of %CoV. It is observed that the upper bound values of all axial stresses remain unchanged, but the upper bound value of shear stress is notably changed as the CoV increases. Additionally, an increase in CoV has the potential to increase the lower bound values of all stress components. The axial strain component in the YY-plane remains nearly unchanged as CoV increases, while the shear strain component experiences the maximum changes. The maximum percentage change is observed in shear components, so it clearly implies that the benchmark slope is much more susceptible to shear failure.

The failure characteristics parameters, nodal failure, and probability of nodal failure are also investigated for different values of %CoV in the spatial probabilistic finite element analysis. The variation of %nodal failure with step-loading (step-load increment of -0.1 kN/m at the top edge of the slope) is investigated for deterministic material values (CoV = 0) and six different values of CoVs; thereafter simulation results are plotted, as shown in Figure 21. The nodal failure is increased with external load increments as well as a rise in CoV values.

The probability of nodal failures versus %CoV is also plotted using a bar chart in Figure 22. The probability of nodal failure is also increased as CoV increases. It is also observed that the nodal failure pattern and probability of nodal failure are not significantly changed for CoV values of 0% and 5%. A sudden change in nodal failure pattern and probability of nodal failure is identified as the %CoV crosses the value of 20%, as shown in Figures 21 and 22.

4.   Discussion

We explore the behavior of heterogeneous elastoplastic slopes and aim to provide insights into the deformations and failure characteristics using the nonlinear spatial probabilistic finite element method under self-weight and step loading sequences considering different levels of material uncertainties.

4.1.   Deformations and failure analysis

During the numerical simulation of the first randomized material properties set for the test bench slope (uniform pdf, 10% CoV, and with the application of 100 load steps), detailed deformation profiles are determined, which provide exact spatial displacement information in both dimensions at nodal points and are cumulatively represented in Figure 8. The deformation profile of the benchmark slope has a maximum deformation of 0.2129 m in the X direction and exhibits a radial deformation pattern near the unconfined slope surface. This radial deformation pattern is generated due to slope confinement and restrained movement conditions. It is also observed from the frequency distribution plot of deformation that the number of nodal displacements in the positive X-direction is very few (<50), and none of the nodes are displaced in the upward Y-direction. In the same way, the maximum deformation of 1.5639 m in the downward y-direction is observed; hence, it presents an irregular tensile-compressive type deformation pattern. This behavior is due to the restraining of Y-direction movement in the bottom and self-weight and step loading in a downward direction only. It is noteworthy that the toe area lies in layer 2 of the slope, which has low cohesion and friction angle of 9.8 kPa and 5^o, respectively. Thus, the slope exhibits the maximum density of failure nodes in layer 2, as presented in Figure 14. The variations of tensile-compressive stress and shear stress spread over the slope surface are observed due to the heterogeneity and spatial randomization of the slope. The optimal stress components are determined by values of -275.09 kPa, -607.86 kPa, 72.50 kPa, and -276.85 kPa in the XX-plane, YY-plane, XY-plane, and ZZ-plane, respectively. The spatial region of higher and positive values of shear stress is mostly confined by horizontally restrained side boundaries, and some parts are near the toe region of the slope, so the common area of layer-2 and effective positive shear stress values is the most prominent slope failure region. Also, the dominant strain component is shear-strain in the XY-plane, with a maximum value of 367240 µε. The dominant shear strain, downward axial stress, and positive shear stress indicate that the slope will undergo sliding in the shear plane with a low value of tensile-compressive failure, so it is reasonably unstable. Although it is validated by the percentage probability of nodal failure value of 40.646% for the benchmark slope problem, which exhibits an unstable slope condition. A total of 232 failure nodes are observed for the 100^th step-load sequence on the testbench slope. The failure region is classified into four zones based on normalized nodal failure frequency, as depicted in Figure 16 and Table 2. The high value of normalized failure frequency represents a more severe failure zone over the testbench slope, as represented in Figure 17. Inclusion of failure frequency information with failure nodes provides potential findings for studying crack propagation by numerical methods.

4.2.   Effects of uncertainty levels

The level of uncertainty around deterministic values of material properties is attributed to CoV. The testbench slope behavior is investigated for the seven different values of %CoV ranging from 0% to 50%. The finite element limit analysis is used to investigate deformations, stress components, and strain components for seven different values of %CoV and is presented in Figures 18–20. The maximum changes in horizontal deformation limits are observed at 50% CoV, and the corresponding percentage changes in the lower and upper limits are observed at -29.52% and 7.17%, respectively. Therefore, the maximum change in deformation is observed at the lower limit of horizontal deformation. The fusion of the deformation pattern in Figure 8 and the boxplot of deformation presented in Figure 18 indicate the enhancement of swelling of the toe area and shrinking of the crown area of the slope surface with an enhancement of uncertainties in material properties.

In the same way, lower limits of all four stress components are effective, varying in proportion to the CoV. However, the maximum %change in stress is observed at 96.57% in the lower limit of the stress component in the XX-plane, for the %CoV changing from 0 to 50%. This is due to the rise in uncertainty in the material properties of the slope, which is confined by the slope's side boundaries, restraining horizontal movement at both side edges. The magnitude of the upper and lower limit values of shear stress components increased with the rise in uncertainties, as shown in Figure 19. It is observed that the rise in upper limit values (positive values) is comparatively larger than the lower limit (negative values) of shear stress, which means the uncertainties enhance the shear sliding or slope instability. This is also validated by the rise in probability of nodal failure with an increment in CoV, as shown in Figure 22. The nodal failure-step loading sequence curve for different values of %CoV provides a potential insight into uncertainty analysis, as shown in Figure 21. The nodal failure patterns are not affected much by the rise in %CoV from 0% to 5%; this indicates that the lower uncertainties (%CoV under 5%) do not affect the failure behavior for the testbench slope. Here, a sudden rise in the nodal failure curve is observed when the %CoV is equal to or more than 30%. Also, the abrupt rise in probability of nodal failure is recorded for %CoV values changing from 20% to 30%, as represented in Figure 22. This study indicates that if the uncertainty level in material properties crosses a certain level, the deterministic solutions cannot be considered for stability analysis, and there is a strict need for the inclusion of uncertainty analysis for correct and effective instability analysis.

4.3.   Limitation and future work

Further refinement of the numerical models and incorporation of more sophisticated soil constitutive models can enhance the accuracy and reliability of slope stability measurements. Conducting field investigations and validation studies can improve the numerical analysis's findings and provide real-world data to improve model calibration and verification. Effective slope stability analysis can also be used to preserve important physical properties and structures ^[68,69]. Additionally, exploring multi-scale approaches that integrate information from laboratory testing, field observations, and numerical simulations can provide a more comprehensive understanding of slope behavior across different spatial and temporal scales. The present modeling and analysis are also used for a set framework for Monte-Carlo analysis and sensitivity analysis of slope properties, including multi-level spatial material variabilities and spatiotemporal randomization, which can be further used in reliability and risk analysis.

5.   Conclusions

We employ modeling and analysis of advanced numerical techniques to simulate the complex behavior of heterogeneous slopes using the MATLAB platform. Nonlinear PFEM accounts for both material and geometric nonlinearities and probabilistic aspects effectively in slope stability analysis. The developed slope model incorporates the quantitative spatial uncertainties in terms of %CoV in slope material properties using PFEM. After that, numerical simulation is successfully performed, and detailed information about deformations, stress patterns, strain patterns, possible pre-failure zones, and failure features of heterogeneous slopes is acquired. Additionally, two novel quantitative parameters are introduced and computed to shed more light on failure investigations: nodal failure and probability of nodal failure in this probabilistic finite element analysis. Overall, this study offers a comprehensive framework for more realistic elastoplastic modeling and a novel approach to analyzing the detailed deformation and failure characteristics of a geotechnical slope. The main conclusions for the present research work are as follows:

● The failure region is mostly found in layer 2 of the benchmark slope, as it has a low value of cohesion and friction angle.

● Considering 10% CoV in all six material properties, the benchmark slope's horizontal deformation profile displays a radial deformation pattern close to the unconfined slope surface, with a minimum and maximum deformation of -0.099 m and 0.2129 m, respectively.

● The testbench slope exhibited an uneven tensile-compressive type deformation pattern due to the maximum deformation of 1.5639 m in a downward y-direction for 10% CoV. The whole Y-direction nodal deformations are in a downward direction only. This behavior is caused by self-weight, step loading in a downhill direction exclusively, and restriction of Y-direction movement at the bottom.

● The extreme value of stress components is -263.6 kPa, -600.9 kPa, and 71.4 kPa in XX-plane, YY-plane, and XY-plane, respectively, and the extreme value of strain components is 8700 µε, 3400 µε, and 36000 µε in XX-plane, YY-plane, and XY-plane, respectively, for the first set of randomized material uncertainties of 10%.

● For every load step, the percentage nodal failure of the test bench slope of 10% CoV is calculated and varies from 39.69% to 40.13% for 100 successive step loading sequences. The failure region is further divided into four zones based on the normalized failure frequency to facilitate a more thorough and comparative failure analysis.

● Maximum deviations are observed in the lower limit of horizontal displacement, while varying material uncertainties exist for the testbench slope. The %changes in the lower limit of the horizontal deformation with respect to its deterministic boundary value are -2.38%, -4.88%, -10.27%, -16.17%, -22.58%, and -29.52, for %CoV values of 5%, 10%, 20%, 30%, 40%, and 50%, respectively. This indicates that the X-direction deformation decreased as the %CoV increased.

● The upper and lower bound values of the shear stress components also increased in magnitude as material uncertainties increased. Furthermore, the increase in the upper limit values (positive values) of shear stress is significantly greater than the decrease in the lower limit values (negative values), indicating that the rise in material uncertainties amplify the shear sliding and slope instability.

● The failure characteristic parameter, the percentage probability of nodal failures, is increased with an increment in %CoV. The percentage probability of nodal failures for the testbench slope is observed at 39.86%, 40.46%, 40.84%, 41.78%, 41.93%, and 42.09% for %CoV values of 5%, 10%, 20%, 30%, 40%, and 50%, respectively.

● The sudden changes between nodal failure curves of different %CoVs are used for defining an acceptable threshold of uncertainties and for considering deterministic solutions. Deterministic solutions for the testbench slope cannot be taken into consideration if the amount of uncertainty exceeds 20%.

It can be summarized that by integrating probabilistic methods and finite element methods into slope stability analysis, engineers can estimate pre-failure precisely, analyze detailed deformations, stresses, and strains efficiently; thus, they can effective decisions, manage risks and mitigation efficiently, and enhance the resilience of slope infrastructure in geotechnical projects.

Author contributions

Peeyush Garg: conceptualization, data curation, formal analysis, methodology, software, writing–original draft, writing–review & editing; Pradeep Kumar Gautam: conceptualization, formal analysis, methodology, validation, writing–original draft, writing–review & editing; Amit Kumar Verma: formal analysis, writing–review & editing, supervision; Gnananandh Budi: formal analysis, writing–review & editing, supervision. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

We would like to extend our sincere gratitude to the editors and reviewers for their invaluable feedback and guidance throughout the review process.

Conflict of interest

The authors declare no conflicts of interest.