A New Insight into the Stability of Precariously Balanced Rocks

Large granitic boulders resting on steep slopes represent considerable safety hazards that largely depend on the location of the contact surface characterized by the impression d, denoting the parallel distance between the contact surface and the original rock surface. On the other hand, this impression reflecting the often convex nature of the contact between boulders and resting platforms, cannot be measured precisely, so Factors of Safety (FoS) computed with this input may have significant uncertainties. Using geometric 3D analysis, here, we present the concept of computing FoS as a function of the impression d, admitting a much more reliable estimate of the actual hazards. Beyond introducing the FoS functions, we also identify all failure modes, some of which have not yet been investigated. We compute the FoS functions for the boulder Pena do Equilibrio, located in Spain. Our computations for FoS against sliding match all earlier results. However, we also compute FoS against toppling and against torsion and show that the latter may be critical. Since our methods are general, this suggests that torsion phenomena, which have been scarcely studied so far, may be relevant to analyze the stability of other natural rock boulders. A novel method to evaluate the safety of large boulders resting on steep slopes. Instead of Factors of Safety (FoS) computed at a fixed geometry with significant uncertainties, FoS functions are introduced, depending on a single geometric parameter. Considering spatial components of possible displacements, all failure modes of the boulder are identified, some of which have never been investigated before. Using geometric 3D analysis, our method is applied to the boulder Pena do Equilibrio, located in Spain and we find that torsional instability, one of the newly identified failure modes is critical. A novel method to evaluate the safety of large boulders resting on steep slopes. Instead of Factors of Safety (FoS) computed at a fixed geometry with significant uncertainties, FoS functions are introduced, depending on a single geometric parameter. Considering spatial components of possible displacements, all failure modes of the boulder are identified, some of which have never been investigated before. Using geometric 3D analysis, our method is applied to the boulder Pena do Equilibrio, located in Spain and we find that torsional instability, one of the newly identified failure modes is critical.


