Stability estimation of slopes with determined slip surface by the MSTAB-3D method based on sliding body limit equilibrium analysis

The method developed for this study, established on the premises of the limit equilibrium flat analysis for a spatial solution, is a modification of the STAB-3D method, previously described by the author. It combines the analyses methods of 2D slices of flat cross–sections with the spatial analyses methodology rooted in a specific breakdown of a landslide sliding body into 3D elements assuming some simplifying solution. However, this method is solely applicable in case of a landslide failure with a stipulated slip surface and with a consistent decline of a determined slide direction. Such a method was developed in the article published earlier, which provided then its basic assumptions and the equilibrium formulations. The following publication thereof, presents overall suppositions for this method as well as its modification involving the resultant forces brought to the equilibrium with the generalized slide direction. Apart from that, a comparative analysis was carried out on the impact of this modification applicability of the obtained results with regard to the STAB-3D method. The algorithm was also presented concerning the modified method with its results being compared to a couple of selected methods LEM (limit equilibrium method). The undertaken analysis reveals that the modified MSTAB-3D method determines stability indicators that are very similar to its earlier version. Moreover, the results occur to be also approximating the values obtained in the course of other methods with regard to the flat cross-section analysis.


Introduction
For decades, attempts have been made to create foundations of the 3D analysis and new analysis methods of spatial slip bodies are being developed. The solutions discovered so far pose many limitations, they sometimes undergo multiple modifications and some cases lack successful outcomes. Among others, Hungr [1], Hutchinson and Sarma [2], Michalowski [3], J. L, Chen and Chemeau [4] and Zhang [5], to name only a few, have been dealing with this issue for over thirty years and have been trying to take the three dimensional character of the stability analysis of slopes into consideration. Cheng et al [6], Griffiths and Marquez [7], Huang and Tsai [8,9], Hungr [10], Kalatehjari et al. [11] and Wei et al. [12] belong to a more recent group of authors investigating that area of knowledge.
Apart from the widely known and utilized flat methods of slope analysis, it is possible to distinguish a collection of numerical methods comprising computer applications stemming from a solution grounded in the theory of continuous media mechanics applied to fragmented materials such as soils. Computer numerical programs of soils analysis are based on the assumptions concerning the following methods: the Finite Difference Method, the Finite Element Method, and the Boundary Element Method. Among numerous publications on 3D analysis, the works presented by e.g. Hicks and Spencer [13] and Li et al. [14] and Nian et al. [15] and Zhang Y. et al. [16] shall be mentioned. Most of them can serve to depict both the stress value resulting in the stability loss as well as the stability checking of slopes or embankments. The most widely applied algorithm among the numerical methods concerns gradual reduction in shear strength and analyzing the change in stability values until the instability moment is reached, i.e. obtaining the factor of safety FS=1 (Shear strength reduction technique-SSR). The other method (Modified shear strength reduction technique MSSR), is a modification of the SSR one which enables to determine a couple of potential slide routes instability lines in case of a complex slope shape as well as varied soil layering, in contrast to the SSR method allowing for depicting a single line only.
A reversion to the proven 2D methods is also observed through incorporation of their strict analytical values 3D and for the sake of results verification. Cheng and Yip, among others, in [17] included a novel approach to the 3D analysis of a slope stability expanding the possibilities of mature and reliable 2D methods applied so far: Bishop's method or Janbu and Morgenstern-Prica's method. He formulated the assumptions for stability in three-dimensional space where it is possible to indicate the main slide direction of soil mass as well as the sectional one. Such a method of analysis is indicated for utilisation in more complex soil layering conditions. Thus, it is possible to apply the familiar and reliable 2D methods implemented for the flat slices for the sake of 3D analysis in the entire sliding body space. Similar conclusions arise from the article written by Loehr et al. [18] which describes the quasithree-dimensional method based on conventional, twodimensional stability analysis of slopes, defined as the resistance-weighted procedure (RW). It is an approximate method resulting in slightly overrated values of stability indices by about 10% with regard to the remaining 3D methods. Nevertheless, Michalowski [19], who conducted studies according to the kinematic approach to boundary analysis, has drawn conclusions determining the dependence of slope stability indices by comparing 2D approach to the 3D one. During estimating a 3D spherical weakness area failure shape, he proved that the 2D analysis method of the same slide data leads to obtain lower safety factor values than in case of the 3D method for variable proportions of dimensions and outer loads affecting the slopes. These discrepancies, depending on the width/height proportion ratio of a sliding body, may reach as much as 50% (with the above ratio of 1/1) or 10% (in case of 5:1 proportion ratio). Huang and Tsai [8] presented an interesting method for the 3D and asymmetrical stability analysis of slopes. It based on the assumptions developed by Hungr [1] and later by Hungr et al. [10]. It was further expanded and complemented by Huang et al. [9] and it renders basic integrals for a spatial analysis of a sliding body within landslides. It takes into account two assumed directions of movement: parallel and perpendicular one. This method enables to verify the pre-defined direction of landslide failures besides providing the classic balance values of strength and torque in 3D conditions. Similar assumptions underlie the MSTAB-3D method. However, due to the limitations of that method for landslides with stipulated movement directions (e.g. for occurring events as well as for the pre-determined land failures in case of slopes with geological features enabling concise determination of the landslide direction), it was possible to simplify that complex assessment. This simplification relies on rendering the 3D column positioning unnecessary, thus limiting the analysis to a whole sliding body calculations in sections parallel to the landslide movement direction. The reason for such a presumption is the fact that in case of the direction being perpendicular to the slide, the resultant force integrals are in equilibrium for the mentally shredded strips of cross-sectional block of elements ( Fig.1 Block Element 3-Dimensional BED) throughout the entire landslide body. It is proven by the on-site observations of land failures having only one direction of slide whereas they are stable perpendicularly to it.
The methods presented above determine the stability conditions of a Landslide Sliding Body (LSB) which may be applied less or more adequately to the processes undergoing the landslide occurrence. The majority of them perceive LSB as a kind of continuum material with constant strength parameters and with the stability features remaining intact if the complex safety indexes as well as specified soil parameters are taken into consideration. Meanwhile, the studies and observations generate results which demonstrate that some landslide failure once activated reveals considerably variable dynamics (from slow motions to stabilization periods and sometimes violent accelerations). That variability stems from the influence of soil factors inducing its decomposition and strength parameters change such as water pressure and buoyancy force or earthquake. Usually, the disturbance of a stable structure within the LSB of the already activated landslide result in perturbation of consolidation and causes the internal soil parameters reductions inside LSB. Therefore, the assumptions of e.g. numerical methods derived from a solution accepting the theory of continuous media mechanics may occur inadequate in case of the already activated landslides. Even after the application of some modification like the SSR method, the simplification may only be partially eliminated. By contrast, the solutions basing on the stability estimation of LSB overall limit equilibrium reflect its stability state in a more reliable way.
Due to the fact that the analysis methods concerning slope stability which have been applied in recent years and constantly modified or improved, still reveal numerous limitations and are very complex, a simplified analysis method of quasi-spatial equilibrium of LSB has been developed, called shortly MSTAB-3D (Modification STABility and 3 Dimensions). It is a modification of the STAB-3D method, previously described in [20]. Its main surmises are also discussed in [21] as well as in the present article. It has been practically evaluated and explored during resuming control over the landslides occurring within the area of Poland. This method, grounded in the familiar and elementary issues of soil mechanics and statics laws, enables to consider all the mentioned statics and strength elements simply and clearly and it proved to be effective in real life practice while inhibiting the development of a couple of active landslides. It explicitly concerns the cases where there is a need to counteract the already occurred, slowly moving landslides being alternatively in states of unstable, temporary balance and cyclical activations.

