Applying SLAM-Based LiDAR and UAS Technologies to Evaluate the Rock Slope Stability of the Grotta Paglicci Paleolithic Site (Italy)

: This study focuses on slope stability and geological hazard analyses at the Italian Paleolithic site of Grotta Paglicci. The site is characterized by a cave that contains rich archaeological and anthropological finds, spanning various Paleolithic periods, and includes faunal remains, lithic artifacts, human burials, ornaments, mobiliary art objects, and unique Paleolithic wall paintings. The study employs a multi-technique approach that includes topographic surveys carried out by the robotic total station and GNSS receivers, photogrammetric acquisitions with an unmanned aerial system, 3D SLAM-based LiDAR mapping, and an engineering geological survey. The collected data allowed for the creation of georeferenced 3D models that were utilized in rock slope stability analysis and modeling. The results of this comprehensive survey highlighted how the bedding and joint discontinuities influence rock stability in both the external and internal areas of the cave. The integrated use of SLAM-based LiDAR and photogrammetry has been proven to be an efficient and essential tool in the evaluation of the structural interactions between the external morphology and the cave, thus allowing the proposal of safety measures that will keep the site accessible for future activities.


Introduction
Natural caves, which are characterized by significant anthropological and geological heritage, are often affected by rock instabilities that compromise their access and human activities in safe conditions.Overcoming these limitations may require the employment of modern survey techniques that allow the assessment and mapping of geohazards.
Underground natural environments have typically complex structures that require investigation with the use of both multiple techniques and technologies.The most commonly utilized methods of spatial data acquisition in complex and extensive caverns are referred to as LiDAR (light detection and ranging) and digital photogrammetry.The application of solely ground-based LiDAR for spatial mapping is documented in several case studies, such as in the Las Caldas and Penade Candamo caves, in Spain [1], the Dachstein cave, in Austria [2], the Wonderwerk cave, in South Africa [3], the Parpello cave, in Spain [4], the Les Fraux cave, in France [5], the Pollera cave, in Italy [6], the Jenolan and Koonalda caves, in Australia [7], the Domica cave, in Slovakia [8], the Elvandi cave, in Ecuador [9], the Perhersk Lavra cave, in Ukraine [10], and the Grotta Giusti cave, in Italy [11].Nevertheless, a major drawback of ground-based LiDAR applications is the prolonged time required for data acquisition, as described in [12] for the ˙Ince giz caves in Turkey.
Other examples of terrestrial laser scanning being applied to hypogeal cavities with an additional geological aim include the Santa Barbara cave, in Italy [13], the Orgnac's cave, in France [14,15], the Märchenhöhle cave, in Austria [16], the Pastora cave, in Spain [17], the caves of the Ardèche-Chartreuse-Vercors-Bauges Département regional parks in France [18], the Algar do Penico cave, in Portugal [19], and the Grotta A cave, in Italy [20].
Most of the latter cited works employ traditional methods of rock slope stability analysis, such as the kinematic and dynamic methods based on rock mass characterization.Few papers include finer analyses such as numerical modeling [22,24].
Nowadays, LiDAR systems, such as the SLAM-based system (simultaneous localization and mapping), may allow greater coverage and speed, greater availability of data in remote or inaccessible areas, versatility regarding terrain, efficient data collection, and rapid data processing [26][27][28][29][30].This acquisition approach overcomes the major drawbacks typical of ground-based terrestrial laser scanning, which are due to longer acquisition times and, in complex and extensive morphologies, the presence of several station points, which are required to avoid shadow areas.
Moreover, in terms of spatial data acquisition, sometimes an additional survey, which includes not only the cave but also the external surrounding areas, can be necessary.This can be performed with similar ground-based technologies or UAS-based (unmanned aerial system) equipment.
In view of these premises, a multi-technique approach was adopted in this work to study the whole geological hazard of Grotta Paglicci, including both the external and the inner hypogeal system.
Grotta Paglicci is a Paleolithic cave located in the municipality of Rignano Garganico (Foggia, Italy; Figure 1).Ongoing archaeological and paleoanthropological research at this site has been conducted for over fifty years by the University of Siena, in collaboration with the local heritage office (Soprintendenza Archeologica, Belle Arti e Paesaggio of the Provinces of Barletta-Andria-Trani-Foggia).Moreover, Grotta Paglicci has never been studied from an engineering geological perspective, and stability analysis has not been performed until now.The actual cave consists of four main chambers (Figure 2).In particular, the Atrium contains an important 12-meter-thick anthropogenic deposit that represents one of the most complete Upper Paleolithic stratigraphic sequences in Europe.Notable findings include lithic artifacts, faunal remains, Gravettian burials, isolated human remains, ornamental objects, and portable art.Furthermore, wall paintings, presumably attributable to the Gravettian period, were also discovered in the innermost chamber of the cave (Chamber 3 in Figure 2).The adopted multi-technique approach consisted of: (i) a topo-cartographic survey, carried out using RTS (robotic total station) and GNSS (global navigation satellite system) receivers, aimed at data georeferencing; (ii) a spatial data acquisition executed by UASphotogrammetry for the external zones, and by SLAM-based laser scanning for the cave system.
Spatial data was integrated with information on rock discontinuities collected during the engineering-geological survey, aimed at characterizing all the joint systems and the rock mass.Moreover, a kinematic stability analysis was carried out for the slope.Finally,   The adopted multi-technique approach consisted of: (i) a topo-cartographic survey, carried out using RTS (robotic total station) and GNSS (global navigation satellite system) receivers, aimed at data georeferencing; (ii) a spatial data acquisition executed by UASphotogrammetry for the external zones, and by SLAM-based laser scanning for the cave system.
Spatial data was integrated with information on rock discontinuities collected during the engineering-geological survey, aimed at characterizing all the joint systems and the rock mass.Moreover, a kinematic stability analysis was carried out for the slope.Finally, The adopted multi-technique approach consisted of: (i) a topo-cartographic survey, carried out using RTS (robotic total station) and GNSS (global navigation satellite system) receivers, aimed at data georeferencing; (ii) a spatial data acquisition executed by UAS-photogrammetry for the external zones, and by SLAM-based laser scanning for the cave system.
Spatial data was integrated with information on rock discontinuities collected during the engineering-geological survey, aimed at characterizing all the joint systems and the rock mass.Moreover, a kinematic stability analysis was carried out for the slope.Finally, 2D finite element method (FEM) numerical modeling enabled us to define the distribution of principal stresses and displacements along several cross-sections.
Since the area, given its anthropological significance, is intended to be opened for future scientific campaigns, the obtained results, due to their high detail and comprehensiveness, are being used to plan future remediation efforts for site safety.In fact, the results concerning kinematic mechanisms and slope stability are fundamental for the main purpose of preventing and mitigating possible hazardous events.These actions will not only secure the safety conditions for personnel involved in future research but also preserve the site in the long term.
The remediation project, which will comprise the installation of support and reinforcement and the monitoring system, is not thoroughly discussed in the present manuscript, but how our findings led to the actual project's technical proposal will be explained in the discussion in Section 6.

Geological Setting
The study area is located at the border between two main geodynamic-structural domains (Figure 1A): the eastern Apulian Foreland, represented by the Gargano Promontory, and the western Bradanic Foredeep basin, known as the Tavoliere delle Puglie [31].
Limestones are the dominant lithology in the area, and their rocky layers exhibit thickness variations ranging from 30 to 200 cm and dip toward the south, varying from horizontal to 50-60 degrees [31].This lithology belongs to the "Monte Calvo del Gargano" member (CBA 1 in Figure 1B), which dates to the Callovian-Valanginian stages (165.3± 1.1 Ma −132.6 Ma) [32,33].This member forms part of the "Calcari di Bari" succession, which is Upper Jurassic-Lower Cretaceous in age (201.4 ± 0.2−66.0Ma).
The Grotta Paglicci study area has a complex structural history that is characterized by the presence of ancient strike-slip faults such as the Rignano Garganico Fault and the Candelaro Fault (Figure 1B).These faults have experienced reactivation with extensional or trans-extensional kinematics and are presently partially buried by Quaternary deposits, contributing to the formation of steep escarpments along the southern and southwestern edges of the Gargano Promontory [31].
The predominant morphogenetic process in the Gargano Promontory is karstification, which varies in distribution and shape depending on lithology, structural characteristics, and tectonic evolution [34,35].Within the study area, there are two primary karst hypogeal cavities: Grotta Paglicci and the easternmost, and least documented, Grotta dei Pilastri, which are interconnected by a relatively small cavity or tunnel.Above this connecting tunnel, a shallow sinkhole is evident on the surface, resulting from a possible collapse of the tunnel roof.Presently, this connecting cavity remains inaccessible and is covered by collapse debris, which includes blocks of varying sizes and fine material.Furthermore, within Chamber 2 of Grotta Paglicci (Figure 2), observations revealed that joints have allowed limited water inflows and enabled the root systems of the surface vegetation to reach the underground space.

