In order to take advantage of a refined 3D view of the study area, we have performed a photographic coverage of the Mt. Peron niche area through a drone flight.

We used a DJI Phantom 4 Pro Plus multi-rotor UAV (unmanned aerial vehicle) equipped with a built-in high-resolution camera. The built-in lens has a field of view of 84° with a focal length of 8.8 mm/24 mm (35 mm equivalent focal length). The camera is equipped with a 1" 20 MP CMOS sensor, and we recorded photographs in RAW format; all the photographs were manually shot with the drone hovering.

We chose to model the three major pillars visible in the niche area (Figure 3), considering that these are the most hazardous features for potential instability. It is noteworthy that the niche area presents an elevation drop of ca. 500 m with an almost vertical rock wall. It is an inaccessible site and remote sensing is the most important data source for this case study. A total of 145 photographs (Figure 4A) were shot following a parallel-axis shooting strategy with the addition of a grid of oblique photographs, especially for the upper sector, where complex rock dihedra are present. Photographs were processed in RAW format and were acquired from 50 to a maximum of 200 m from the target surface, with most of the images shot at ca. 100 ± 50 m. We obtained an optimal overlapping of the photographs for the 3D model building, with the modeled area covered by at least nine photographs for each sector of the point cloud (Figure 4B).

We used the SfM photogrammetry software Metashape (Agisoft, 2022) to produce a 3D point cloud from the UAV photography. The camera alignment was obtained; thanks to the co-registration of 81,306 tie points, and the internal accuracy of the model was tested over 16 control points manually checked over the different photos, with an RMS reprojection error of 0.189 m. The dense point cloud was built without down-sampling of the original images (i.e., ultra-high-quality setting in Metashape) and with a mild depth filtering.

We were not able to position ground control points (GCPs) on the site due to the inaccessibility of the surveyed cliff. We then relied on the direct georeferencing of the model by means of the photographs’ geotagging, as recorded by the built-in navigation—grade GNSS receiver. We are aware that this is a source of uncertainty in the model reconstruction; so, after the point cloud generation, we analyzed the accuracy of the model positioning through an ex-post 2D validation approach. We provided an estimate of the accuracy in the absolute georeferencing of the virtual outcrop by comparing the orthophoto derived from our virtual outcrop with a high-resolution orthophoto-mosaic derived from an aerial coverage of the area from the Veneto Region (AGEA 2015; last accessed in 2022; courtesy of Veneto Region). The reference orthophoto coverage has a nominal resolution (GSD) of ca. 30 cm (specifically 50 cm in mountain regions and 20 cm in plain sectors).

For fracture digitizing and orientation measurements, we used the CloudCompare software (CloudCompare, 2022). CloudCompare is a non-commercial processing software for 3D point clouds and meshes and is available under the GNU General Public License. It includes many functionalities and is widely used for its versatility and the availability of many science-oriented plugins developed by the user community. We performed data digitization through the plug-in Compass (Thiele et al., 2017) which provides an assisted interpretation procedure for linear fracture detection. Planar surfaces can be directly measured on the point cloud, if well exposed, by clicking on the fracture surface; the calculated orientation is the best-fitting planes passing through the points included in the circular window around the clicked point. Alternatively, the estimation of the structure orientation is based on its intersection with a non-flat surface (i.e., the virtual outcrop surface). Linear traces of discontinuities on the virtual outcrop can be digitized in CloudCompare by manually selecting the points including the discontinuity trace. Then, a computer-assisted procedure helps in selecting all the points in the cloud connecting two near points. The software uses a least-cost-path algorithm to follow these discontinuity traces along the surface between the user-defined start and end points and then calculates the best-fitting plane to estimate the orientation. The least-cost-path procedure, applied on points included in a spherical search, is based on user-defined scalar values, associated with each point in the cloud (e.g., color similarity, brightness, or presence of high-curvature regions—edges). We, here, adopted an approach based on the brightness of the points (i.e., the absolute value of the associated color in a greyscale rendering), since fractures are visually associated with hairline shadows in the acquired images. Finally, once the fracture trace on the modeled surface is defined (Figure 5), the best fitting plane is calculated (i.e., the plane where most of the points on the fracture lie) together with its orientation (i.e., dip direction and dip).

Alternatively, if most of the detected points of the fracture trace are almost collinear and the best-fitting plane is not adequately constrained, the fracture can be manually digitized by selecting only three points of the fracture and proceeding with the classical solution of the three-point problem.