Description of MSTAB-3D method 2.1 Overall assumptions of the method
MSTAB-3D method belongs to a new approach on the issue of LSB stability through the prism of 3D spatial view in relation to landslides of structural type (Structural Slope Failure -SSF) after the determination of the slip surface (when the soil layering is explicitly specified according to its slide surface having a consistent decrease value and verified direction of the slide). That method reflects all the main factors occurring in nature which significantly affect the analysed stability of LSB such as: hydrostatic and groundwater flow pressure, external stresses, restraining forces induced by the application of supportive equipment and anchoring as well as a possibility of considering the forces resulting from the impact of vibrations generated by e.g. vehicle transportation. It enables to designate the factor of safety of either steady slopes or embankments which reveal symptoms of stability loss or active land failures. It may also be utilised for establishing the extent of landslide movement within stabile rock masses zones.
With regard to the earlier version elaborated in detail in the article [20], this method was slightly modified. The divergence between both methods assumptions lies in a different summing mode of resultant forces in the discretized section slices. By contrast to the more simplified, earlier version, the vector aggregation was applied both for restraining forces as well as for sliding ones responsible for the landslide failure. Simultaneously, the resultant forces were approximated to a common direction line plotting them on a mutual, differentiated direction. Such a line was chosen by assigning a general decrease line of the landslide slip surface in the slice which was the closest to the main axis of the landslide sliding body. The other hypotheses of MSTAB-3D are in analogy to the earlier ones in STAB-3D version and they are shortly delineated beneath including the introduced changes and uncovering its new algorithm.
The so-called "MSTAB-3D" method belongs to the Limit Equilibrium Methods (LEM). It only accounts for structural type of slide failures (SSF) after determination of the slip surface and it allows for factor of safety calculation for the entire three-dimensional sliding body in a simplified manner. The simplification of the analysis derives from a reflection of the fact that a sliding body as a whole is internally balanced in the direction perpendicular to the slide. The internally appearing pressures and sliding forces perpendicular to the slide direction are excluded from the stability analysis as they are in balance. On the other hand, the equilibrium in the parallel direction is analysed where the overall, calculated result is always different from zero. Hence, the division of the landslide sliding body LSB is concluded into elements discretized in vertical flat sections parallel to the slide direction (Fig.1).   Fig. 1. The view of a slope and sliding body with its division: 1-outer boundary of LSB sliding body, 2cross-sections (CS1÷ CS8) parallel to the slide direction, with BED elements lying between them, 3layering of the advantageous slip surface.
The hypothesis of the method allows to suggest that all potential stresses stemming from self weight as well as from the outer impact (buoyancy force, hydrostatic and water flow pressure, car and rail transportation load, dynamic and seismic vibrations, buttresses, anchoring, supporting constructions etc.) in the form of resultant forces applied to each discretized integral elements (BED) incorporated by the entire sliding body are being taken into consideration.
Moreover, additional suppositions were also accepted and they concerned the following: a) a sliding body was divided into K BED elements, by creating m flat, vertical cross-sections CS parallel to the slip direction, what is shown in Fig.1. b) each CS cross-section was discretized into n flat, vertical slices BFE (Fig. 3). c) maximum overall angle of inclination g was designated for the slip surface (Fig. 3) in CS sections, which constituted a reference point for balancing the sliding forces as well as the restraining ones both within each CS and BED, d) it was assumed that all vectors of resultant forces Zk (Fig.2), operating perpendicularly to the slide direction in each BED element, after adding vector values in the entire LSB sliding body, result in 0 values.
e) it was assumed that shear resistance is allocated to the vertical walls of outer BED and its resultant force value is higher than the resultant values of the operating shear stresses and therefore, the sliding body stability is affected by resultants of the sliding forces and restraining ones for the entire sliding body, f) on the sliding coating of the contact surface (the bottom of BED), soil shear resistance shall be observed which is dependent upon soil cohesion as well as from soil angle of international friction of this coating which induces the force responsible for BED elements integration, g) internal interactions between BFE slices are balanced within the cross-section CS, thus, they are not taken into account in the stability analysis, h) the resultants of the sliding forces as well as the restraining ones derived from the internal forces balance present in the cross-sections CS (Fig. 4)  It must also be highlighted here that this methodology premises shall not include research on landslides with round-cylindrical nature of slide surface (circle sliding lines in 2D cross-section) which usually prevail in non-cohesive, homogenous soils. The method discussed in this article accounts only for the cases when an explicit slide layering was formed at the interface of soil layers. It is supposed to be the weakest area of balance which is responsible for the sliding body shape. It is very improbable, however not impossible, for such a sliding coating to be created of itself in a cylindrical shape. Such an occasion would impose a different research technique to be applied since the condition of a rotational movement of a sliding body should be taken into account. The applied MSTAB-3D method complies with a guideline rendering that the sliding coating takes on a shape different from a round-cylindrical one and only the condition of balance concerning the sliding body for the displacement over such a layer becomes significant in the study. Since this issue belongs to one of the basic assumptions of this method, it was possible to exclude stability torque analysis from the presented method.