Anthropological and Archaeological Setting
The importance of Grotta Paglicci as an archaeological site was detected for the first time in 1955 by Professor Raffaello Battaglia (University of Padova); preliminary fieldworks, directed by the Natural History Museum of Verona, in 1961 and 1963 brought to light a Paleolithic deposit.From the seventies onward, Grotta Paglicci was studied by the University of Siena, always in collaboration with the local heritage office [36].The site is composed of the present-day cave and a rock shelter that was part of the same underground system before the collapse of the ceiling (Figure 2).The importance of the site lies in the long stratigraphic sequence found both in the cave and in the rock shelter; indeed, the deposit covers a time span that embraces middle Paleolithic phases and an Upper Paleolithic sequence.It represents one of the most complete deposits in Europe and it can be considered as a reference for the Mediterranean area [37].A 12-meter-thick stratigraphy was found in the Atrium (the yellow star in Figure 2).Its bottom, referable to an ancient phase of the middle Paleolithic, partially correlates with the deposit of the rock shelter (the blue star in Figure 2) [38].The Upper Paleolithic sequence spans from the Marginally Backed Bladelets Aurignacian (about 40,000 years before the present (BP)) to the Final Epigravettian, dated to about 13,000 years BP.In the aforementioned chronological time span, the sequence covers all the Upper Paleolithic phases known in southern Italy.Archaeological layers yielded a large number of faunal remains, lithic artifacts, human remains (among which two complete Gravettian burials), portable art objects, and a portion of a wall painting that had collapsed from the vault [39][40][41][42][43][44][45][46][47][48][49][50].In addition, one of the earliest pieces of evidence of flour production [51] and the only evidence of Paleolithic dogs found in Italy to date (the yellow star in Figure 2) [52] were discovered.Some layers of this sequence can be correlated with a thinner deposit located just outside the present-day cave entrance (the purple star in Figure 2) [53] and with another stratigraphy found in the inner part of the hypogean system between Chambers 1 and 3 (the green star in Figure 2) [54].The layers are generally well-anthropized, except in some contexts, both in Atrium and Chamber 1, where evidence of hyena dens was detected (the yellow and green stars in Figure 2) [53][54][55].Moreover, some Upper Paleolithic parietal paintings (the only example discovered in Italy to date) are still preserved in the inner part of Chamber 3 (the red stars in Figure 2) [56,57].

Materials and Methods
To fully characterize the morphologies of both the cave and the external surrounding areas, an integrated approach including topographic, UAS photogrammetric, and SLAMbased LiDAR surveys was adopted.Moreover, for the analysis of the rock slope's stability, a kinematic stability study and 2D FEM numerical modeling were carried out, and analysis was also based on results from the in situ engineering-geological survey.

Topographic Survey
The topographic survey, conducted utilizing a Leica TM GS15 dual-frequency geodetic GNSS (global navigation satellite system) receiver and a Leica TM Nova MS50 robotic total station, was carried out to georeference the data and make them measurable.To reach this goal, the coordinates of various types of targets were measured, in particular (i) black/white targets for RTS scanning, (ii) yellow/black targets for the UAS photogrammetric survey, and (iii) circular-shaped optical targets for SLAM-based LiDAR scanning (Figure 3).
The coordinates of the targets located on the ground were measured using the GNSS receiver operating in NRTK (network real-time kinematic) mode.Measurements were executed, with acquisition times ranging from 20 s to a few minutes, to achieve positional accuracy values of about ±1 cm for planimetry and about ±1.2 cm for altimetry.The correction service used in NRTK is the HxGN TM Smartnet, while Leica TM Infinity software (Version 3.4.2) [58] was used to visualize and validate, in the laboratory, the measured points.Using the Convergo application and related geodetic grids, the points' elevation was transferred from the WGS84 ellipsoid to sea level.The GNSS data also allowed us to georeference the RTS survey, including points placed on the rock walls and in inaccessible areas.The relative accuracy of points measured by the RTS, considering that prisms were not utilized, can be estimated within the range of 5-8 mm.
The survey led to a total set of 30 measured targets that were then used to georeference the outputs obtained from the other techniques in the chosen reference system (i.e., ETRF2000/UTM33N).

UAS Photogrammetric Survey
The photogrammetric survey was performed using a DJI TM Mavic2 Pro drone equipped with a Hasselblad L1D-20C, with a 1-inch sensor and 20-megapixel resolution, combined with a 10.26 mm lens.The study area was surveyed by a total of 3 different flights: the first one was executed in automatic mode, with nadiral shooting, using UgCS software (Version 4.18) [59] for flight planning.Two additional flights were carried out in manual mode, with the direction of photo acquisition almost orthogonal to the rock slopes, and the camera tilting from 45° to 65°.In this way, it was possible to obtain complete coverage of the external part of the study area from optimal points of view.
A total of 550 images were acquired in an area of 1.9 ha, with a ground spatial resolution of about 1 cm.The flights were planned to ensure sufficient overlap and sidelap between successive images (i.e., 75%).The photos were externally aligned thanks to a series of ground control points (GCPs) surveyed by the GNSS survey.GCPs related to the nadiral shots were materialized using 12 artificial yellow/black targets (Figure 3A,D).Among them, 9 were used as GCPs and 3 as check points (CP).In addition, 14 points were measured by the RTS on the vertical rock walls (Figure 3B,C): 11 were used as GCPs and 3 were used as CPs.The software used in the UAS photogrammetric data processing is Agisoft Metashape TM (Version 2.0.2) [60], which exploits structure-from-motion and multiview stereo techniques to align the frames-using a bundle block adjustment-and create a 3D dense point cloud.

UAS Photogrammetric Survey
The photogrammetric survey was performed using a DJI TM Mavic2 Pro drone equipped with a Hasselblad L1D-20C, with a 1-inch sensor and 20-megapixel resolution, combined with a 10.26 mm lens.The study area was surveyed by a total of 3 different flights: the first one was executed in automatic mode, with nadiral shooting, using UgCS software (Version 4.18) [59] for flight planning.Two additional flights were carried out in manual mode, with the direction of photo acquisition almost orthogonal to the rock slopes, and the camera tilting from 45 • to 65 • .In this way, it was possible to obtain complete coverage of the external part of the study area from optimal points of view.
A total of 550 images were acquired in an area of 1.9 ha, with a ground spatial resolution of about 1 cm.The flights were planned to ensure sufficient overlap and sidelap between successive images (i.e., 75%).The photos were externally aligned thanks to a series of ground control points (GCPs) surveyed by the GNSS survey.GCPs related to the nadiral shots were materialized using 12 artificial yellow/black targets (Figure 3A,D).Among them, 9 were used as GCPs and 3 as check points (CP).In addition, 14 points were measured by the RTS on the vertical rock walls (Figure 3B,C): 11 were used as GCPs and 3 were used as CPs.The software used in the UAS photogrammetric data processing is Agisoft Metashape TM (Version 2.0.2) [60], which exploits structure-from-motion and multi-view stereo techniques to align the frames-using a bundle block adjustment-and create a 3D dense point cloud.

SLAM-Based LiDAR Survey
The survey of the internal morphology of Grotta Paglicci, conducted by employing the SLAM-based LiDAR mapping technique, was aimed at rebuilding the geometries of the cave.The system used is an Emesent™ Hovermap ST model, which is able to create a 3D dense point cloud that is useful in simultaneously locating, mapping, and navigating within the enclosed chambers, thanks to an inertial navigation system and without the aid of GNSS.
This LiDAR system was set in double-return mode, with an acquisition speed of up to 600,000 points/s, in order to potentially guarantee high levels of ground detail.This technology characteristic proved to be very effective for the surveying of Grotta Paglicci, which is characterized by large internal chambers with a total approximate length of 60 m, confined spaces, corridors, ups and downs, and numerous shaded areas.
For the internal survey of the entire cave, a total of 3 runs were executed.The coregistration of the 3D dense point clouds was performed using the iterative closest point (ICP) algorithms within Trimble TM Realworks software (Version 12.4) [61].
In addition, the Emesent™ Hovermap ST was equipped with an action camera that recorded videos simultaneously with the runs' acquisition; this was particularly useful for the attribution of RGB colors to the interior point cloud of the cave.To ensure adequate lighting conditions during the laser scanning acquisitions, the cave chambers were fully illuminated using LED lamps (Figure 4).