Finally, poles to fracture planes are statistically analyzed in order to calculate the average orientation of each fracture. For this statistical analysis, we used the software Wintensor v.5.8.8 (Delvaux and Sperner, 2003). First, we clustered the poles to the fracture planes through a K-means algorithm (MacQueen, 1967): the number of clusters is fixed a priori by the operator, and the calculator finds the centroids of the clusters that minimize the total intracluster variance and maximize the intercluster distance. A number of clusters were fixed with a visual inspection of the data, aided by density contouring of fracture poles. The output includes the orientation of the fracture sets, the Fisher’s dispersion factor (i.e., K value, Fisher, 1953), and the mean cone angle (MCA) enclosing the data distribution.

To evaluate the stability of the three rock pillars under investigation, the distinct elements method (DEM) was used (Cundall and Strack, 1979). This method is widely used to simulate static and dynamic problems involving dry or wet granular media of any size and shape (Gabrieli et al., 2013), from powders (Martin et al., 2003) to sands (Gabrieli et al., 2009), gravels, and fractured rock masses (Spreafico et al., 2016; Calvetti et al., 2019). Each body is treated as a distinct rigid element having six degrees of freedom, and it communicates with the other bodies through specific contact laws. The algorithm is explicit and solves the equations of motion of all bodies in the time domain for both translation and rotation through a finite difference integration scheme.

At each time step: i) the list of contacts between the various blocks is updated, ii) each contact provides a normal and tangential force that is transferred to the center of mass of each block, iii) the body forces and possibly other external forces are added to the resulting force, and iv) the positions and velocities of the blocks are updated, integrating over time.

In this case, the method is particularly suitable for the intrinsically fractured nature of the rock mass along predefined planes and for the ability to manage multiple contacts, fractures, and detachment problems between the blocks, even with large deformations.

The complexity of the rock mass and the three prisms and, above all, the presence of the joint planes (that will be called K1), which are oblique with respect to the other three (K2, K3, and S0), made it necessary to use a 3D numerical modeling. Moreover, the presence of several sub-vertical bedding planes (called S0) means that these three pillars cannot be assimilated to three non-fractured blocks that detach but rather to a cluster of adjacent blocks that could potentially detach, slide, and overturn.

The software used for this analysis is 3DEC by Itasca (v. 5.2, 2016; Itasca, 2016), which can generate polyhedral 3D blocks starting from the deterministic and statistical description of the planes that delimit them (Discrete Fracture Network—DFN). It also internally implements a method for evaluating the safety coefficient of the rock mass, similar to the strength reduction method used in finite element codes.

By means of the SfM-VsM workflow, we obtained a sparse point cloud from camera alignment with a mean key point size of 2.2 pixels, a maximum reprojection error of 0.56 pixels, and a mean tie point multiplicity of 8.5. The estimated relative accuracy in the camera position has been evaluated through an application of the Metashape software, and we found an average accuracy of 1.5 m (Figure 4C). The dense point cloud, calculated using the MvS processing, consists of 149,814,051 points, after trimming, with an average spacing of ca. 5 cm (Figure 4D).

We performed a 2D georeferencing to check if the positioning errors were mainly due to shifting and/or tilting or an internal deformation of the model. This procedure is a simplification of a real 3D georeferencing of the model but, since accurate GCPs are not available, it still constitutes an exploration tool that can provide the relative distribution and orientation of errors projected on a horizontal plane. The accuracy of a zero-order shift (i.e., a simple shift on the horizontal plane) was 5.65 m and of a second-order transformation (i.e., implying internal deformation of the model) of 3.18 m. Conversely, assuming a first polynomial transformation (i.e., implying shift scale and rotations only), we obtained a root mean square error of 1.47 m (a maximum error of 4.83 m), with greater errors aligned along the y-axis. We can thus assume that the 3D georeferencing of the virtual outcrop can suffer from small misorientations, especially around the x-axis of rotation, but also those internal distortions are not the main source of error. If we assume that internal distortions are relatively small and consider the maximum error (4.83 m) as due to the model misorientation, this will result in ca. 1–2 degrees of orientation errors for fractures as large as those digitized in the outcrop (i.e., tens of meters long).

Considering that our scope is to perform a first-order stability analysis of a rock cliff where no other detailed measurements are available, we, here, assume that the estimated error is small enough to be accepted for the successive analysis.