Determination of MSTAB-3D method preliminary parameters
According to the method under question, prior to the equilibrium analysis, preliminary and preparatory works shall be conducted for the sake of specifying the calculation model. In case one encountered a slope or embankment where the probability of land failure is significant or the landslide is already active longitudinally, a suitable geotechnical and engineering recognition must be carried out basing on the undertaken drillings or potential exposition of soil layers. The following step is connected with determination of a potential zone of an anticipated land failure requiring stability analysis. At this phase, the application of georadar technology is helpful as it is intended to screen out the soil density variations and it may be supported by the Geographic Information System (GIS) and explicitly discussed by Xie et al. [22] and Shen et al. [23] through the prism of applicability issues in such studies. These are very useful systems especially in the zones with active landslides where it is not possible to conduct any drillings or excavations.
A contours plan of a slope as well as the slide plate shown in Fig.1 are prepared for the studied area. Delineation of the surface waters flow conditions and groundwater interaction within a landslide sliding body and its bedding constitute vital elements of the preparatory recognition. It enables to consider buoyancy, hydrostatic and waters flow pressure which can also be outlined in a form of contours plan. Geometrical parameters are reported both for the sliding body itself and within its direct surroundings. The designed spatial model is intersected with vertical, parallel planes which mostly are located at equal distances to each other along entire potential failure surface of the sliding coating ( Fig.  1 -cross-sections CS1÷CS8). The following step involves plotting of the slope or embankment shape, land failure line and ground water table line onto these planes (Fig.  3). The enumerated stages are the preliminary bases for the static analysis of the slope or embankment with the MSTAB-3D method.