SLAM-Based LiDAR Survey
The survey of the internal morphology of Grotta Paglicci, conducted by employing the SLAM-based LiDAR mapping technique, was aimed at rebuilding the geometries of the cave.The system used is an Emesent™ Hovermap ST model, which is able to create a 3D dense point cloud that is useful in simultaneously locating, mapping, and navigating within the enclosed chambers, thanks to an inertial navigation system and without the aid of GNSS.
This LiDAR system was set in double-return mode, with an acquisition speed of up to 600,000 points/s, in order to potentially guarantee high levels of ground detail.This technology characteristic proved to be very effective for the surveying of Grotta Paglicci, which is characterized by large internal chambers with a total approximate length of 60 m, confined spaces, corridors, ups and downs, and numerous shaded areas.
For the internal survey of the entire cave, a total of 3 runs were executed.The coregistration of the 3D dense point clouds was performed using the iterative closest point (ICP) algorithms within Trimble TM Realworks software (Version 12.4) [61].
In addition, the Emesent™ Hovermap ST was equipped with an action camera that recorded videos simultaneously with the runs' acquisition; this was particularly useful for the attribution of RGB colors to the interior point cloud of the cave.To ensure adequate lighting conditions during the laser scanning acquisitions, the cave chambers were fully illuminated using LED lamps (Figure 4).

Engineering-Geological Survey and Rock Mass Classification
The engineering-geological survey was carried out in the rock outcrops close to the entrance of the Grotta Paglicci, considering accessible and safe areas (Figure 5), both inside and outside the cave, which were considered important for characterizing the joints and bedding and to classify the rock mass.The collected data refer to the limestones of the "Monte Calvo del Gargano" member and includes discontinuities orientation, spacing, persistence, aperture, the type and thickness of infilling, roughness (using Barton's comb and JRC-joint roughness coefficient-classes), weathering, humidity, and indirect uniaxial compressive strength (using Schmidt's hammer-rebound R-value or sclerometer).

Engineering-Geological Survey and Rock Mass Classification
The engineering-geological survey was carried out in the rock outcrops close to the entrance of the Grotta Paglicci, considering accessible and safe areas (Figure 5), both inside and outside the cave, which were considered important for characterizing the joints and bedding and to classify the rock mass.The collected data refer to the limestones of the "Monte Calvo del Gargano" member and includes discontinuities orientation, spacing, persistence, aperture, the type and thickness of infilling, roughness (using Barton's comb and JRC-joint roughness coefficient-classes), weathering, humidity, and indirect uniaxial compressive strength (using Schmidt's hammer-rebound R-value or sclerometer).
In the selection of survey sites, in addition to areas of logistical and anthropological interest for future excavation and analyses, different slope orientations were also considered in such a way as to be later adopted to assess their stability.Due to the external morphology of the area, which is characterized by high and stepped vertical and sub-vertical rock walls, in addition to the presence of vegetation and unstable blocks, it was decided to conduct additional surveys by photointerpreting georeferenced photos and 3D data from the UAS photogrammetry and laser scanning analyses.The photointerpretation, carried out with CloudCompare software (Version 2) [62], included the measurement of discontinuity orientation (Figure 6), spacing, and persistence.All the collected data were processed and were used to identify the main discontinuity systems and to characterize the rock mass according to the rock mass rating In the selection of survey sites, in addition to areas of logistical and anthropological interest for future excavation and analyses, different slope orientations were also considered in such a way as to be later adopted to assess their stability.
Due to the external morphology of the area, which is characterized by high and stepped vertical and sub-vertical rock walls, in addition to the presence of vegetation and unstable blocks, it was decided to conduct additional surveys by photointerpreting georeferenced photos and 3D data from the UAS photogrammetry and laser scanning analyses.The photointerpretation, carried out with CloudCompare software (Version 2) [62], included the measurement of discontinuity orientation (Figure 6), spacing, and persistence.
In the selection of survey sites, in addition to areas of logistical and anthropological interest for future excavation and analyses, different slope orientations were also considered in such a way as to be later adopted to assess their stability.Due to the external morphology of the area, which is characterized by high and stepped vertical and sub-vertical rock walls, in addition to the presence of vegetation and unstable blocks, it was decided to conduct additional surveys by photointerpreting georeferenced photos and 3D data from the UAS photogrammetry and laser scanning analyses.The photointerpretation, carried out with CloudCompare software (Version 2) [62], included the measurement of discontinuity orientation (Figure 6), spacing, and persistence.All the collected data were processed and were used to identify the main discontinuity systems and to characterize the rock mass according to the rock mass rating All the collected data were processed and were used to identify the main discontinuity systems and to characterize the rock mass according to the rock mass rating or RMR method [63], the Q-slope method [64], and the slope mass rating or SMR method [65].
These empirical rock engineering methods are commonly utilized for the definition of rock mass characteristics, the choice of appropriate support and reinforcement, and the performance assessment of natural or man-made slopes.
The RMR method allows the characterization of the rock mass and, from the empirical formulas proposed in [63] when laboratory tests are not available, to derive the rock mass friction angle (φ), cohesion (c), and deformation modulus (E).Moreover, the data gained from the RMR method represent a fundamental input for the failure criterion adopted in numerical modeling (see Section 4.5.2).
The SMR method is specifically designed to assess the stability of rock slopes, focusing on parameters relevant to their stability, including slope height, slope angle, orientation and the persistence of discontinuities, the presence of water, and rock mass strength.
The Q-Slope method also considers the environmental and geological conditions that can affect slopes over a long period of time and allows the identification of the slope limit dip angle for stability.In addition, it is possible to derive the required and acceptable slope dip angles equal to a 1%, 15%, 30%, and 50% probability of failure (PoF).
In view of their peculiarities, the reason behind the choice of utilizing all these three methods resides in the purpose of comprehensively characterizing the slope stability of Grotta Paglicci.

Rock Slope Stability Analysis
In order to assess the slope stability, we chose to apply two methods: kinematic stability analysis, a conventional and traditional analytical method, and 2D FEM numerical modeling.
Kinematic stability analysis focuses on geometric criteria, friction angle, and possible failure mechanisms, while 2D FEM modeling considers the stress and deformation behavior of the slope under various loading conditions.When analyzing slope stability, FEM analysis is widely utilized to investigate continuous rock masses that are affected by discontinuities [66,67].

Kinematic Stability Analysis
The kinematic stability analysis is based on the Markland test [68], which represents a stereographic technique used to analyze the possibility of failure in rocky slopes due to the following kinematic mechanisms: planar sliding, wedge sliding, and direct, oblique, and flexural toppling [69].
The analysis was carried out using RocScience TM Inc. Dips software (Version 8.022) [70], which is based on the methods of Goodman [71] and Hudson and Harrison [72].This software allows one to analyze and visualize orientation-based geological data and to perform kinematic stability analysis, which considers the slope orientation and the discontinuities' attitude and properties (i.e., friction angle).
The data collected by the engineering-geological in-situ survey, and the photo interpreted with the 3D dense point cloud, provided an optimal database for the statistical kinematic stability analysis, which included the rock mass friction angle and the different slopes.This choice was made while considering the average slopes' natural orientation and their location in proximity to pivotal areas such as the entrance of Grotta Paglicci and the external rock shelter (Figure 2).The selected natural slopes were also utilized in the SMR method calculations.