On the virtual outcrop, we measured 159 discontinuity planes that have been clustered into four sets of fractures (i.e., S0, K1, K2, and K3; Figure 6; Table 1) with 18 data resulting as non-classified. The first cluster, S0, represents the bedding, corresponding to sub-vertical planes striking ca. WSW-ENE and located in the forelimb of the anticline, in the hanging wall of the BL thrust. The cluster K1 represents a principal sub-vertical fracture system, NW-SE-striking, that cuts through the bedding and, together with S0, isolates the main pillar and rock prisms on the analyzed wall (Figure 6C). Clusters K2 (NE-SW striking high-angle fractures) and K3 (low-angle to sub-horizontal planes) represent two secondary and less pervasive discontinuity systems.

The four detected fracture sets have been compared with field-based measurements, taken at the base of the rock cliff and along a rappel descent on the same sector of the three pillars and during a field survey in areas close to the sides of Mt. Peron (Figure 7).

Field-measured fracture sets were coded as S0 (for bedding) and as J1 to J6 classes. Our clusters S0 and K1 from the UAV survey show good agreement with S0 and J1, respectively (Figure 7).

The point cloud of the rock mass obtained through the SfM-VsM workflow was used to describe in a simplified way the external surface of the rock mass in the DEM model (Figure 8A). It should be noted how a higher detail in the description of the model surface has no benefit in this analysis: the three prisms are quite big (the smallest is 100 × 15 × 20 m) and almost regular and higher detail of the surface would be functional to the stability analysis of detachment of the smallest blocks (metric and sub-metric) but not to an overall stability analysis; moreover, it would greatly increase the computational cost of the model. On the contrary, more detail is needed in the definition of the joints because they affect the level of fracturing of the rock mass and have a great impact on the stability analysis.

First, the planes delimiting the three prisms from the rest of the rock mass have been added (namely S0 on the backplane and K1 and K2 for the other sides). The rock mass behind and at the base is assumed to be stable and fixed in the model. Subsequently, the measurements on the joints coming from the virtual and conventional discontinuity mapping have been used to create a DFN. For this purpose, we used WinTensor and the DFN module of 3DEC by Itasca (v.5.2, 2016; Itasca, 2016). We adopted a probabilistic joint orientation according to a Fisher distribution (Fisher, 1953) with a dispersion factor k calculated on the outcrop data (Table 1). For the fracture size distribution, we assumed circular fractures of finite size following a power law with a scaling exponent equal to 3.2 (Priest and Hudson, 1981; Priest, 1993; Bonnet et al., 2001). The density of the joints was expressed as p10, in terms of the number of fractures per unit of length; alpha was estimated along three scanlines in the point cloud. In 3DEC, fractures of the DFN have been generated until a target value of p10 is reached for each joint set. Although as evidenced by other authors (Elmo and Stead, 2010; Dershowitz et al., 2014), other parameters such as the areal intensity (p21) or the volumetric intensity (p32) can provide a more realistic description of the fracture density of the rock mass, considering the fairly regular shape of the pillars with some faces parallel to the joint planes, we decided to use the p10 value as the DFN generation parameter, taking care to choose scanlines orthogonal to these faces and thus maximizing its information content.

The final DFN model with the four discontinuity sets S0, K1, K2, and K3 used in the simulations is depicted in Figure 8C while their orientations with the poles of the external surfaces of the three pillars superimposed are represented in Figure 8B. Note how the three pillars share the same orientation of the south-facing face, almost parallel to the bedding plane S0, while the other subvertical surfaces confining the three dihedrals belong to the K1 plane subset. This particular fairly regular configuration of the three vertical prisms being analyzed, with these two well-exposed surfaces forming an angle between 70° and 55°, has allowed a sufficiently detailed statistical reconstruction of the families of the internal joints.