Fundamental correlations for the MSTAB-3D method
Tab.1 presents a detailed algorithm of proceedings aiming to delineate the safety factor FS according to the MSTAB-3D method. Following the preliminary preparations discussed in point 2.2., it starts from distributing a distinguished model of a sliding body into spatial elements BED with the width of zk. They are limited by CS sections which are parallel to the drop direction of the slip surface and they constitute vertical side walls for each of the elements. However, the crosssections are divided into flat slices with vertical walls and are called Basic Flat Elements BFE, by analogy to the manner presented e.g. in [24].

k = K k < K
Stability analysis starts from defining the forces operating on a basic BFEi element as well as their balance allowing for resultant forces delineation. One of the major forces present in each case is the soil weight Pi that ( Fig. 3 and 4 respectively) may be calculated according to the following formula: That force is located in the centre of gravity of BFEi, and i is a value approximate to a unit weight of soil which forms the analysed slice. External load forces such as: abutments, anchoring or soil cover stresses etc. are represented here by an integral vertical value Di and a horizontal one Vi. To be able to reflect the additional external vertical force Di operating on a slice, an aggregate Gi force value is calculated, in line with the following formula: (3) Within a slice, underneath the ground water table, it is necessary to take its buoyancy into consideration which is determined by its force Wi w calculated by a formula (4) in which w represents the unit weight of water.
Hence, the total resultant vertical force Gi w that operates on the slice can be determined by the formula (5): Additionally, also a horizontal external force Vi may interact with the BFEi slice which may be divided into its integrals: Vpi -parallel to the slide course (inclined at an i angle to the horizontal line) according to the formula (6) and perpendicular to the slide course Vsi according to the formula (7): The hydrostatic water pressure value is determined by a resultant force Wi, which is calculated on the basis of the following formula: The value resultant to all the forces within the slice and perpendicular to the slide course Ni is computed by the formula below: Calculation of Ni force which affects the intensity of friction enables to determine the value of shear soil resistance of the contact coating in the formula (12).
Apart from the forces described in formulae (2 ÷ 9) the forces presented in Fig. 4  The next proceeding involves indication of the resultant sliding force Hi o by aggregation of all the total forces operating in a parallel direction to the slide course (10). Then it is plotted on the slide course along the overall angle of inclination of the slide surface g in the formula (11). Consequently, the resultant sliding force Ti o is determined as a sum of all the integral forces operating along the direction of slide (12), and then it is plotted onto the sliding course along the overall angle of inclination of the slide surface g in the formula (13).
The sliding and restraining forces calculated on the basis of the (11) and (13) correlations for single slices of BFEi allow for generating the total sliding and restraining forces for entire flat cross-sections CSi. In order to receive the values of sliding and restraining forces Hk and Tk representing the entirety of BED elements with the width of zk (Fig. 2b), one needs to work out the mean values from two adjacent cross-sections CSj and CSj+1 (Fig. 2a) The result value in a form of the stability factor obtained in that way is characteristic for the state of the entire landslide sliding body. It is also possible to designate stability factors in a similar way for both elements of a landslide sliding body: the selected ones as well as separately distinguished flat sections CS.