Numerical Modeling
The geotechnical modeling of Grotta Paglicci was carried out utilizing RS2 TM software (Version 11.019) [73], in addition to the results of the kinematic stability analysis, in order to compute the rock mass principal stresses, probable displacements, strength factors, and strains.In this case study, the rock mass can be considered as a slightly jointed rock mass, with three main joint systems plus the bedding.Generally, in the past, it was always preferable to utilize discrete element methods for these discontinuous rock masses' media, but recently, it has proved effective to utilize FEM numerical modeling as an equivalent method [67].
The co-registration between the external and internal 3D point clouds of the cave allowed the evaluation of the morphological characteristics needed for the selection of the cross-sections applied in the numerical modeling.The cross-sections were created with CloudCompare freeware software, using the 3D point cloud.After the production of the relative photoplanes, cross-sections were digitized in CAD environment and, subsequently, imported into RS2™ software (Version 11.019) [73] in DXF format.The geometric features, derived from the topographic and geological data, were processed along the AA', BB', and CC' cross-sections, the traces of which are shown in red in Figure 7. Their orientation is not related to the slopes considered in the SMR calculations and kinematic stability analysis but was chosen to analyze the stability conditions along pivotal directions inside the cave.
to compute the rock mass principal stresses, probable displacements, strength factors, and strains.In this case study, the rock mass can be considered as a slightly jointed rock mass, with three main joint systems plus the bedding.Generally, in the past, it was always preferable to utilize discrete element methods for these discontinuous rock masses' media, but recently, it has proved effective to utilize FEM numerical modeling as an equivalent method [67].
The co-registration between the external and internal 3D point clouds of the cave allowed the evaluation of the morphological characteristics needed for the selection of the cross-sections applied in the numerical modeling.The cross-sections were created with CloudCompare freeware software, using the 3D point cloud.After the production of the relative photoplanes, cross-sections were digitized in CAD environment and, subsequently, imported into RS2™ software (Version 11.019) [73] in DXF format.The geometric features, derived from the topographic and geological data, were processed along the AA', BB', and CC' cross-sections, the traces of which are shown in red in Figure 7. Their orientation is not related to the slopes considered in the SMR calculations and kinematic stability analysis but was chosen to analyze the stability conditions along pivotal directions inside the cave.The 2D modeling process involved several stages: firstly, a control point grid was created, based on the acquired topographic data.This control point grid was then used to create a surface mesh representing the cave's general shape.The mesh was subsequently subdivided into smaller elements in order to increase the accuracy in correspondence with the areas of greatest interest.Once the mesh was created, the appropriate geotechnical parameters, including the strength of the rock material, the internal friction angle, and the material cohesion, were applied to represent the geological features of the cave.Assigning geotechnical properties to a 2D numerical model is necessary to accurately determine how The 2D modeling process involved several stages: firstly, a control point grid was created, based on the acquired topographic data.This control point grid was then used to create a surface mesh representing the cave's general shape.The mesh was subsequently subdivided into smaller elements in order to increase the accuracy in correspondence with the areas of greatest interest.Once the mesh was created, the appropriate geotechnical parameters, including the strength of the rock material, the internal friction angle, and the material cohesion, were applied to represent the geological features of the cave.Assigning geotechnical properties to a 2D numerical model is necessary to accurately determine how the material deforms, how it resists shear stresses, and how it interacts with water and other materials, reflecting real-world conditions as closely as possible.Furthermore, it is crucial for obtaining reliable predictions of slope stability, deformation, and failure mechanisms.
All the parameters that were assigned to the model material are summarized in Table 1.The values related to unit weight (MN/m 3 ), Poisson's ratio, and tensile strength were obtained by consulting previous local technical documentation involving the same lithology [74,75].Meanwhile, the values for cohesion (c), friction angle (φ), and Young's modulus (E) were derived from theoretical calculations proposed by the authors of [63].The remaining data concerning material behavior, porosity value, and static water mode were chosen based on conditions observed during the engineering-geological survey.It must be mentioned that the models are entirely included within the limestones of the "Monte Calvo del Gargano" member.Throughout the numerical simulations, the RS2 TM software (Version 11.019) [73] calculated the behavior of the slope under different conditions, such as structural loading, ground movements, and simply the self-weight of the materials.Since there were no direct measurements of in situ stress conditions, the analysis was carried out considering only the intrinsic characteristics of the material and the discontinuities involved, thus relying exclusively on the material's self-weight (i.e., body force only, as shown in Table 1).The green color refers to the limestones of Monte Calvo del Gargano member as shown in the numerical models.
Figures 8-10 show the geotechnical models created with RS2™ software (Version 11.019) [73] using the geometric information derived from the AA', BB', and CC' sections, respectively.
GeoHazards 2024, 5 467 of 484 the material deforms, how it resists shear stresses, and how it interacts with water and other materials, reflecting real-world conditions as closely as possible.Furthermore, it is crucial for obtaining reliable predictions of slope stability, deformation, and failure mechanisms.
All the parameters that were assigned to the model material are summarized in Table 1.The values related to unit weight (MN/m 3 ), Poisson's ratio, and tensile strength were obtained by consulting previous local technical documentation involving the same lithology [74,75].Meanwhile, the values for cohesion (c), friction angle (φ), and Young's modulus (E) were derived from theoretical calculations proposed by the authors of [63].The remaining data concerning material behavior, porosity value, and static water mode were chosen based on conditions observed during the engineering-geological survey.It must be mentioned that the models are entirely included within the limestones of the "Monte Calvo del Gargano" member.Throughout the numerical simulations, the RS2 TM software (Version 11.019) [73] calculated the behavior of the slope under different conditions, such as structural loading, ground movements, and simply the self-weight of the materials.Since there were no direct measurements of in situ stress conditions, the analysis was carried out considering only the intrinsic characteristics of the material and the discontinuities involved, thus relying exclusively on the material's self-weight (i.e., body force only, as shown in Table 1).
Figures 8-10 show the geotechnical models created with RS2™ software (Version 11.019) [73] using the geometric information derived from the AA', BB', and CC' sections, respectively.The green color refers to the limestones of Monte Calvo del Gargano member as shown in the numerical models.
The created models include the main discontinuities and their orientation, spacing, and persistence, as computed by combining data from the traditional field survey with the interpretation of geomatic data.The latter is derived from both UAS photogrammetry and SLAM-based laser scanning and refers to both external areas and internal cave chambers.The Barton-Bandis failure criterion [76], which has been widely utilized and extensively studied and validated for various rock types and geological conditions, was chosen for the present slope stability analysis.The selected failure criterion permits the characterization of the discontinuities' shear strength, based on parameters such as the JRC, which is directly derived from the engineering-geological survey.In addition, the joint wall compressive strength (JCS), the residual friction angle, and the normal and shear The created models include the main discontinuities and their orientation, spacing, and persistence, as computed by combining data from the traditional field survey with the interpretation of geomatic data.The latter is derived from both UAS photogrammetry and SLAM-based laser scanning and refers to both external areas and internal cave chambers.The Barton-Bandis failure criterion [76], which has been widely utilized and extensively studied and validated for various rock types and geological conditions, was chosen for the present slope stability analysis.The selected failure criterion permits the characterization of the discontinuities' shear strength, based on parameters such as the JRC, which is directly derived from the engineering-geological survey.In addition, the joint wall compressive strength (JCS), the residual friction angle, and the normal and shear stiffnesses were, respectively, derived using empirical formulas from [77][78][79].
The required engineering-geological characteristics assigned to the discontinuities are summarized in Table 2.The FEM analysis for the models created along the three cross-sections was performed according to a certain number of iterations and a maximum tolerance value up to convergence.

SLAM-Based Point Cloud Coregistration and Georeferencing
Figure 11 shows the nadiral view of the georeferenced 3D point cloud obtained from the coregistration of the three SLAM-based LiDAR surveys.The orange dots (EMn in Figure 11) represent the location of the 4 GCPs (the circular-shaped optical target in Figure 3E) used to georeference the entire survey in the ETRF2000/UTM33N absolute system.The georeferencing accuracy was estimated at about 6 cm.
The colors refer to the different discontinuity systems as shown in the numerical models.

SLAM-Based Point Cloud Coregistration and Georeferencing
Figure 11 shows the nadiral view of the georeferenced 3D point cloud obtained from the coregistration of the three SLAM-based LiDAR surveys.The orange dots (EMn i Figure 11) represent the location of the 4 GCPs (the circular-shaped optical target in Figur 3E) used to georeference the entire survey in the ETRF2000/UTM33N absolute system.Th georeferencing accuracy was estimated at about 6 cm.

UAS Photogrammetric Data Processing
Structure-from-motion and multi-view stereo algorithms allowed the alignment o the images through the identification of 1.99 million recognizable tie points.The root mea square error (RMSE) obtained in the process of exterior orientation was about 2.9 cm fo GCPs and 1.9 cm for CPs.The final 3D dense point cloud comprises 170 million point that were interpolated to generate a continuous model in raster format (i.e., a digital dens elevation model-DDEM) with a spatial resolution of 1.8 cm/pixel.Based on the DDEM the images were ortho-corrected and mosaiced (in both the frontal and nadiral views) a a ground sampling distance (GSD) of 1 cm/pixel (Figure 12).

UAS Photogrammetric Data Processing
Structure-from-motion and multi-view stereo algorithms allowed the alignment of the images through the identification of 1.99 million recognizable tie points.The root mean square error (RMSE) obtained in the process of exterior orientation was about 2.9 cm for GCPs and 1.9 cm for CPs.The final 3D dense point cloud comprises 170 million points that were interpolated to generate a continuous model in raster format (i.e., a digital dense elevation model-DDEM) with a spatial resolution of 1.8 cm/pixel.Based on the DDEM, the images were ortho-corrected and mosaiced (in both the frontal and nadiral views) at a ground sampling distance (GSD) of 1 cm/pixel (Figure 12).

Engineering-Geological Data
The study is characterized by the presence of various discontinuities with highly variable orientations.The information from traditional engineering-geological surveys, integrated with that from geomatic data photointerpretation, permitted the characterization of the rock mass thanks to the abundance of deterministically measured discontinuities (i.e., 364) uniformly distributed across the entire study area, even at different elevations from the ground.The high number of measured data has enabled the identification of the primary discontinuity systems, as shown by the stereographic representation in Figure 13.

Engineering-Geological Data
The study is characterized by the presence of various discontinuities with highly variable orientations.The information from traditional engineering-geological surveys, integrated with that from geomatic data photointerpretation, permitted the characterization of the rock mass thanks to the abundance of deterministically measured discontinuities (i.e., 364) uniformly distributed across the entire study area, even at different elevations from the ground.The high number of measured data has enabled the identification of the primary discontinuity systems, as shown by the stereographic representation in Figure 13.