Geological Background and Problem Statement
Granite boulders are geomorphic structures, typically displayed in a wide variety of shapes, ranging from ellipsoidal, relatively slender blocks through cuboids to almost perfect rock spheres. Their size can also span from several centimeters in diameter to tens of meters. Many of these structures are formed as a result of a long, two-stage, insitu physico-chemical weathering process of granitic rock masses (Twidale 1982). Internal fracture networks represent preferential paths for water percolation, which causes relevant degradation of part of the rock mass, turning it into a set of rounded, sound rock blocks embedded in decomposed granite. If this material is eroded, the remaining rock mass emerges as a group of piled boulders. At this stage, and as a consequence of their rounded shape, these boulders tend to freestand on relatively small contact areas, sometimes close to precarious equilibrium conditions (Fig. 1a,c) (Domokos et al. 2012).
Other examples of these balanced rocks are represented by allochthonous blocks (also known as 'erratic blocks') transported by glaciers (Fig. 1b) or those derived from rockfall events (Fig. 1d). These shapes, resulting from fracture (Domokos et al. 2022) and abrasion in a glacial environment, typically display sharp edges, unlike the rounded shapes resulting from weathering.
The interest in studying the stability of precariously balanced rigid bodies arose in the 1960s after some earthquakes took place in Chile. At this time, some slender structures were found to remain surprisingly stable after the occurrence of the seismic events, whereas others, apparently more stable, were seriously damaged. These observations were noticed by Housner (1963), who first proposed a set of equations to model the rocking motion and damping of a rigid block (inverted pendulum). After this seminal approach, several authors have shown interest in studying the aforesaid problem, improving the existing solutions over a number of years (Yim et al. 1980;Mohammad et al. 1980;Ishiyama 1982;Makris and Roussos 2000;Prieto and Lourenço 2005).
In a parallel manner, from a rock engineering perspective, the first approach considering the study of the toppling mechanism in a 2D rectangular block and a rock slope was developed by Ashby (1971) and later improved by other authors (Hoek and Bray 1974;Sagaseta 1986). These results were based on static assumptions and did not consider a seismic trigger when dealing with the stability against toppling. Coombs et al. (1976) found in precariously balanced rocks (PBRs)-a real example of the problem of a rigid block susceptible to toppling-appropriate elements to be used as natural seismoscopes, helping in constraining estimates of peak ground motions. Later, Brune et al. (1996) developed an analytical-numerical methodology for estimating the ground acceleration necessary to topple such a PBR, relating this value with the existence of past seisms.
Some recent works started analyzing the mechanism of toppling in blocks with shapes closer to reality, i.e., Alejano et al. (2010) examined the sliding and toppling stability of a huge granitic boulder through a limit equilibrium approach and geometrical information gathered with LiDAR. The influence of rounded corners in toppling for both single and groups of blocks was studied by Alejano et al. (2015Alejano et al. ( , 2018 in the laboratory and on-site, respectively. Vann et al. (2019) developed a practical reliability-based approach to evaluate the stability of precariously balanced boulders, which accounted for the 3D effects of real natural contact surface on the stability. Almost at the same time, a limitequilibrium-based approach was proposed by Pérez-Rey et al. (2019) to evaluate the factor of safety for sliding and toppling mechanisms of a large granitic boulder.
The limitations of a 2D stability analysis approach for PBRs have been noted by various authors, leading to performing comparisons between 2 and 3D calculations (Haddad et al. 2012). Domokos et al. (2012) formally investigated, through 3D models of pebbles, how the static equilibria of rocks exist at two scales, explaining the phenomenon of rocking stones. Later, Domokos et al. (2022) determined the relation between the number of macro-and micro-equilibria and described their coevolution on a surface prone to abrasion.
The advance and affordability of 3D point-cloud acquisition and management lead some authors to use 3D point clouds as useful tools to analyze the stability of PBRs. Zábranová (2020) demonstrated through the registration of eigenfrequencies and numerical modeling based on a 3D point cloud that the use of a 2D approach for estimating the stability against toppling of PBRs might not be appropriate. In the same article, the authors also emphasized the relevance of local irregularities on the contact area when estimating the stability of these structures. Muñiz-Menéndez et al. (2020) presented a methodology to analyze the sliding and toppling stability of a large granitic boulder using 3D discrete numerical modeling. Following this line of thought, Pérez-Rey et al. (2021) presented a physical modeling approach for analyzing toppling stability based on laboratory tilt tests, limit equilibrium computations, and a 3D-printed version of a real PBR obtained from a 3D point cloud. In this article, the relevance of the 3D component in the stability of boulders was also highlighted. Recent work presented by Saifullah and Wittich (2022) also resorted to an LiDAR-based 3D point-cloud model of a PBR to analyze the toppling stability response by means of experimental and 3D discrete numerical modeling approaches. They also studied the influence of the contact normal stiffness on the probability of toppling.
The study of PBR stability has been the subject of research over several years, as it can be derived from the already referred works. These structures have represented natural tools for paleoseismic analyses, as they are reliable indicators of maximum historical seismic activity. Their stability assessment is of great relevance in the field of civil engineering due to the risk their failure may pose in infrastructures. Additionally, precariously balanced rocks are particularly interesting geomorphic structures, representing important natural heritage elements with potential touristic interest.
The safety of boulders at a counterintuitive balanced position is a delicate problem. In a similar way to a wide variety of civil engineering structures, the critical question is a quantitative safety measure expressing the resistance against toppling, sliding, and any possible movements of the block. In the case of artificial structures, uncertainty is mainly attributed to the spatial and temporal distribution of the actions (loads) and the material performance (strength). At the same time, the geometry is assumed to be known with a marginal error. This strategy is well reflected in engineering practice, where safety (or partial) factors apply to loads and material qualities but not to the geometry. In the case of boulders, such an approach is controversial: the geometry and, in particular, the exact contact surface entail a significant uncertainty. One possible approach to address this difficulty is applying complex probabilistic methods that treat the contact surface as a random function. However, such an approach relies on so many additional assumptions that any numerical outputs might be questioned from the point of view of practical engineering. In this paper, by resorting to a geometric 3D analysis of a large granitic PBR (Pena do Equilibrio boulder, located in Spain and previously shown in Fig. 1a), we present the concept where the geometric uncertainty about the contact surface is expressed by the single variable d, denoting the impression of the boulder, which is the parallel distance between the contact surface and the original rock surface (see Fig. 2). Although several Factors of Safety (FoS) may sensitively depend on the magnitude of d, this value is hard to measure or estimate with high accuracy. To accommodate this situation, we propose to use FoS functions, which are FoS computed as functions of d, instead of constant values. These functions offer a broad and balanced view of safety and appeal to the engineer's intuition without involving complicated probabilistic tools. They admit an immediate, intuitive view of safety, as they are plotted against the geometric parameter d, on which the geometry of the contact surface sensitively depends. In the case of failure modes with a strong dependence on the exact geometry of the contact surface (e.g., toppling), FoS functions admit a much more reliable estimate of the actual 1 3 hazards potentially induced by this structure. Needless to say, the computation of FoS functions relies on a fully 3D, point-cloud-based numerical analysis of the mechanics of the boulder's stability.
Beyond introducing the FoS functions and calculating them for this boulder, we also identify all failure modes, some of which had not yet been investigated. Our computations for FoS against sliding match all earlier results (Pérez-Rey et al. 2019). However, we also compute FoS against crushing, toppling, and against torsion/rotation and show that the latter may be the most critical mode of failure. Given that the methods presented in this paper are general, torsion, which has not been studied before, may also be critical in the case of other PBRs.

Main Hypotheses
Our approach is illustrated in Fig. 2, and it is based on the following assumptions: (a) the contact surface is planar; (b) its slope is known.
In this approach, the only unknown parameter is the impression depth d identifying the location of the contact surface with respect to the original rock surface. We use d as a single geometric parameter to characterize the variation of the safety factor, which we compute as a function of d. This computation relies on the following additional assumptions: (c) We assume that the contact area was created, while the boulder material was crushed. Crushing might happen under freeze-thaw cycles; beyond this contribution, the presence of water on the contact surface is not considered further, because it is presumably dry during most of the boulder's lifetime. (d) We assume that from the 3D dataset available on the existing rock surface, the missing (crushed) part can be reconstructed by extrapolation (see Fig. 2b).

The Considered Failure Modes
As noted earlier, the contact area, and hence the equilibrium is rather sensitive to the value of d. Our assumptions enable us to compute the contact area as a function of d (Fig. 2c). We identify four failure modes with corresponding FoS functions: sliding (FoS s ), toppling (FoS t ), crushing (FoS c ), and rotation around the contact normal (Fig. 3). In the absence of an apparent scale, the latter is given as the critical torsion T crit needed to dislodge the boulder. Safety is computed for two loading cases: pure gravity and gravity combined with a seismic action, where both an average and an extraordinary earthquake are considered. Our method is illustrated via the example of Pena do Equilibrio in Spain (Fig. 2). This boulder with a mass of 365.5 t appears to be on the verge of dislodging. The details of the mechanical assumptions and the computational method are summarized in the following sections.

Mechanical Considerations
We study the equilibrium of a homogeneous boulder B standing on a slope S (Fig. 4). The slope is assumed to be a plane with a dip β. Nonetheless, 0 ≤ β ≤ π/2 holds. We investigate the question under which conditions the currently  The boulder B on a slope (ramp) S with a dip. The boulder's weight W acts at the center of gravity G. Red arrows denote the admissible movement components that vanish in the case of equilibrium observed static balance position of the boulder B may be maintained under gravitational and seismic loads. We identify four physically relevant scenarios (failure modes) in which the boulder may get dislodged, namely: sliding, toppling, crushing, and rotation about the contact normal. These failure modes are illustrated in Fig. 3.
We aim to determine the factors of safety (FoS) for each of these cases, where FoS = 1 means that any incremental load could move the boulder.

Conditions of Equilibrium
For a horizontal slope (β = 0), equilibrium under the gravitational force is obtained if the boulder rests on one of its balancing points (Domokos et al. 2009(Domokos et al. , 2012. In geometric terms, the balancing or equilibrium points are those surface points of the convex hull of B at which the surface normal passes through the centroid G. In this case, the gravitational force is balanced by a vertical contact force, so equilibrium under self-weight is attained without friction. Under seismic action, either the dynamic response (i.e., the motion of the boulder) is analyzed, or in a simplified setting, a static load is used to characterize the effect of the seismic action. This latter approach is widely adopted in structural engineering, and we follow this interpretation in this paper (Chopra 1995).
If the dip angle differs from zero, friction becomes essential in maintaining the equilibrium even under pure gravity. We assume that the contact between the boulder and the slope follows Coulomb's law in dry friction with a coefficient of static friction μ s . Note that Coulomb's law relates the normal and the tangential forces independently of the contact area (Popov 2017).
To simplify the analysis, we assume that the boulder is locally convex in the contact's neighborhood. Strict convexity in the neighborhood of the contact is equivalent to point contact between the boulder and the slope. In mechanical terms, such an assumption is equal to an unlimited compressive strength of the material (both for the boulder and the slope). The unlimited compressive strength assumption can be justified for smaller objects, such as pebbles (Domokos et al. 2012), but the performance of the material should be included in a realistic model for boulders.

Role of the Material
Let f d denote the compressive strength of the boulder and assume that the compressive strength of the slope's material exceeds f d . Limiting the compressive strength eliminates the possibility of a point contact between the boulder and the slope, as it dictates a contact surface with a non-vanishing area. We assume that a polygon denoted by C encloses the contact surface with area A (Fig. 5). Nonetheless, for a contact transferring normal force N, the area is bounded below via Eq. 1 In our work, we assume a perfectly plastic contact. In other words, the boulder and the slope are assumed to be rigid bodies as long as the stress is below f d . As we observe (and record) the current, loaded geometry of the boulder, this assumption is justified.

Equilibirum Under Gravity
In the case of pure gravity, solely a vertical load W is acting at the G centroid of the boulder. Here, W = mg, where m is the mass of the boulder, and the gravitational acceleration g with a magnitude of g = 9.81 m/s 2 . The normal and tangential components of W with respect to the slope are denoted to W N and W T , respectively. Equilibrium is demonstrated by excluding all possibilities of rigid-body motions. Following Fig. 4, let u and ρ denote the vectors of displacements and rotations at point P, where P is the vertical projection of the centroid G onto the slope S. In the local basis [xyz] axis z and x in respect coincide with the normal and gradient of S. Based on our assumptions, W N is a compressive force, and we exclude cohesive stresses in the contact. Hence, W N should be balanced by compressive stresses denoted to (x, y) . The equilibrium conditions read the following: u z = 0, if σ ≤ f d in all points of the contact polygon C, -u y = 0, because the load W has no y-component, -u x = 0, if friction hinders sliding. As the frictional force should exceed W T , it follows that (2) tan ≤ s , Fig. 5 Under a seismic action with maximal value H and gravity W, the resultant of the loads pass through at a point inside ellipse E, as the direction of the horizontal H is arbitrary -ρ x = 0 and ρ y = 0, if W passes through the interior of the contact polygon C (i.e., P ∈ C) . Furthermore, there exists a normal stress distribution that balances W N and σ ≤ f d holds in all points of the contact polygon C, -ρ z = 0 holds if ρ x = 0 and ρ y = 0 also hold. (The proof is given in Appendix A.)

Equilibrium Under Seismic Actions
Solely, the horizontal component of the seismic action is considered; the associated horizontal force H is given via H = αW, where α is the seismic coefficient characteristic for the area. Note that the direction of the horizontal H is arbitrary (Fig. 5).
The conditions of equilibrium accordingly. ρ x = 0, and ρ y = 0, if the resultant force passes through the interior of the contact polygon C (i.e., in Fig. 5, the ellipse E is inside C). Furthermore, there exists a normal stress distribution that balances W N and σ ≤ f d holds in all points of the contact polygon C. Note that due to the dynamic nature of the earthquake, we do not require an admissible normal stress distribution for the resultant load W + H, only for W. This simplification is in agreement with experiences on natural and artificial structures in seismically active areas (Li et al. 2009;Alhajj Chehade et al. 2021). -ρ z = 0 holds if ρ x = 0 and ρ y = 0 also hold. (The proof is given in Appendix A.)

Failure Modes
Considering the conditions of equilibrium above, four failure scenarios can be identified. They are the following (and after each, we give the corresponding FoS): (a) sliding (FoS s ), (b) toppling (FoS t ), (c) crushing of the material (FoS c ), and (d) torsional moments acting onto the boulder.
(3) sin + cos cos − sin ≤ s , While in cases (a), (b), and (c), the Factor of Safety might be computed as the ratio of the resistance of the structure and the design action, in the last case, there is no apparent scale. Therefore, instead of a safety factor, we compute the critical torsional moment T crit necessary to dislodge the boulder.
As we have two load cases (pure gravity and gravity combined with a seismic action), in principle, eight safety factors might be associated with the boulder. However, crushing is assumed to be identical for the two load cases. In the lack of information on seismic torsion, T crit is only calculated for the pure gravity load case. Table 1 summarizes the computed safety measures.
Observe that a simplistic, planar model neglects failure mode (d) entirely and considers just an a-priori fixed orientation of the earthquake; hence, it might overestimate FoS t,e .

Safety Factor Functions
Pérez-Rey et al. (2019) determined the safety of the Pena do Equilibrio at a fixed contact surface obtained via pointcloud processing of LiDAR data. The reconstruction of the contact surface is not without ambiguities, and the computed safety factors are somewhat sensitive to the reconstruction of the contact polygon C. This latter depends primarily on the depth d of the impression, i.e., the distance between the plane of contact (i.e., the slope S) and the plane parallel to S and being tangent to the crushed (not anymore existing) part of the boulder. To illustrate the sensitivity mentioned above of the FoS s , we aim to compute the safety factor functions. These give the value of the safety factors w.r.t. the impression d.

Reconstruction of the Contact Surface
The original rock surface of the Pena do Equilibrio is impossible to measure below its contact plane, while data close to the ground are noisy; therefore, it was necessary to restore the entire rock surface by extrapolating from the clean part of the 3D LiDAR dataset. After meshing the recorded point cloud and choosing the orientation and dip of the slope S were fixed. Next, the 3D mesh was sliced at 1 mm intervals with planes parallel to the slope S. We found that the obtained cross-sections constitute a single polygon only if their distance from S was greater than or equal to 10 cm. It shows that the neighborhood of the contact on B is known inaccurately; hence, the part of the mesh closer to S than 10 cm was discarded. The arising hole in the mesh was patched by approximating the missing cross-sections with scaled copies of the hole's boundary, which is a closed, non-intersecting, planar polygon (Fig. 6a). To obtain the number of copies and the scaling factors, we fitted an oblate spheroid of radii r M and r m to approximate the boulder's surface. Here, r M > r m , and we constrained the oblate spheroid as follows: 1. Its center of gravity coincides with the boulder's center of gravity G. Some straightforward calculations yield r m = 1.961 m and r M = 4.161 m. This is consistent with previous results that measured a z G = 1.95 m distance between the boulder's center of gravity and its contact plane (Pérez-Rey et al. 2019).
The oblate spheroid protrudes to the opposite side of the contact plane S by 1.1 cm in the direction of its r m axis. In our model, we assume that the boulder also protrudes 1.1 cm below the contact plane and its tangent plane at the furthest point is parallel to the slope S. From now on, the impression d is measured from the tangent plane. Now, we need to rebuild the mesh around the hole. Let C 0 denote its boundary, which is a polygon in the plane at a d = 11.1 cm distance. Let A(z) denote the area enclosed by the circle formed as the intersection of the oblate spheroid and the d = z plane. We approximate the unknown crosssection of the rock surface with the d = z plane as C 0 scaled by a factor of √ A(z) A(11.1cm) . For example, the contact polygon C in the d = 1.1 cm distance is C 0 scaled by a factor of √ A(1.1cm) A(11.1cm) . We patched the hole with a new mesh by calculating the cross-sections using this technique with planes placed at a one-millimeter distance (Fig. 6b). Finally, faces between subsequent polygons were added (Fig. 6c, d).

Computation of the Safety Factor Functions
The safety against sliding is independent of the contact surface. Based on Eq. (3) Here, the factor FoS s associated with pure gravity is obtained via substituting = 0 into the formula above.
In the case of toppling failure, Pérez-Rey et al. (2019) computed the safety factor based on the critical dip of the slope needed to induce the toppling of the boulder. In Appendix B, we show that an equivalent formulation obtains FoS t by computing the support function h(.) of the convex hull of the contact polygon C (Fig. 7b). In specific, the general formula for FoS t,e reads where selects a point along the circumference of the contact polygon (see Appendix B for details). Here again, safety under static loads can be obtained by substituting = 0 . As the support function h(.) depends on the selection of the contact polygon, FoS t is a function of the impression d. Similarly, the safety against toppling in the case of seismic action depends on the impression d. The details of computing FoS t (4) FoS s,e = cos − sin sin + cos s. We also compute FoS c the safety against crushing using a f d = 100 MPa ultimate strength for the granite material of B and express this safety in terms of strength reserve. Here, an ideally plastic stress distribution in compression is computed over the contact polygon C for the normal force W N acting at point P (Fig. 7f). For the details, see Appendix A. If the plastic resistance N R exceeds the W N magnitude of W N , then the boulder is safe against crushing, and the safety factor reads Finally, the torsional moment around the axis normal to the slope and passing through the center of mass of the boulder needed to dislodge the rock is computed. In the case of point contact, the equilibrium is always unstable, because, in the lack of torsional resistance, an arbitrarily small torque can rotate B, which lowers the potential energy (i.e., the height of G) of the boulder. For an extended contact surface, friction stabilizes the equilibrium. We find that the safety against torque is linked to the safety against sliding, and we express safety in terms of a critical torque T crit (again, details are given in Appendix A).

Results and Discussion
We demonstrate the computation of the safety factor functions on the Pena do Equilibrio in Spain. The mass of the boulder is approximately m = 365.5 t (Pérez-Rey et al. 2019), and the dip of the slope is found to be β = 27 O . The coefficient of static friction reads μ s = 2/3. According to the Spanish National Act, in the region, the seismic coefficient α 1 = 0.032 on average and α 2 = 0.065 as indirectly estimated for exceptional earthquakes.
Regarding sliding, we find FoS s = 1.31 and FoS s,e1 = 1.21 and FoS s,e2 = 0.75, respectively. The first two agree with the values given in Pérez-Rey et al. (2019), and the last one, FoS s,e2 shows that in the case of an exceptional earthquake that acts close to parallel with the gradient of the slope, the safety of the boulder is lost.
In the case of toppling, the safety is slightly below the values of Pérez-Rey et al. (2019). In specific, at d = 12 mm, which provides a contact area identical to the area considered in Pérez-Rey et al. (2019), we find that the boulder might topple under an average earthquake around the axis denoted in red in Fig. 7c. The difference is attributed to the fact that in the new method, all possible toppling axes are considered.
The numerically computed safety factor functions for sliding, toppling, crushing, and the critical torsion function are summarized in Fig. 8

Conclusions
Spatial analysis of the precarious balance of granite boulders is essential in estimating their safety. A complete verification requires that the analysis includes rotation around the local normal of the contact beyond sliding and toppling. We find that a considerably small torsional moment is sufficient to drive the boulder out of its balanced position. Although pure gravity does not induce such a torsional moment, either Fig. 7 Steps of the computation for a contact polygon C at d = 14 mm impression. The polygon was obtained from the extrapolated surface (a) and its convex hull (b). The axis of toppling under gravity (green) and seismic action (red), respectively (c). The ellipses represent possible positions of point P for a generic (d) and an exceptional earthquake (e). The compressed part of the contact surface assumes a plastic stress state (blue polygon, f) artificial actions or even a compound seismic action might produce rotation around the axis normal to the slope. This observation underscores the importance of spatial analysis as pure planar models miss to grab the rotational phenomenon.
Considering the significant weight of the boulders, the crushing strength at the interface also affects the boulder's safety. In the example of the Pena do Equilibrio in Spain, we demonstrate that, as opposed to classical problems in structural engineering, the uncertainty associated with their geometry, especially the shape of the contact surface with the foundation of the boulder, should be accounted for in the analysis. By introducing safety factor functions, we expressed the uncertainty of the contact geometry in a oneparameter-dependent manner. We find that the safety of the Pena do Equilibrio cannot be verified under an extreme seismic event.
The interplay between rotation around the normal of the contact surface and seismic action calls for further, more detailed analysis in the future. Similarly, the case of nonplanar or partial contact between the boulder and its foundation also provides a promising future research direction, besides considering the existence of rock bridges that would ultimately improve the strength properties of the contact surface.
Appendix A See Fig. 9.   Fig. 8 The safety factor functions with reference to the impression d: a sliding, b toppling, c crushing, and d the critical torsional moment T crit. In panels (a) and (b), the results for pure gravity are given in purple; for the seismic action, yellow and brown correspond to an average and an extraordinary earthquake, respectively. The vertical green lines indicate the case of d = 12 mm impression. The horizontal blue lines correspond to FoS = 1 (where it is relevant)

About Rotational Equilibrium and Computation of T crit
Here, we demonstrate that the safety of the boulder on a slope against toppling implies safety against torsion around axis z and compute the T crit critical moment that dislodges the body. First, assume the slope dip is critical, i.e., tan β = μ s holds. If there exists a distribution of the σ(x,y) compressive normal stresses over the contact polygon C, such that at each point of the contact surface, the boulder is in equilibrium. Observe that the direction of the shear stress at each point of C is parallel to the axis x, i.e., the slope gradient. Equation (A4) yields that the T torsional moment vanishes, because The obtained result also shows that the shear capacity of the contact is entirely exploited to prohibit sliding. Therefore, the shear resistance of the surface cannot contribute to hindering torsion. It means that even a slight torsion can dislodge the boulder.
If the slope dip is not critical, i.e., tan β < μ s holds, then the shear capacity of the contact might contribute to hindering torsion. The computation of the T crit torsional resistance is a delicate question. Here, we aim to approximate a lower bound on T crit . First, a normal stress distribution with a minimal area is sought. Based on classical considerations (Timoshenko 1956), we assume that the contact is without cohesive forces and the compressive σ(x,y) normal stresses are in a plastic state. Let C * ⊆ C denote the plastic zone under the compressive force acting at P (Fig. 9(a)). Assuming a line separating the compressed and non-stressed regions, the solution of the equilibrium equations is unique (or there exists no solution). Let A* denote the area of C*. To find the bound on T crit , we assume that the points of C* with a high distance from P contribute entirely to resisting sliding, and the shear stress over the circular area in the vicinity of P resists torsion. This choice is motivated by minimizing the arm of stress in the torsional moment. It means that C* is separated into two disjoint sets: the shear stresses in C S * hinder sliding, and the shear stresses in C T * hinder torsion (Fig. 9 b).
The compressive resistance of the contact reads which value is used to compute the FoS c factor in Eq. (5). It is straightforward to show that the lower bound on T crit requires a minimal A T * area of C T *. This minimum is attained if the normal stress acting on C T * is maximal, i.e., it equals the compressive strength f d . From the equilibrium equations, these considerations yield Let R T denote the radius of the circle with area A T *. The torsion resisted by the shear stresses is found after the integration is performed

Appendix B
See Fig. 10.
(15) T crit = 2 3 s f d R 3 T . Fig. 9 Computing the lower bound on the critical torsion T crit . The compressed zone for the contact C in the plastic state (a). The shear stresses associated with sliding (blue) and torsion (red (b)) Funding Open access funding provided by Budapest University of Technology and Economics. The work of Gábor Domokos was supported by Nemzeti Kutatási, Fejlesztési és Innovaciós Alap, K134199, and the National Research, Development, and Innovation Fund of Hungary, TKP2021-BME-NVA-02.

Data availability
The computational data used in this study will be made available on request.

Conflict of Interest
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.
(18) FoS t,e = min h( ) cos + sin 1 z G tan . Fig. 10 Analysis of the contact against toppling failure under pure gravity (a) and gravity combined with seismic action (b)