For the constitutive model of the joints, a generalization of the Mohr–Coulomb shear strength envelope was also assumed with a normal tensile limit value while rock blocks were simulated as non-deformable bodies having a unit weight of 25.5 kN/m³. The mechanical parameters used in the model were estimated based on the field measurements obtained by Di Giusto (2012) and from samples taken at the toe of the vertical slope. The same parameters were, therefore, extended to the other three pillars under investigation. The joint compressive strength was estimated with the Schmidt hammer (Deere Miller, 1966; Aydin, 2008) while the RQD index was determined from the density of joints (Palmström, 1982). Spacing conditions of discontinuities and groundwater conditions have been estimated quantitatively and qualitatively by geol. Di Giusto during the field survey, providing a final RMR estimate of 58 (Barton 1974; Bieniawski, 1989). On the basis of additional sclerometer tests and roughness measurements in the niche at the toe of the first pillar on weathered and unweathered rocks, the shear strength of joints was evaluated according to Barton and Choubey (1977). The basic friction angle was assumed to be equal to 32° for this rock type. The cohesion value attributed to the joint planes usually takes into account several aspects including the real cohesion of the intact rock and the fraction of the rock bridges on these surfaces (Jennings, 1970; Calvetti et al., 2019). In the absence of this information, we decided to quantify the joint cohesion by carrying out a small stability analysis on the large rock slab of known volume, placed at the top of the biggest pillar and adhered to the sub-vertical wall, with a procedure similar to that reported by Sampaleanu (2017). In fact, for this high-slope configuration, the friction angle of the joint plays a marginal role compared to cohesion and allows a better evaluation of the latter. It has also been hypothesized that this value maybe the same in the joints of the other pillars. This assumption represents a precautionary estimate of the real cohesion value since it is likely that the bedding plane S0 on which the stability of the slab was evaluated represents the plane with the greatest degree of alteration and lowest normal stress and therefore the lowest cohesion value. Due to the difficulties in getting samples and the lack of other specific geomechanical tests, the joint stiffnesses were not calculated, but realistic values have been assumed according to literature (Fairhurst, 2014). We used a fairly low value to compensate for the rigid block assumption and the reduction of joint stiffness for tension mechanisms which is expected for subvertical joints. In any case, it must be recalled that the stiffness of the joints has little or no influence on the evaluation of the safety factor which is instead greatly influenced by the strength parameters (Wines, 2020).

The values reported in Table 2 are therefore assumed for the joints.

For the evaluation of the stability of the rock mass, the strength reduction method, which consists in looking for the shear resistance reduction factor of the joints (friction angle and tangential cohesion only) that leads to the limit equilibrium condition (Matsui and San, 1992), was adopted. The calculation started by applying gravity with a fictitious high joint shear strength and reaching a stable condition. After that, the joint properties were updated, and a bracketing approach was used in the factor of safety calculation by narrowing step by step the factor span between a stable and unstable solution (Dawson et al., 1999).

The stability analysis was performed on each of the three pillars separately, i.e., by fixing two pillars each time and solving the equilibrium of the remaining one.

The factor of safety (FOS) and its associated instability mechanism were assessed for each pillar to identify those in more precarious conditions and not to evaluate a safety factor in an absolute sense.

It is emphasized that the simplifications introduced in the model (e.g., the assumption of rigid blocks, the discontinuity planes imposed with the DFN, and the parameters calculated as lower limits) do not allow a refined numerical analysis but to determine on which pillars to focus attention for further surveys. It should be noted that, for a generic FOS greater than 1, the corresponding instability mechanism is obtained by dividing the tangent of the friction angle of the joints and the shear cohesion by the same FOS and observing the kinematics of the blocks in these conditions. The displacements obtained with this method are, therefore, not the real displacements of the mass but highlight the type of rupture mechanism. The smallest pillar (V1) showed a stable condition with a FOS = 1.6 associated with a planar sliding instability mechanism and a tendency to overturn toward the south-east partially sliding on the bedding plane (Figure 9A).

The intermediate pillar V2 has a safety factor of 1.7, slightly higher than the first. This result can be explained by the fact that despite having a height similar to the smallest one but a width approximately double the bedding plane, it has the surface facing south-west and orthogonal to the bedding with a lower dip value (84° instead of 90°), and therefore, overall, a center of mass that makes the pillar more stable. The failure mechanism, in this case, is of the sliding type on the bedding plane with partial rotation at the toe (Figure 9B).

Finally, the last pillar V3 has a safety factor of just above unity (FOS = 1.02). Its southeast face has a slope angle on average higher than the other two prisms (78° instead of 73°) and more than the bedding dip. This is also evidenced by the 30 × 40 × 3 m prismatic rock slab, slightly bent at the base, that appears to be suspended on its top (Figure 10) and that was used for the calibration of the joint cohesion. Even removing these slabs having a considerable thickness and whose stability is mostly controlled by the joint cohesion, the safety factor remains low, and the mechanism is a complex failure controlled by the S0 bedding surfaces (Figure 9C). This can be explained by the particular geometry of the base of the prism which has an external surface with an inclination close to the vertical (88°) (Figure 10). The precarious conditions of this huge pillar are also confirmed by the evident detachment of blocks in recent times, as highlighted by the weathered rock surfaces at the foot.