Engineering-Geological Data
The study is characterized by the presence of various discontinuities with highly variable orientations.The information from traditional engineering-geological surveys, integrated with that from geomatic data photointerpretation, permitted the characterization of the rock mass thanks to the abundance of deterministically measured discontinuities (i.e., 364) uniformly distributed across the entire study area, even at different elevations from the ground.The high number of measured data has enabled the identification of the primary discontinuity systems, as shown by the stereographic representation in Figure 13.Following the method proposed by [63] and explained in Section 4.4, a value of the basic rock mass rating (RMRb) equal to 67 was calculated for the rock mass under study.This value indicates good quality (i.e., "Class II").An indicative estimate of the shear   Following the method proposed by [63] and explained in Section 4.4, a value of the basic rock mass rating (RMRb) equal to 67 was calculated for the rock mass under study.This value indicates good quality (i.e., "Class II").An indicative estimate of the shear strength parameters was provided by the formulas reported in [63], which include the rock mass friction angle (φ), cohesion (c), and Young's modulus (E).The calculated values for this case study are presented in Table 4.
Table 4. Shear strength parameters and the deformation modulus (E), calculated according to [63].

Q-Slope Method
The empirical assessment of the rock mass conditions was also conducted by calculating the Q-slope value, as proposed by [64] and explained in Section 4.4, which resulted in a value equal to 0.81.Considering an average dip of 80 • for the external slope, when applying the graph shown in Figure 15, the slope results in unstable conditions and is, therefore, prone to possible failure.Moreover, Table 5 presents the results obtained from the application of the empirical formulas, as proposed by [64], to estimate the required slope dip angle of the outside topography in the absence of stabilization measures, accepting a certain probability of failure (PoF) equal to 1%, 15%, 30%, and 50%, respectively.
GeoHazards 2024, 5 472 of 484 strength parameters was provided by the formulas reported in [63], which include the rock mass friction angle (φ), cohesion (c), and Young's modulus (E).The calculated values for this case study are presented in Table 4.
Table 4. Shear strength parameters and the deformation modulus (E), calculated according to [63].

Q-Slope Method
The empirical assessment of the rock mass conditions was also conducted by calculating the Q-slope value, as proposed by [64] and explained in Section 4.4, which resulted in a value equal to 0.81.Considering an average dip of 80° for the external slope, when applying the graph shown in Figure 15, the slope results in unstable conditions and is, therefore, prone to possible failure.Moreover, Table 5 presents the results obtained from the application of the empirical formulas, as proposed by [64], to estimate the required slope dip angle of the outside topography in the absence of stabilization measures, accepting a certain probability of failure (PoF) equal to 1%, 15%, 30%, and 50%, respectively.The application of the slope mass rating method proposed by [65] has allowed the additional characterization of the rock mass by considering the geometric relationships between the observed discontinuities and the natural slopes.The analysis of the rock mass quality was verified by considering planar, toppling, and wedge failure phenomena.In Table 6, the average orientations of the discontinuity systems (Sn), plus the selected natural slopes (Vn), which are shown in Figure 16A, are reported.The application of the slope mass rating method proposed by [65] has allowed the additional characterization of the rock mass by considering the geometric relationships between the observed discontinuities and the natural slopes.The analysis of the rock mass quality was verified by considering planar, toppling, and wedge failure phenomena.In Table 6, the average orientations of the discontinuity systems (Sn), plus the selected natural slopes (Vn), which are shown in Figure 16A, are reported.6, considered for the application of the SMR method (A).Examples of the kinematic stability analysis executed on slope V3 through a stereographic projection (Wulff equal-angle methodlower hemisphere): planar sliding (B), wedge sliding (C), and direct and oblique toppling (D).S1, S2 and S3 indicate the different discontinuity systems as shown in Table 6.6, considered for the application of the SMR method (A).Examples of the kinematic stability analysis executed on slope V3 through a stereographic projection (Wulff equal-angle method-lower hemisphere): planar sliding (B), wedge sliding (C), and direct and oblique toppling (D).S1, S2 and S3 indicate the different discontinuity systems as shown in Table 6.
Table 7 shows the SMR results for the three analyzed slopes.The minimum SMR value for planar failure is 57, which corresponds to a "Class III" rock mass of normal quality, which is partially stable, showing sliding along planes or wedges and requiring systematic interventions.The minimum SMR value for toppling is 82 which corresponds to a "Class I" rock mass of very good quality, completely stable, with no instability phenomena, and with no need for stabilization interventions.The SMR values, when analyzing the wedge failure, are generally "good" (Class II) and, in one case, "very good" (Class I), indicating overall stability and no need for stabilization interventions.However, in the case of the intersection between joint systems S1 and S2 (on slope V1), as well as between S1 and S3 (on slopes V1 and V2), the SMR values show a sudden decrease to 22 and 32, respectively.These SMR values correspond to a "Class IV" rock mass of poor quality with instability phenomena along the wedges, which requires extensive remedial measurements.

Statistical Kinematic Stability Analysis
The statistical kinematic stability analysis was based on the discontinuity systems and slopes presented in Table 6 and Figure 16A.The value of the friction angle was set to 38 • accordingly, with the results shown in Table 4 and confirmed by the available literature for the study area [74].In Table 8, which summarizes the results obtained from the statistical kinematic stability analysis, wedge sliding and toppling (both direct and oblique) come out as the predominant possible failure mechanisms.Planar sliding only occurs in the case of slope V3, and it is related to the discontinuity system S1, which acts as the sliding surface.A flexural toppling phenomenon has not been observed, although it is usually considered to be fairly probable with this kind of geological setting.Figure 16B-D shows examples of the kinematic stability analysis executed on slope V3 through stereographic projections for planar sliding, wedge sliding, and direct and oblique toppling, respectively.Table 8. Results from the statistical kinematic stability analysis executed along slopes V1, V2, and V3, as shown in Table 6.

Numerical Modeling
The cross-sections obtained by processing the topographic and geological data, enriched by assigning geotechnical properties to the rock material and discontinuities (i.e., joints and bedding), allowed the creation of geotechnical models that provided a detailed virtual representation of the Grotta Paglicci site.The obtained results allowed for an evaluation of rock mass behavior, providing insights about the localization of stresses along the considered cross-sections.The trend of the field differential stress σ 1 -σ 3 around the cave is shown in Figures 17-19.This value is crucial to understanding the behavior and the structural response of a geotechnical model since it allows us to calculate both the shear and tensile stresses in the presence of jointed rocks and cavities.
GeoHazards 2024, 5 475 of 484 joints and bedding), allowed the creation of geotechnical models that provided a detailed virtual representation of the Grotta Paglicci site.The obtained results allowed for an evaluation of rock mass behavior, providing insights about the localization of stresses along the considered cross-sections.The trend of the field differential stress σ1-σ3 around the cave is shown in Figures 17-19.This value is crucial to understanding the behavior and the structural response of a geotechnical model since it allows us to calculate both the shear and tensile stresses in the presence of jointed rocks and cavities.Regarding the AA' cross-section (Figure 17), the highest values of differential stress reach up to 39 MPa and are primarily localized at the intersection between the bedding, joints, and cave surface.Specifically, the range of maximum values varies from approximately 18 to 39 MPa.
Values of differential stresses along the BB' cross-section (Figure 18) result much lower than those shown in Figure 17, with maximum values of 1.31 MPa in areas where the S2 joint system intersects the cave surface.In this case, the differential stress ranges from a minimum of 0.87 to a maximum of 2.90 MPa at the edges of the cave chambers.
Similar values can be observed along cross-section CC' (Figure 19), where the recorded maximum values of differential stress range from a minimum of 0.72 MPa to a maximum of 1.44 MPa.
In Figure 20, the total displacement values, indicated by the red colors, stand at approximately 5 cm and are located in the roof connecting the Atrium and Chamber 1, where the S1 discontinuity intersects with the limestone bedding, creating potentially unstable geometries.This critical condition was not encountered along cross-sections BB' and CC', which present values of lower than one millimeter; Figure 21, as an example, shows the displacements along cross-section CC'.Regarding the AA' cross-section (Figure 17), the highest values of differential stress reach up to 39 MPa and are primarily localized at the intersection between the bedding, joints, and cave surface.Specifically, the range of maximum values varies from approximately 18 to 39 MPa.
Values of differential stresses along the BB' cross-section (Figure 18) result much lower than those shown in Figure 17, with maximum values of 1.31 MPa in areas where the S2 joint system intersects the cave surface.In this case, the differential stress ranges from a minimum of 0.87 to a maximum of 2.90 MPa at the edges of the cave chambers.
Similar values can be observed along cross-section CC' (Figure 19), where the recorded maximum values of differential stress range from a minimum of 0.72 MPa to a maximum of 1.44 MPa.
In Figure 20, the total displacement values, indicated by the red colors, stand at approximately 5 cm and are located in the roof connecting the Atrium and Chamber 1, where the S1 discontinuity intersects with the limestone bedding, creating potentially unstable geometries.This critical condition was not encountered along cross-sections BB' and CC', which present values of lower than one millimeter; Figure 21, as an example, shows the displacements along cross-section CC'.Regarding the AA' cross-section (Figure 17), the highest values of differential stress reach up to 39 MPa and are primarily localized at the intersection between the bedding, joints, and cave surface.Specifically, the range of maximum values varies from approximately 18 to 39 MPa.
Values of differential stresses along the BB' cross-section (Figure 18) result much lower than those shown in Figure 17, with maximum values of 1.31 MPa in areas where the S2 joint system intersects the cave surface.In this case, the differential stress ranges from a minimum of 0.87 to a maximum of 2.90 MPa at the edges of the cave chambers.
Similar values can be observed along cross-section CC' (Figure 19), where the recorded maximum values of differential stress range from a minimum of 0.72 MPa to a maximum of 1.44 MPa.
In Figure 20, the total displacement values, indicated by the red colors, stand at approximately 5 cm and are located in the roof connecting the Atrium and Chamber 1, where the S1 discontinuity intersects with the limestone bedding, creating potentially unstable geometries.This critical condition was not encountered along cross-sections BB' and CC', which present values of lower than one millimeter; Figure 21, as an example, shows the displacements along cross-section CC'.