Comparison of MSTAB-3D to other selected methods
For the sake of verification of the discussed method, the results obtained in the course of a few methods application were collected in table 2 and 3 and presented for comparison. A well documented active landslide failure at the opencast working levels in KWB Bełchatów coal mine presented in Fig. 5 was the study subject. Initially, calculations were carried out for a few methods of flat landslide section analysis. Here, limit equilibrium methods LEM were chosen that correspond to the discussed method under question but at various levels of simplification. The findings presented in Table 2 below stemming from those calculations highlight a slight divergence in the value of stability factor, between 0,5÷0,7%, for Masłow-Berer's method, Szachunianc's method and the earlier version of STAB-3D. Whereas, the result value obtained in the course of applying the Large Blocks Body method based on misleading, simplifying assumptions is safer and lower by over 7%. The second attempt reflected generating a comparison of the STAB-3D and MSTAB-3D methods in a spatial approach for analysing a sliding body as a complete entirety. Its results are outlined in Tab. 3. The compilation clearly indicates that the result value obtained for a landslide sliding body with the use of the MSTAB-3D is lower by 0,2286% than the one generated by the STAB-3D method. Therefore, modification of the method which approximates real internal correlation within the landslide sliding body is practically equally accurate as its earlier developed version. On the other hand, the comparison of results deriving from the analysis of the flat approach generates increased indices values with regard to its spatial solutions counterpart.

Conclusions
Research shows that the LEM methods used till now offer the most reliable results especially with regard to the longitudinal, active landslides. Consequently, there is a need to develop and introduce methods adapting the solutions of LEM 2D methodology into the spatial approach 3D. The following article outlines a solution called the MSTAB-3D method which is a modification of the STAB-3D method described in [20].
The method under question may be characterised with the possibilities and benefits such as: (i) it enables to carry out an analysis of balance between the major forces affecting the landslide sliding body with the use of a non-complex algorithm of calculations, (ii) it allows for determination of the safety factor concerning an intact slope or a landslide failure which is already activated with a possibility of taking the elements supporting its stability into account. However, that method may only be applied in case of a rock mass topography characterised by determinate slip surfaces and the slide direction is explicitly determined. Such a limitation does not pose any exclusion criteria as far as its practical applicability is concerned due to the widespread pervasiveness of such cases as e.g. in flysch soils or in contact surfaces of rock layering, silt or clayey soils. It is based on the analysis of 2D cross-sections intersecting the entire landslide sliding body into spatial elements enclosed within the vertical planes which are parallel to the direction of the sliding plate descent. Moreover, these sections are divided into slices with vertical walls. This assumption permits not only to disregard all the forces which are perpendicular to the direction of slide in a landslide sliding body, what was proved accurate by Hungr [1] but also to indicate sliding and restraining forces present only for one direction. This hypothesis leads to a simplification of the issue without any danger of committing a practically inaccuracy in the sphere of precise solutions generating.
The MSTAB-3D method differs from the STAB-3D one with the calculation mode of the safety factor for the entire landslide sliding body, which involves approximating the resultant sliding and restraining forces within the slices to a generalised sliding line. This line represents an overall inclination of the slip surface of the least stabile section. For the sake of results verification obtained during the MASTAB-3D methodology application, calculations concerning an exemplary slope were conducted which underwent a landslide failure and the occurrence was later documented explicitly. The comparison of the discussed method with the other LEM methods and with its earlier counterpart has enabled to underlie a conclusion that both methods the old version and its modification identify lower values of safety factor in 3D analysis then in 2D study of the most disadvantageous section. Nevertheless, the divergence between both methods in 3D analysis is minute and is about 0,2%. Moreover, the modified MSTAB-3D method introduces significantly homogenous values of the safety factor with regard to the remaining analysed ones (it constitutes a value of 0,5÷0,7%) while comparing their results for 2D section analysis. It is also possible to discern that analysing the landslide failures through the 3D prism applied e.g. in the MSTAB-3D method, accedes to obtain results reflecting the spatial character of them which is considerably different (about 10%) from the analysis of a selected single flat section.