Discussion
The integration of several geomatic techniques, together with data from a traditional engineering-geological survey, allowed the production of a large number of results, which are pivotal for the site's geological characterization and the identification of critically hazardous conditions in relation to slope stability.Furthermore, the combined use of UAS photogrammetry and SLAM-based laser scanning proved to be efficient in terms of data accuracy and completeness.In previous papers, as reported in [12][13][14][15][16][17][18][19][20]80], the use of ground-based LiDAR has proven to be very reliable but also, simultaneously, time-consuming in terms of data acquisition and processing.With the advent of SLAM technology, which can be used manually or integrated into UASs or ground robots, laser scanning has undergone further improvements in terms of the number of acquired points, acquisition time, and completeness [26][27][28][29][30].
In the present case study, the use of UASs and ground robotic equipment could not be used within the cave, due to its complexity (tortuous and narrow connections between the various chambers and the presence of speleothems).Conversely, the use of portable SLAM-based LiDAR permitted us to map the four chambers of Grotta Paglicci with an optimal level of resolution with spatial accuracy-thanks also to a thorough topographic survey-and in a relatively short acquisition time (less than one hour).
The RMSEs obtained from the adopted methodologies demonstrate the success of the alignment and georeferencing of the 3D dense point clouds created by UAS photogrammetry and laser scanning.The possibility of having a merged global model has provided an optimal tool, not only to map the visible joints, both internal and external, but also to verify their correspondence in terms of persistence.This aspect should be considered as one of the main features of this case study from a stability analysis and stress computation perspective.Indeed, thanks to the moderate thickness of the cave vault, it was possible to model the effective length of discontinuities, leading to a reliable numerical analysis.
The engineering-geological characterization indicates a rock mass of good quality, but critical issues can occur when the slope dips exceed 60°.These slopes, combined with unfavorable discontinuity attitudes, can lead to possible block and wedge failures and overall stability issues on various slope orientations.These results aligned with field observations, where residual sliding planes and various blocks from past gravitational events were identified at the bottom of the rocky walls.Since the presented methods of rock mass classification are empirical, they were also supported by statistical kinematic stability analysis, executed along three natural rock slopes of average dips.The results of such stability analyses indicate wedge sliding and direct and oblique toppling as the predominant mechanisms, while planar sliding is only possible along slope V3, which appears to be the most strongly affected by critical phenomena.The interactions among the dip direction of this slope, the S1 discontinuity system, and bedding facilitate several possible failing mechanisms.In the case of wedge sliding, S1-S2 and S1-S3 are the most common joint intersections creating instabilities, although there is not a predominant combination in the case of direct and oblique toppling.The outcome of the present

Discussion
The integration of several geomatic techniques, together with data from a traditional engineering-geological survey, allowed the production of a large number of results, which are pivotal for the site's geological characterization and the identification of critically hazardous conditions in relation to slope stability.Furthermore, the combined use of UAS photogrammetry and SLAM-based laser scanning proved to be efficient in terms of data accuracy and completeness.In previous papers, as reported in [12][13][14][15][16][17][18][19][20]80], the use of groundbased LiDAR has proven to be very reliable but also, simultaneously, time-consuming in terms of data acquisition and processing.With the advent of SLAM technology, which can be used manually or integrated into UASs or ground robots, laser scanning has undergone further improvements in terms of the number of acquired points, acquisition time, and completeness [26][27][28][29][30].
In the present case study, the use of UASs and ground robotic equipment could not be used within the cave, due to its complexity (tortuous and narrow connections between the various chambers and the presence of speleothems).Conversely, the use of portable SLAM-based LiDAR permitted us to map the four chambers of Grotta Paglicci with an optimal level of resolution with spatial accuracy-thanks also to a thorough topographic survey-and in a relatively short acquisition time (less than one hour).
The RMSEs obtained from the adopted methodologies demonstrate the success of the alignment and georeferencing of the 3D dense point clouds created by UAS photogrammetry and laser scanning.The possibility of having a merged global model has provided an optimal tool, not only to map the visible joints, both internal and external, but also to verify their correspondence in terms of persistence.This aspect should be considered as one of the main features of this case study from a stability analysis and stress computation perspective.Indeed, thanks to the moderate thickness of the cave vault, it was possible to model the effective length of discontinuities, leading to a reliable numerical analysis.
The engineering-geological characterization indicates a rock mass of good quality, but critical issues can occur when the slope dips exceed 60 • .These slopes, combined with unfavorable discontinuity attitudes, can lead to possible block and wedge failures and overall stability issues on various slope orientations.These results aligned with field observations, where residual sliding planes and various blocks from past gravitational events were identified at the bottom of the rocky walls.Since the presented methods of rock mass classification are empirical, they were also supported by statistical kinematic stability analysis, executed along three natural rock slopes of average dips.The results of such stability analyses indicate wedge sliding and direct and oblique toppling as the predominant mechanisms, while planar sliding is only possible along slope V3, which appears to be the most strongly affected by critical phenomena.The interactions among the dip direction of this slope, the S1 discontinuity system, and bedding facilitate several possible failing mechanisms.In the case of wedge sliding, S1-S2 and S1-S3 are the most common joint intersections creating instabilities, although there is not a predominant combination in the case of direct and oblique toppling.The outcome of the present analysis shows how all the slopes contouring the external area of Grotta Paglicci are prone to kinematic instabilities; therefore, deeper dynamic slope stability studies are in progress to define the most vulnerable walls and to plan further actions of reinforcement and support.
Numerical modeling allowed for a detailed study of the cave's stability and its underground structures.A fundamental aspect addressed through numerical modeling is the assessment of solutions to prevent potential cave-ins or damage to the rock formations inside the cave.Similar research in underground environments [22,81] has proved that the utilization of 2D numerical modeling is a highly effective tool for assessing structural stability and proposing adequate conservation strategies.
Numerical analysis has allowed for the identification of areas where critical issues or failures may occur and could lead to episodes of instability.To critically assess the obtained differential stress values, the method described in [82] can be used; it allows the calculation of the threshold value of the differential stress responsible for the creation of new failures in a slightly fractured rock mass.
The rock mass of Grotta Paglicci can be considered slightly fractured, as confirmed by its rating value (an RMRb equal to 67, Class II, and good quality) and the scarce number of joints intersecting the modeled cross-sections.Moreover, since the uniaxial compressive strength of the limestone belonging to the "Monte Calvo del Gargano member" has an average value of 37 MPa [75], the rock mass may exhibit criticalities in areas where the differential stress exceeds a value of 12.12 MPa according to Equation (1), as proposed by [83,84]: The differential stresses shown in Figures 18 and 19 indicate very low values (i.e., lower than 3 MPa), while Figure 17 indicates that the highest values reach up to 39 MPa in the Atrium.Referring to the information provided in [85], if the differential stress is less than four times the tensile strength the rock, extensional failure may occur.Conversely, if the differential stress is more than four times the tensile strength of the rock, shear failure may occur.The characteristic tensile strength of calcareous rocks, as deduced from various laboratory tests conducted on intact rock specimens, generally falls within the range of 1-10 MPa, with values typically ranging around 5 ± 2 MPa [77,[86][87][88][89][90].According to the information provided in [85], the threshold value to discriminate the different types of failure is equal to 28 MPa.When considering 12.12 MPa as the failure threshold, as presented in Equation (1), it can be concluded that shear failure may occur when σ 1 − σ 3 > 28 MPa, while extensional failure may occur when σ 1 − σ 3 < 28 MPa.The results shown in Figure 17 indicate some zones with possible shear failures that, up to now, have not occurred, probably due to the presence of undetermined rock bridges along the joints.Meanwhile, Figures 18 and 19 indicate only the possible extensional failures that have a low probability of occurring, due to their very low differential stress values.
In general, the stress distribution in various chambers shows how the forces are located on their side walls, resulting in low stress values in the roof.This stress distribution can be explained by comparing the geometrical structure of the hypogeal chamber to a dome: this tridimensional geometry discharges the vertical lithostatic weight, similarly to a horseshoe tunnel's geometry, along the lateral walls [91,92].Nevertheless, the high stress values shown in Figure 17 are not confirmed in Figures 18 and 19, where the cross-sections have been built approximately orthogonal to the AA'.The high stress values shown in the AA' cross-section, especially at the entrance of the cave (i.e., 39 MPa in Figure 17), can probably be explained by the fact that the left lateral wall of the Atrium dome has a thickness that is insufficient to distribute the stresses developed by the roof material's weight.This hypothesis could be confirmed by calibrating the model stresses through the execution of in situ measurements, as described in [82].
The evaluation of displacements along the various sections has provided information about potential criticalities at the transition between the Atrium and Chamber 1. Specifically, after analyzing the results in section AA', it is evident that possible instability phenomena can develop at the interaction between the S1 joint system, the bedding, and the cave roof (Figure 20).This critical result is confirmed by the kinematic stability analysis, which indicates the possibility of planar sliding at the intersection between the S1 joint system and slope V3.The analyzed situation, indeed, satisfies both the conditions for the occurrence of planar sliding: namely, the daylight of the joint on the slope and the dip exceeding the friction angle [68].Thanks to this type of analysis and depending on the type of failure that may develop, it was possible to identify those zones where it will be necessary to install safety measures to avoid risks for the personnel present on site for various activities (e.g., archaeological, anthropological, and geological investigations, and maintenance).
A specific project regarding these safety measures has already been proposed to the dedicated agencies, following the technical recommendations suggested by our findings.In particular, the cross-sections deriving from the topographic data and 3D models were utilized to measure the vertical wall thickness, in order to evaluate the maximum depth for the threaded rod that should be installed on the external multi-step slope.The engineeringgeological data and the resulting kinematic stability analysis led to the definition of areas prone to hazardous instabilities in which the only solution is scaling the rock mass.In addition, an in situ monitoring system will be installed, powered by a solar panel, which will comprise electrical crackmeters and inclinometers with live remote alert systems.

Conclusions
In the ever-evolving field of geospatial data acquisition and 3D modeling, to meet the growing demand for accuracy, speed, and versatility in mapping and surveying, experts are increasingly turning to a multi-technique survey approach.This cutting-edge strategy, as used in this work, leverages the synergies of two powerful technologies: UAS photogrammetry and SLAM-based LiDAR mapping.Each of these techniques offers unique strengths and capabilities, and, when combined, they create a comprehensive tool that can tackle a wide range of surveying and mapping challenges.This holistic approach enables specialists to streamline workflows, enabling, simultaneously, a high level of data accuracy, detail, and reduced fieldwork time and costs.
The comprehensive geomatic multi-technique approach conducted at Grotta Paglicci has significantly advanced our understanding of the complex morphology of the area and it has provided essential data for a thorough engineering-geological site characterization.By leveraging the adopted methodologies, this study achieved a detailed site characterization, unveiling crucial insights into slope stability and potential geological hazards in a challenging cave environment, and overcame the limitations posed by narrow passageways and the presence of speleothems.Merging the external and internal 3D point clouds enabled a holistic view of the cave system, offering a unique perspective on its geological and engineering-geological characteristics and joint persistence.This study's significance lies in its blend of innovation and traditional methods, marrying advanced technology with empirical insights.The resulting comprehensive slope stability analysis not only elucidates the cave's current state but also lays the foundation for future conservation endeavors.
By integrating state-of-the-art tools, this research enriches our understanding of Grotta Paglicci's geological complexities.It establishes a robust framework for heritage preservation that can be replicated in other similar caves or natural shelters.
Additionally, the present multi-technique approach can be used for both static and multitemporal surveys and can also be applied in large artificial underground caves, while achieving an optimal balance between the required survey precision and accuracy, the volume of data to be processed over large areas, and the effectiveness of stress and displacement monitoring.
In conclusion, the presented study exemplifies the power of technological integration, demonstrating the synergy between the traditional approach and contemporary advancements.As we navigate the delicate balance between preservation and exploration, this research can be used to bring to the fore sustainable conservation practices in the realm of archaeological and natural heritage.

Figure 1 .
Figure 1.Tectonic scheme of the study area, modified from [31]: in (A), the red box indicates the position of Grotta Paglicci.Geological map of the study area (B) with 25-meter contour line intervals; the letters in (B) indicate: a-slope cover deposit, b-present alluvial deposit, b2-eluvial-colluvial deposit, bn-terrace alluvial deposit, CBA1-Monte Calvo del Gargano member, CBA2-Borgo Celano member.The inset map, (C), shows the location of Grotta Paglicci in the Gargano Promontory.

Figure 2 .
Figure 2. Cave system of Grotta Paglicci.The archaeological and anthropological findings are indicated by colored stars: the yellow star marks the 12-meter-thick stratigraphy (Middle and Upper Paleolithic sequence); the red star marks the wall paintings (Upper Paleolithic); the blue star marks the external rock shelter (Lower/Middle Paleolithic); the purple star marks the Upper Paleolithic sequence outside the present-day entrance; the green star marks the Upper Paleolithic deposit.

Figure 1 .
Figure 1.Tectonic scheme of the study area, modified from [31]: in (A), the red box indicates the position of Grotta Paglicci.Geological map of the study area (B) with 25-meter contour line intervals; the letters in (B) indicate: a-slope cover deposit, b-present alluvial deposit, b 2 -eluvial-colluvial deposit, b n -terrace alluvial deposit, CBA 1 -Monte Calvo del Gargano member, CBA 2 -Borgo Celano member.The inset map, (C), shows the location of Grotta Paglicci in the Gargano Promontory.

Figure 1 .
Figure 1.Tectonic scheme of the study area, modified from [31]: in (A), the red box indicates the position of Grotta Paglicci.Geological map of the study area (B) with 25-meter contour line intervals; the letters in (B) indicate: a-slope cover deposit, b-present alluvial deposit, b2-eluvial-colluvial deposit, bn-terrace alluvial deposit, CBA1-Monte Calvo del Gargano member, CBA2-Borgo Celano member.The inset map, (C), shows the location of Grotta Paglicci in the Gargano Promontory.

Figure 2 .
Figure 2. Cave system of Grotta Paglicci.The archaeological and anthropological findings are indicated by colored stars: the yellow star marks the 12-meter-thick stratigraphy (Middle and Upper Paleolithic sequence); the red star marks the wall paintings (Upper Paleolithic); the blue star marks the external rock shelter (Lower/Middle Paleolithic); the purple star marks the Upper Paleolithic sequence outside the present-day entrance; the green star marks the Upper Paleolithic deposit.

Figure 2 .
Figure 2. Cave system of Grotta Paglicci.The archaeological and anthropological findings are indicated by colored stars: the yellow star marks the 12-meter-thick stratigraphy (Middle and Upper Paleolithic sequence); the red star marks the wall paintings (Upper Paleolithic); the blue star marks the external rock shelter (Lower/Middle Paleolithic); the purple star marks the Upper Paleolithic sequence outside the present-day entrance; the green star marks the Upper Paleolithic deposit.

Figure 3 .
Figure 3. Stages of the GNSS survey (A); the RTS survey (B); an example of a black/white target for RTS scanning (C); a yellow/black target for the UAS photogrammetric survey (D); a circular-shaped optical target for SLAM-based LiDAR scanning (E).

Figure 3 .
Figure 3. Stages of the GNSS survey (A); the RTS survey (B); an example of a black/white target for RTS scanning (C); a yellow/black target for the UAS photogrammetric survey (D); a circular-shaped optical target for SLAM-based LiDAR scanning (E).

Figure 4 .
Figure 4. Stages of the SLAM-based LiDAR survey at Grotta Paglicci.

Figure 4 .
Figure 4. Stages of the SLAM-based LiDAR survey at Grotta Paglicci.

Figure 5 .
Figure 5. Traces of scan lines showing where the engineering-geological surveys were conducted.

Figure 6 .
Figure 6.Example of various discontinuity orientations as measured on the slope, shown in terms of dip and dip direction, with a right-hand rule (A).Detail of the 3D dense point cloud, with the planes for joint orientation measurement highlighted in green (B).

Figure 5 .
Figure 5. Traces of scan lines showing where the engineering-geological surveys were conducted.

Figure 5 .
Figure 5. Traces of scan lines showing where the engineering-geological surveys were conducted.

Figure 6 .
Figure 6.Example of various discontinuity orientations as measured on the slope, shown in terms of dip and dip direction, with a right-hand rule (A).Detail of the 3D dense point cloud, with the planes for joint orientation measurement highlighted in green (B).

Figure 6 .
Figure 6.Example of various discontinuity orientations as measured on the slope, shown in terms of dip and dip direction, with a right-hand rule (A).Detail of the 3D dense point cloud, with the planes for joint orientation measurement highlighted in green (B).

Figure 7 .
Figure 7. New topographic map of the study area from the processing of photogrammetric and SLAM-based LiDAR mapping techniques.Red letters and dotted lines represent the cross-sections used in the 2D numerical modeling.

Figure 7 .
Figure 7. New topographic map of the study area from the processing of photogrammetric and SLAM-based LiDAR mapping techniques.Red letters and dotted lines represent the cross-sections used in the 2D numerical modeling.

Figure 8 .
Figure 8. Geotechnical AA' cross-section; the green and yellow lines represent the main identified joints, while the gray lines indicate the limestone bedding.

Figure 8 .
Figure 8. Geotechnical AA' cross-section; the green and yellow lines represent the main identified joints, while the gray lines indicate the limestone bedding.

Figure 8 .
Figure 8. Geotechnical AA' cross-section; the green and yellow lines represent the main identified joints, while the gray lines indicate the limestone bedding.

Figure 9 .
Figure 9. Geotechnical BB' cross-section; the yellow and magenta lines represent the main identified joints, while the gray lines indicate the limestone bedding.

Figure 9 .Figure 10 .
Figure 9. Geotechnical BB' cross-section; the yellow and magenta lines represent the main identified joints, while the gray lines indicate the limestone bedding.GeoHazards 2024, 5 468 of 484

Figure 10 .
Figure 10.Geotechnical CC' cross-section; the yellow and magenta lines represent the main identified joints, while the gray lines indicate the limestone bedding.

Figure 11 .
Figure 11.Nadiral view of the georeferenced 3D point cloud obtained from the entire SLAM-base LiDAR survey.The orange dots (EMn) represent the location of the GCPs used for point clou georeferencing.

Figure 11 .
Figure 11.Nadiral view of the georeferenced 3D point cloud obtained from the entire SLAM-based LiDAR survey.The orange dots (EMn) represent the location of the GCPs used for point cloud georeferencing.

Figure 12 .
Figure 12.Nadiral view of the georeferenced 3D dense point cloud (A); the blue flags represent the location of GCPs used for the images' exterior orientation.Nadiral orthophotomosaic (B); detail of the frontal orthophotomosaic (C).

Figure 13 .
Figure 13.Stereogram of the discontinuities measured during the in situ engineering-geological survey and by the photointerpretation of geomatic data within CloudCompare software (Version 2) [62].Data are shown in terms of dip and dip direction, with a right-hand rule, using the Schmidt equal-area method (lower hemisphere).Gray shades represent different density concentrations, Figure 14 displays the map of the main joints, categorized by systems, as interpreted from the geomatic data (i.e., 3D dense point cloud, DDEM, and orthophotomosaics).

Figure 12 .
Figure 12.Nadiral view of the georeferenced 3D dense point cloud (A); the blue flags represent the location of GCPs used for the images' exterior orientation.Nadiral orthophotomosaic (B); detail of the frontal orthophotomosaic (C).

GeoHazards 2024, 5 470 of 484 Figure 12 .
Figure 12.Nadiral view of the georeferenced 3D dense point cloud (A); the blue flags represent the location of GCPs used for the images' exterior orientation.Nadiral orthophotomosaic (B); detail of the frontal orthophotomosaic (C).

Figure 13 .
Figure 13.Stereogram of the discontinuities measured during the in situ engineering-geological survey and by the photointerpretation of geomatic data within CloudCompare software (Version 2) [62].Data are shown in terms of dip and dip direction, with a right-hand rule, using the Schmidt equal-area method (lower hemisphere).Gray shades represent different density concentrations, Figure 14 displays the map of the main joints, categorized by systems, as interpreted from the geomatic data (i.e., 3D dense point cloud, DDEM, and orthophotomosaics).

Figure 13 .
Figure 13.Stereogram of the discontinuities measured during the in situ engineering-geological survey and by the photointerpretation of geomatic data within CloudCompare software (Version 2) [62].Data are shown in terms of dip and dip direction, with a right-hand rule, using the Schmidt equal-area method (lower hemisphere).Gray shades represent different density concentrations.

Figure 14
Figure 14 displays the map of the main joints, categorized by systems, as interpreted from the geomatic data (i.e., 3D dense point cloud, DDEM, and orthophotomosaics).

Figure 14 .
Figure 14.Map of the joints as interpreted from the geomatic data.

Figure 14 .
Figure 14.Map of the joints as interpreted from the geomatic data.

Figure 15 .
Figure 15.Stability chart of the Q-slope method (modified from [64]); the black cross indicates the critical stability conditions for an average slope angle of 80°.

Figure 15 .
Figure 15.Stability chart of the Q-slope method (modified from [64]); the black cross indicates the critical stability conditions for an average slope angle of 80 • .

Figure 16 .
Figure 16.Nadiral orthophotomosaic showing the orientations of the three slopes, V1, V2, and V3, in Table6, considered for the application of the SMR method (A).Examples of the kinematic stability analysis executed on slope V3 through a stereographic projection (Wulff equal-angle methodlower hemisphere): planar sliding (B), wedge sliding (C), and direct and oblique toppling (D).S1, S2 and S3 indicate the different discontinuity systems as shown in Table6.

Figure 16 .
Figure 16.Nadiral orthophotomosaic showing the orientations of the three slopes, V1, V2, and V3, in Table6, considered for the application of the SMR method (A).Examples of the kinematic stability analysis executed on slope V3 through a stereographic projection (Wulff equal-angle method-lower hemisphere): planar sliding (B), wedge sliding (C), and direct and oblique toppling (D).S1, S2 and S3 indicate the different discontinuity systems as shown in Table6.

Figure 17 .
Figure 17.Differential stresses resulting from the geotechnical modeling along the AA' cross-section (A).Details of the obtained differential stress values (B).

Figure 18 .
Figure 18.Differential stresses resulting from the geotechnical modeling along the BB' cross-section (A).Detail of the obtained differential stress values (B).

Figure 17 .
Figure 17.Differential stresses resulting from the geotechnical modeling along the AA' cross-section (A).Details of the obtained differential stress values (B).

Figure 17 .
Figure 17.Differential stresses resulting from the geotechnical modeling along the AA' cross-section (A).Details of the obtained differential stress values (B).

Figure 18 .
Figure 18.Differential stresses resulting from the geotechnical modeling along the BB' cross-section (A).Detail of the obtained differential stress values (B).Figure 18. Differential stresses resulting from the geotechnical modeling along the BB' cross-section (A).Detail of the obtained differential stress values (B).

Figure 18 .
Figure 18.Differential stresses resulting from the geotechnical modeling along the BB' cross-section (A).Detail of the obtained differential stress values (B).Figure 18. Differential stresses resulting from the geotechnical modeling along the BB' cross-section (A).Detail of the obtained differential stress values (B).

Figure 19 .
Figure 19.Differential stresses resulting from the geotechnical modeling along the CC' cross-section (A).Detail of the obtained differential stress values (B).

Figure 20 .
Figure 20.Displacement values resulting from geotechnical modeling along cross-section AA'.

Figure 19 .
Figure 19.Differential stresses resulting from the geotechnical modeling along the CC' cross-section (A).Detail of the obtained differential stress values (B).

GeoHazards 2024, 5 476 of 484 Figure 19 .
Figure 19.Differential stresses resulting from the geotechnical modeling along the CC' cross-section (A).Detail of the obtained differential stress values (B).

Figure 20 .
Figure 20.Displacement values resulting from geotechnical modeling along cross-section AA'.Figure 20.Displacement values resulting from geotechnical modeling along cross-section AA'.

Figure 20 .
Figure 20.Displacement values resulting from geotechnical modeling along cross-section AA'.Figure 20.Displacement values resulting from geotechnical modeling along cross-section AA'.

Figure 21 .
Figure 21.Displacement values resulting from geotechnical modeling along cross-section CC'.

Figure 21 .
Figure 21.Displacement values resulting from geotechnical modeling along cross-section CC'.

Table 1 .
Engineering-geological characteristics of the modeled rock material.

Table 1 .
Engineering-geological characteristics of the modeled rock material.

Table 2 .
Engineering-geological characteristics of the discontinuity systems.

Table 3 lists
the engineering-geological parameters, organized according to the various discontinuity systems.

Table 3 .
Engineering-geological parameters of the identified discontinuity systems.

Table 3 lists
the engineering-geological parameters, organized according to the various discontinuity systems.

Table 3 .
Engineering-geological parameters of the identified discontinuity systems.

Table 6 .
Average values of the dip direction and dip of discontinuity systems and slopes used to verify the possible criticalities.

Table 6 .
Average values of the dip direction and dip of discontinuity systems and slopes used to verify the possible criticalities.

Table 7 .
Results of the SMR method application with reference to the three analyzed slopes, V1, V2, and V3, as shown in Table6.
* Superscript stars indicate the resulting critical values.