High‑Resolution 3D FEM Stability Analysis of the Sabereebi Cave Monastery, Georgia

This study assesses the static stability of the artificial Sabereebi Cave Monastery southeast of Georgia's capital, Tbilisi. The cliff into which these Georgian-Orthodox caverns, chapels, and churches were carved consists of a five-layered sequence of weak sedimentary rock—all of which bear a considerable failure potential and, consequently, pose the challenge of preservation to geologists, engineers, and archaeologists. In the first part of this study, we present a strategy to process point cloud data from drone photogrammetry as well as from laser scanners acquired in- and outside the caves into high-resolution CAD objects that can be used for numerical modeling ranging from macro- to micro-scale. In the second part, we explore four distinct series of static elasto-plastic finite element stability models featuring different levels of detail, each of which focuses on specific geomechanical scenarios such as classic landsliding due to overburden, deformation of architectural features as a result of stress concentration, material response to weathering, and pillar failure due to vertical load. With this bipartite approach, the study serves as a comprehensive 3D stability assessment of the Sabereebi Cave Monastery on the one hand; on the other hand, the established procedure should serve as a pilot scheme, which could be adapted to different sites in the future combining non-invasive and relatively cost-efficient assessment methods, data processing and hazard estimation. Strategy application to a real-life geoarchaeological case study Demonstration of versatile FEM model usage for different geotechnical Failure potential underground compound of seven caves and sub-caves


Introduction
The Sabereebi Monastery is an artificially carved cave complex in the Kakheti Region in Eastern Georgia, some 60 km southeast of the capital Tbilisi (41° 28′ 15.17″ N, 45° 33′ 39.42″ E; Fig. 1a-b). Seven caves with their associated sub-caves are located inside one of the numerous northwest-southeast trending ridge outcrops, which are characteristic of the Iori Upland-an intermountain Plateau within the Iveria Plain raising to 1000 masl from the depression between the Greater and the Lesser Caucasus (Fig. 1b). The Iori Upland displays a hilly relief exposing anti-and synclinal hummocks and ridges of tectonoerosive origin with relative heights of up to 300 m above the surrounding flatlands and lengths of several kilometers (Fig. 2;Javakhishvili et al. 2019;Tielidze et al. 2019a, b). The uplift rate is estimated at 2-3 mm/year (Gobejishvili 2011). Lithologic sequences are gently inclined dipping towards the northeast and consist of continental and marine sandstone, conglomerate, limestone, marl, and clay dating to Miocene up to Pleistocene ages; Quaternary sediment covers in the form of alluvial, diluvial, and proluvial deposits are frequent (Gamkrelidze 1992;Gobejishvili and Tsereteli 2012;Tielidze et al. 2019b;Tsereteli 1964).
The general basin morphologies of the region are asymmetric as a result of reverse faulting with ridges displaying low-gradient northeast-facing flanks but quasi-vertical cliffs on their southwest faces which transition to more gentle gradients at the cliff foot. The latter are frequently dissected by local small-scale channels which bundle earth and debris flows and, thus, play a key role in erosive processes. In this context, it is conspicuous that ridge concavities predisposing such dissecting channels across the cliffs are often accompanied-if not yet conditioned-by the channel network on the northeast-facing flank of the ridge.
Although drained by the Mtkvari River (also: Kura River)-the main drain in the Intra-Caucasus-Depression towards the Caspian Sea -and the Iori River as left side first-order tributary, the Iori Upland is characterized by an arid and moderate to warm continental half-desert steppe climate, strong winds, distinct temperature contrasts reaching seasonal differences of up to 60 °C, and daily air temperature averages of 10-12 °C (Gagua and Mumladze 2012;Kordzakhia 1964;Mumladze and Lomidze 2012;Tielidze et al. 2019b, c). Precipitation rates are low, amounting to application list for the UNESCO World Heritage (UNESCO 2021). Non-invasive permanent alongside temporary monitoring systems such as ground-based synthetic aperture radar (GBSAR) interferometry Margottini et al. 2016a), close-range digital orthophotogrammetry (Frodella et al. 2020;Kirkitadze et al. 2015Kirkitadze et al. , 2016Spizzichino et al. 2017), differential Global Positioning System (DGPS) methods (Okrostsvaridze et al. 2016), 3D terrestrial laser scanning (TLS; Margottini et al. 2015bMargottini et al. , 2016aMargottini et al. , 2017 and infrared thermography (Frodella et al. 2020;Spizzichino et al. 2017) have been employed to detect displacements, fracture formation, and drainage paths. Moreover, different stationary systems have been installed for the assessment of local meteorological conditions, microclimate, microtremors, and ambient seismic noise Basilaia et al. 2016) as the region is seismically active. In 1283, about two thirds of the holed cliff collapsed during the Samtskhe Earthquake (M s ≈ 7.0), and further severe damage resulting from strong earthquakes within the cavities and across the slope is reported over time (Godoladze et al. 2016;Korzhenkov et al. 2017). Geomechanical tests, single block and local failure analysis, including landslide risk estimation, were conducted (Boldini et al. 2018;Margottini et al. 2012Margottini et al. , 2015aMargottini et al. , 2016bSpizzichino et al. 2017), and Margottini et al. (2016c) estimated the rockfall potential via a mapping approach across the cave city of Vardzia.
Despite the various high-resolution methods applicable to such geoarchaeological sites, none of them resulted, though, in a complete 3D geotechnical numerical model-a shortcoming that we tackled in this work via a multi-stage strategy covering the entire range from 3D point cloud measurements to a comprehensive finite element method (FEM) model. We chose the Sabereebi Cave Monastery as a case study due to the prevailing Fig. 2 Part of Iori Upland with its hilly relief exposing tectonoerosive hummocks and ridges. The Sabereebi Cave Monastery is located in the cliff in the right lower corner (picture as part of the drone survey in 2019). The picture is taken roughly towards the north; scale relations are given by Fig. 1a need for detailed geotechnical assessment in order to ensure the intactness and conservation of the site on the one hand, and due to the manageable size of the structure of interest regarding numerical modeling on the other hand. It should be noted that this study serves as a pilot scheme, which could be adapted in the future for different sites-e.g., for the cave city of Vardzia or the Uplistsikhe and David Gareja Monastery Complexes (Fig. 1b)-as the overall geotechnical problem and the desired type of resulting information would remain the same.

Data
Digital raw data consist predominantly of high-resolution point clouds acquired by a RIEGL© VZ-1000 3D Terrestrial Laser Scanner in July 2019. In total, scans were taken from 21 stations, of which three were located outside the cliff and the other 18 at different stations inside the seven caves and sub-caves ( Fig. 1a) to capture the maximum of carved structures within the rock mass. At each station outside the cliff, two series of scans were conducted: first, a general rotation by 360° around the vertical axis of the laser scanner in low-resolution (i.e., 0.02° horizontally and 0.02° vertically), and second, a rotation segment focusing on a cliff panorama in high-resolution (i.e., 0.01° horizontally and 0.009° vertically) using the same vertical axis. Within the caves and sub-caves, only general rotations in low-resolution were sufficient to assess the dimensions of their voids. All scans were adjusted to each other using the Iterative Closest Point (ICP) Method (Besl and McKay 1992) implemented in the RIEGL© in-house software RiSCAN PRO and subsequently underwent Octree Sampling (CloudCompare 2022) with an octree level of 21 in CloudCompare to simplify and homogenize the composite of all individual scans in terms of exact position within the World Geodetic System 1984 (WGS 84) and the Universal Transverse Mercator (UTM) Zone 38 N.
The backward inclined top of the cliff (Fig. 2) was assessed via unmanned aerial vehicle (UAV) orthomosaic photogrammetry in November 2019 using a DJI© Mavic 2 Pro Drone equipped with a HASSELBLAD© L1D-20c Camera (Fig. 1a). In total, the mapped area of about 180,000 m 2 was covered by 243 pictures with an approximate overlap of 70-80%, ensuring a proper digital terrain model (DTM) generation with Agisoft Metashape (Agisoft LLC 2021). Orientation was provided via control points, which could be re-identified on the photographs, combined with DGPS and a system of continuously operating reference stations (CORS), allowing for sub-centimeter resolution within the horizontal coordinates. The obtained DTM was likewise adjusted to the geographic reference used for the laser scanner data.
Moreover, several series of photographic materials dating to the years 2018, 2019, and 2020 were used to track damage over time. Rock samples were collected solely from fallen debris close to the caves entrances, as it is not suggested to actively extract samples out of a geoarchaeological site.

Methodology and Model Establishment
As mentioned above, this work satisfies the urgency of a comprehensive static assessment of the Sabereebi Cave Monastery but simultaneously presents a new strategy of establishing 3D numerical models from point clouds. In this section, particular attention is, therefore, paid to our methodological reverse engineering approach of data processing and the creation of different types of FEM models relying thereon.

Mesh Generation from Laser Scanner and Drone Data
The first major part of digital data processing targets the individual point clouds acquired by the laser scanner at different locations inside the seven caves. Foremost, the individual cave and sub-cave point clouds were cleaned in MeshLab from artifacts (i.e., outliers mainly caused by acquisition errors) not belonging to the respective geometries without accidentally disturbing the high level of detail, and the resulting point clouds were resampled to adjust the level of resolution and reduce the noise using Poisson-Disk Sampling (Corsini et al. 2012). An essential step for the (re-)construction of a surface mesh bridging zones of point lack is the computation of a preferably unidirectional normal vector field (i.e., either outwards from or inwards to the cave center), which was also carried out in MeshLab before creating a preliminary mesh through the Screened Poisson Method (Kazhdan and Hoppe 2013).
As manual cleaning "from outside through the cave walls" is not trivial, the preliminary meshes still contained different types of artifacts that revealed themselves as drop-shaped notches, pockets, or tunnels (Fig. 3a). To bypass this discrepancy, we designed a morphological filtering pipeline (Fisher et al. 1997) in Houdini, which is based on topological error and artifact removal, hole filling, reprojection for detail recovery, convex-hull buffer re-meshing, and stencil Boolean Subtraction (Fig. 3b). A detailed description of this procedure is given by Domej and Pluta (2020).
The last step consisted of computing non-uniform rational basis-splines (NURBS) surfaces in Geomagic Design X, which envelop the individual caves; the latter were then exported to Standard-for-the-Exchange-of-Product-Data (STP) files to obtain self-contained and closed cave volumes (Fig. 4).
For clarification, it should be noted that the processing of the digital data started out with as many separate cave and sub-cave point clouds as laser scanner stations existed inside the caves. During data handling, those were gradually combined into the seven caves with separate entrances, as shown in Fig. 1a.
The second major part of digital data processing concerns the unified point clouds acquired by the drone and laser scanner from its three stations outside the caves. Here, the first step was to remove automatically some high-roughness features-such as, e.g., shrubs, artificial structures, machinery, or even people-in order to avoid the generation of artifact spikes, sharp or creased edges, self-intersections, and/or non-manifold geometries during the meshing process (Leong et al. 1996). Points of the concerned features were eliminated employing principal component analysis (PCA; Martinez and Kak 2001) at different scales to train a support vector machine (SVM; Cortes and Vapnik 1995) implemented in CloudCompare in its plug-in named CANUPO Suite (French: Caractérisation de Nuages de Points; Brodu and Lague 2012). The resulting cleaned point cloud was then meshed in MeshLab via the ball pivoting algorithm (Bernardini et al. 1999). Artificially created unconformities and elongated features could be removed manually from the mesh with the software Blender. Eventually, the mesh was exported to a Standard Triangle Language (STL; Fig. 5) file with GeoMagic Wrap to obtain a coherent surface representing the slope topography.
At this stage, the previously generated single STP files representing the caves were merged with the slope surface by Boolean Addition-i.e., transforming the selfcontained and closed cave volumes into open and concave appendices to the slope surface. After import to Geomagic Wrap, both types of objects, i.e., the slope surface and the envelopes of the caves-appear as triangulated surfaces ( Fig. 6a) with an excess volume of the caves emerging through the openings in the slope surface. The intersection between the slope surface and the cave envelopes was carried out by manually selecting the edges that both geometries share and deleting the excess volumes per cave (Fig. 6b).
As a final step, the combined computer-aided design (CAD) object was converted from a triangulated surface into a sum of NURBS surfaces via Edelsbrunner's (2003) method of "wrapping finite sets in space" in Geomagic Wrap and exported as a single solid object readable in GTS NX.

Assessment of Photographic Material
The Sabereebi Cave Monastery has been revisited regularly since 2011, and structural changes in its facades and interior and local landsliding activity are documented by photographic imagery. As the goal for this study was a numerical model depicting the physical state of the cliff at the moment of the data acquisition by the laser scanner and the drone, we focused on photograph series and separate pictures of the drone survey taken shortly before and during the field survey in 2019 to identify current damage and endangered features, which either recently had collapsed or are likely to collapse in the near future.
One of the three most prominent features is a collapsed pillar at the entrance to cave 2 of about one meter in height, which was still standing until November 2018 but not to be found again in 2019 ( Fig. 7a-b; NACHPG 2019). Furthermore, several fallen boulders are lying at the entrance to cave 3 with diameters between half and one meter (Fig. 8); all of them previously had been part of the floor of the small sub-cave above the large entrance to cave 3.
Not directly part of an excavated cave-yet endangering at least two-is the jointed rock spur next to the entrance to cave 6 and just above cave 7 (Fig. 9). Brownish weathering affects the bulge over roughly two meters in height, and  Closed mesh envelopes for all caves derived from cleaned and resampled point clouds. The first of the triplet tiles per cave shows the used point cloud (in blue) and zones of point lack that were closed using the Screened Poisson Method (Kazhdan and Hoppe 2013). The other two triplet tiles show the respective cave from in-and outside the cliff (pictures as MeshLab exports). The approximate numbers of vertices for the creation of mesh envelopes are given by the second triplet tiles (in thousands); d doors, e large entrances, t tunnels, and w windows. One cavern in cave 3 was too sparsely covered by the point cloud to create a volume; hence, the few existing points (in yellow) were discarded from the dataset slope-parallel fractures could reasonably let expect toppling to be imminent.

Laboratory Tests of Rock Samples
Two series of geotechnical laboratory tests were considered during the geotechnical property characterization of the exsitu rock samples of sand-and siltstone, as those two layers are of particular interest. From top to bottom, the sandstone layer hosts roughly two-thirds of the cave volumes, whereas the siltstone layer encloses the lower third (Fig. 1a).
First, Bergamini (2020) conducted classic rock-mechanical experiments such as laser-granulometric analyses ( Fig. 10a), mercury porosimetry tests (Abell et al. 1999; Fig. 10b), uniaxial compressive strength (UCS) tests ( Fig. 11) on samples in their original form and after nanosilica treatment, of which the latter could play a role during future reinforcement work, and micro-computer tomography (micro-CT; Fig. 12a-b) analyses. Moreover, X-ray diffraction (XRD; Table 5 in the appendix) analyses and X-ray fluorescence (XRF; Table 6 in the appendix) analyses shed light on geochemical composition and mineral bonding within the rock samples, which also could serve for future studies involving the aspect of water seepage and saturation. The here presented study, however, focuses on unsaturated conditions, and evaluated parameters are given for dry residual and peak conditions (Table 1). Having tested the siltstone samples (i.e., A, B, and E) and the sandstone samples (i.e., C, D, and F), it appeared that those initially identified as siltstone during the field survey in 2019 revealed themselves, though, as sandstone according to the granulometric analyses. Accordingly, we assumed geotechnical parameters of the sand-and siltstone layers as very similar, although visually, there remains a clear difference in the cliff lithology ( Fig. 1a).
The second geotechnical test series ( Fig. 11) consisted of three UCS tests on untreated (i.e., without nano-silica treatment) samples under dry conditions (A1, A2, and A3) and three Brazilian Tests under the same conditions (A1, A2, and A71).
To compromise between the results of Bergamini (2020) and those of this study, we adopted the values of 150 kPa and 42° for the cohesion and the friction angle for sandand siltstone in all models. The estimated tensile cut-off of 50 kPa reflects the overall high degree of rock fracturing. The specific weight for sand-and siltstone in all FEM models correspond to those of Bergamini (2020; Fig. 12a  from the geotechnical tests in this study (i.e., from UCS test samples A1, A2, and A3, and from Brazilian Test samples A1, A2, and A71; Fig. 11).
Since no samples were collected for the layers consisting of topsoil, conglomerate, and clay, suitable values for dry residual and peak conditions were estimated from literature in comparable geologic and geoclimatic It should be noted that most of the listed values per parameter have a range and also differ between peak and residual conditions. Due to the need for single values for the numerical models, and based on experience, we chose a representative value for each parameter per layer, which is indicated in italic in Table 1.

Models of the Sabereebi Cave Monastery
Synthesizing all geometric, geologic, and geotechnical information, we established several FEM model series in GTS NX to study the static behavior of the cliff under gravity only-i.e., in an unsaturated environment. Models either cover the entire cliff or separate caves as so-called "box models". As caves 5, 6, and 7 lie close to each other, those three are represented only by one box model (Fig. 13). All model series share an orientation, in which the y-axis is pointing northwards, rigid boundaries at the bottom, and lateral model borders, which are fixed horizontally (i.e., in x-and y-directions having a free z-component; Fig. 13).
Also, the lithologic sequence of topsoil, conglomerate, sandstone, siltstone, and clay from top to bottom is present in all models-although in different fractions of the considered volume. As very little structural geological data are available for the cliff, we solved several classic three-point problems (Davis et al. 2011) from coordinates available in the 3D Fig. 10 Grain size (a) and pore size (b) distributions of the rock samples collected during the field survey in 2019. The samples initially identified as siltstone (i.e., A, B, and E), though, plot to the grain sizes of sandstone. Pore diameters are measured only up to 100 µm due to an instrumental restriction Fig. 11 Results of the uniaxial compressive strength (UCS) tests (A1, A2, A3, and mean Aµ c of them) and the Brazilian Tests (A1, A2, A71, and mean Aµ t of them); results are corrected according to the specimen shape effect (Obert et al. 1960). The dashed blue circles (A*, B*, and C*) refer to three representative results obtained by Bergamini (2020). Cohesion and friction angle for sand-and siltstone in the models are taken as 150 kPa and 42° as a compromise between the results of Bergamini (2020) and those of this study. The estimated tensile cut-off of 50 kPa reflects the overall high degree of rock fracturing. All depicted samples have a diameter of about 3 cm

Fig. 12
Micro-computer tomography (micro-CT) analyses of sandstone (a) and siltstone (b) as slice intersections and horizontal cross-sections (pictures as Avizo exports). Porosities of the sandstone samples (i.e., C, D, and F) amount on average to 38%, those of the siltstone samples (i.e., A, B, and E) on average to 40% (Bergamini 2020); porosities from samples tested in this study (i.e., UCS test samples A1, A2, and A3, and Brazilian Test samples A1, A2, and A71) reach on average 43%. All porosity estimations rely on a rock density reference of 2700 kg/m 3 Table 1 Geotechnical parameters used for all models The properties of sand-and siltstone are derived from laboratory tests; those for topsoil, conglomerate, and clay are estimated from literature (Gasbarrone 2005;Margottini and Spizzichino 2020;Shafiei and Dusseault 2008). The notations "µ" and "e" stand for mean value and estimation, respectively; values originate from Bergamini (2020). Elasto-plastic materials account for all six parameters, whereas purely elastic materials account only for the Young's modulus, the Poisson ratio, and the specific weight geometry that correspond to bedding planes. On average, the azimuth and dip turned out to be 34° and 7°, with both values fitting well to the on-site measurements. As the lithological sequence appears almost perfectly parallel in reality, bedding planes within the models were assumed as such (Fig. 13).
In this study, we consider several distinct model series in order to account for different geomechanical scenarios and interests of investigation (Table 2). Differences between the model series relate to their extent and the assigned material behaviors. We distinguish five series in total, of which the second and the third series are assessed together due to their similarity. For simplicity, model series are named with their respective abbreviations hereafter; their particularities are introduced in the following sections, respectively.
All models are continuum models, as a comprehensive kinematic assessment of fracture sets just started recently at the site, and exhaustive joint data are expected for a followup study.
Common to all models of the entire cliff are the tetrahedral meshes with a mesh size of 0.5 m within the black-framed cave zones (Fig. 13) and the tetrahedral slope mesh with a mesh size of 10 m, which adjust towards the embedded finer cave zone meshes. The meshes of the separate box models are solely tetrahedral with a mesh size of 0.5 m due to meshability issues.
Properties given in Table 1 are attributed to the mesh elements -i.e., not to geometric layers-and, hence, bedding planes do not appear as clear cut in the models. Calculations rely on the fully implicit Newton-Raphson Algorithm (Ben-Israel 1966) with convergence criteria, as shown in Table 3. Fig. 13 Layouts for the entire cliff models and the box models per cave (picture as GTS NX export). Lithological properties are applied to different mesh sections (i.e., not to different geometries); hence, lithological units follow the parallel bedding planes compromising on tetrahedral mesh elements. All box models have tetrahedral meshes and a mesh size of 0.5 m; else, they respect the same modeling conditions as applied to the entire cliff models

MC-Models
The pure Mohr-Coulomb Models rely on purely elastoplastic material behavior and target the entire cliff along with the seven individual caves and sub-caves, intending to visualize large-and small-scale slope deformation in-and outside the caves. Keeping all involved materials under Mohr-Coulomb conditions, models focus on large-scale slope deformation and reveal the setting of a classic landslide mechanism. Figure 14a represents the entire cliff after the settlement under gravity. The equivalent plastic strain pattern shows that deformation is mainly constrained to zones  (Table 2) after the settlement under gravity (a) and after automatic strength reduction (b) of ~ 33% applying to cohesion and friction angle (pictures as GTS NX export). Results show the formation of a classic superficial landslide mecha-nism with a rotational component due to overburden predominantly within the clay layer. Maxima at model boundaries seem to be artifacts. Color coding for equivalent plastic strain is auto-scaled; the one for total translation is equally spaced beneath the convex topography part of the cliff hosting the caves and, hence, a result of overburden; moreover, a sub-surficial deformation tendency is characterized by slightly higher equivalent plastic strain. This overall landslide mechanism becomes more apparent after automatic strength reduction of about 33% applying to cohesion and friction angle simulating the degradation of the involved material in a uniform way (Fig. 14b); the tensile stress is not affected by the reduction, as it is often ignored for geomaterials (MIDAS 2021). Equivalent plastic strain patterns show higher absolute values and display a sub-superficial landslide mechanism with a slight rotational component starting at the siltstone-clay transition and being particularly pronounced beneath the entrance pedestals of caves 2, 3, and 4 and likewise around the rock spur between caves 6 and 7 (Fig. 9).
Similarly, the total translations emphasize a classic landslide mechanism; cross-sections through cave 2 show predominantly vertical settling under gravity that propagates into the more gently inclined cliff foot, whereas the total translation concentrates in that very part of the slope after automatic strength reduction.
It is interesting to note that-like the gravitational overburden effect-also the classic landslide mechanism seems to be constrained by the convex topography part of the cliff; zones just west and east of the cliff hosting the caves are less affected by both phenomena.
According to a previous field study (NACHPG 2019), failures caused by deeper structural mechanisms-such as classic rotational landsliding-are expected and partly proven by general slope-parallel fracture sets (Figs. 9, 15a-d). Causes might be, for instance, seasonal and humidity-sensitive clay swelling or pre-existing fracture sets.

MC-EL5-Models
The Mohr-Coulomb Models with an elastic layer 5 prohibit the formation of a landslide mechanism by imposing a purely elastic material behavior on the lowest lithological layer (i.e., on the clay layer) and, thus, virtually disabling deformation within that layer and allowing for deformation analysis on a smaller scale around and within the caves.
Very similar to this model series are the Mohr-Coulomb Models with an elastic slope ( Table 2); instead of defining only the clay layer as purely elastic, also the upper four layers (i.e., topsoil, conglomerate, sandstone, and siltstone) surrounding the irregular black-framed cave zones (Fig. 13) are kept elastic to focus on plastic deformation analysis around and within the caves. As results are very similar, both model series are addressed to as MC-EL5-Models hereafter.
Equivalent plastic strain patterns under gravity throughout the five box models (Figs. 16, 17) exhibit particularly vulnerable zones around pillars, thin walls, door frames, windows, niches, and arches. The thereby detected zones fit well to those identified as critical or prone to fail during the field surveys of the past years. Also, the pedestal of cave 6 enlists to the above-mentioned strain affected architectural features, representing the only large-scale floor zone experiencing equivalent plastic strain ( Fig. 17; second of the triplet tiles of caves 5, 6, and 7)-probably as a result of the distinct weakness of the rock spur (Fig. 9).
A confirmative picture is drawn by the stress distribution within the entire cliff in terms of tensile and compressive stress, respectively, as the smallest and greatest principal components of the stress tensor ( Fig. 18a-b). Tensile stress zones concentrate mainly above the conglomerate-siltstone transition where the cave roofs conclude the voids inside the cliff. It is worth noting that tensile stresses seem particularly pronounced above cave entrances except the one above cave 3, as several large boulders had fallen recently from that location, releasing the built-up stress (Fig. 8). Compressive stresses primarily affect cave walls, floors, and pillars within the sandand siltstone layers. As a particular feature emerges again the rock spur between caves 6 and 7, which experiences strong compressive stresses from its top-including the conglomerate layer-to the very bottom beneath cave 7 in the clay layer. This latter fact is noteworthy since only the representation of compressive stress gives information about the mechanical state of the slope foot beneath cave 7; in most other models, cave 7 is not well represented as it is located within the clay layer and, thus, mostly considered elastic ( Fig. 17; second of the triplet tiles of caves 5, 6, and 7)-i.e., not deformable. Tensile and compressive stresses obtained from the box models are shown in Figure 23 in the appendix.  (Table 2; pictures as GTS NX exports). p pillar. Color coding is autoscaled and different per cave Fig. 17 As the caption for Fig. 16, but for caves 3, 4, 5, 6, and 7

MC-EL5-mSRM-Models
The Mohr-Coulomb Models with an elastic layer 5 manually applying the Strength Reduction Method (SRM; Dawson et al. 1999;Griffiths and Lane 1999) likewise target the entire cliff along with the five box models. However, they aim to analyze deformation around small-scale features in and around the caves focusing on weathering effects. The manual reduction of strength parameters applies to cohesion and tensile stress in ten equally spaced steps reducing both by 70% of their initial values (Fig. 19), as Castellanza et al. (2018) and Ciantia et al. (2018) suggest after conducting experimental weathering tests on soft rock. The friction angle is assumed constant, as short-term weathering processes are unlikely to have a significant reducing impact on it (Ciantia et al. 2015a). Given the geologic and geoclimatic environment, we assume only short-term exposure to extreme weather events; long-term weathering effects-including permanent water saturation (Ciantia et al. 2015b)-are not considered in this model series.
Being a further development of the MC-EL5-Models, the MC-EL5-mSRM-Models share the same static state under  (Table 2) after the settlement under gravity (pictures as GTS NX exports) and show high tensile stresses around the cave roofs, whereas compressive stresses affect cave walls, floors, and pillars in particular. Color coding is auto-scaled Fig. 19 Values of cohesion and tensile stress in ten steps as used in the Mohr-Coulomb Models with an elastic layer 5 manually applying the Strength Reduction Method (SRM; Dawson et al. 1999;Griffiths and Lane 1999). SRF Strength Reduction Factor; it is, by definition, the ratio between the initial and the reduced value of an entity (MIDAS 2021) gravity; thus, the equivalent plastic strain pattern shown in Fig. 20a is representative for both model series. Similarly to Figs. 16 and 17, pillars, thin walls, door frames, windows, niches, and arches are identified as particularly vulnerable zones as a result of the concentration of tensile and compressive stresses governed by the initial parametric conditions (Table 1), defining a Mohr-Coulomb Envelope (Figs. 11,19 at SRF c = SRF σ = 1). Thereto we compare a representation of the entire cliff after a manual strength reduction of 70% (Fig. 19 at SRF c = SRF σ = 3.33, Fig. 20b). Absolute values of equivalent plastic strain increase around the mentioned vulnerable zones as a result of the shrinkage of the Mohr-Coulomb Envelope-i.e., as a combined consequence of the reduction of cohesion and tolerated stresses due to weathering effects and higher tensile and compressive stress concentrations in already weakened zones. Theoretically, reducing the friction angle could also alter and shrink the Mohr-Coulomb Envelope, but-as mentioned above-it is assumed constant in this model series.
We conclude, therefrom, that such "non-uniform inwards" weathering effects-as documented, e.g., by Ciantia et al. (2018) for similar cave structures-affect primarily free faces and underling subsurfaces within the sand-and siltstone layers in which the majority of cave voids and architectural features are located and which, for this reason, deserve special attention during preservation measurements.

EL-Models
The purely elastic box models are used to assess the stress regimes and safety factors of pillars (Figs. 16 and 17, Fig. 23 in the appendix, Table 2). The structure named "p2b" is connected to the cave wall through a thin rock volume but evaluated as if it was a self-standing pillar; hereafter, it is referred to as "semi-pillar". As the models are based on purely elastic material behavior for all lithologic layers (Fig. 13), this model series does not include automatic or manual strength reduction.
Safety factor evaluations are obtained for each pillar in two steps. First, average stresses per horizontal pillar crosssection (in kPa) are retrieved from the models at about 1 m above the surrounding cave floors. Second, we use Obert and Duvall's (1967) equation for a general safety factor of a pillar  (Table 2) after the settlement under gravity (a) and after manual strength reduction (b) of 70% applying to cohesion and tensile stress (pictures as GTS NX exports). Results show a strain pattern affecting exposed cave walls uniformly and, therefore, could be attributed to weathering horizons predominantly within sand-and siltstone. Maxima at model boundaries seem to be artifacts. Color coding is auto-scaled in which SF represents the safety factor, σ UCS the uniaxial compressive stress obtained via UCS tests (i.e., Aµ c of Fig. 11), and σ mod the average axial stress obtained via numerical modeling. The uniaxial compressive stress (σ UCS ) in Eq. 1 should be measured on drill-core rock samples with a diameter-to-height ratio of 1. As Fig. 11 shows, samples A1, A2, and A3 of the UCS tests slightly divert from this ratio and obtained uniaxial compressive stresses were corrected using the following relationship: 1960), in which R UCS represents the uniaxial compressive stress obtained via UCS tests from a sample with a ratio between its diameter (d) and height (h) other than 1 but in the allowed ratio range of 0.5 to 2. This concept of the safety factor evaluation, according to Obert and Duvall (1967), has proven effective in a comparable underground study involving pillars by Castellanza et al. (2018) and de Silva and Scotto di Santolo (2018). It should be noted that the concept is similar to-yet should not be confounded with-the Strength Reduction Factor (SRF; Fig. 19) used for manual strength reduction in the MC-EL5-mSRM-Models.
Confirming interpretations of deformation patterns in terms of equivalent plastic strain (Figs. 16, 17) and tensile and compressive stress (Fig. 23 in the appendix), we consider pillars 4 and 5 and the semi-pillar p2b as most unstable due to their fragile structure and low safety factors (Table 4). According to the rankings of horizontal cross-sectional area and safety factor, it appears, however, that both parameters influence each other only partially as follows: • area ranking: p5 < p2b < p4 < p2a < p1 < p3 • safety factor ranking: p2b < p5 < p4 < p3 < p1 < p2a Pillars 4 and 5, and the semi-pillar p2b, share the condition of being thinner compared to other architectonical features; moreover, they bear a relatively heavy overburden due to their rather isolated emplacement within a greater excavated space. In contrast, pillars 1 and 2a are thicker and-at the same time-more homogeneously surrounded by rock volumes of similar properties leading to a more stable stress distribution within the respective pillars. This multifactorial influence on safety factors is particularly well exemplified by pillar 3, having by far the greatest cross-sectional area, featuring the greatest extent of surrounding rock volume, but still ranging among the lower safety factors of all considered pillars.
It is essential to consider that absolute values of safety factors depend on the uniaxial compressive stress value (Fig. 11) used in Eq. 1 with smaller uniaxial compressive stresses entailing smaller safety factors and vice versa; the relative safety factor ranking-i.e., the ranking of failure potential of the individual pillars-remains the same. Moreover, absolute safety factors below 1 theoretically indicate the failures of the six considered pillars. These results are consistent with the respective equivalent plastic strain patterns (Figs. 16,17), but not with the fact that-in reality-all pillars are still standing. The discrepancy arises from two aspects: (i) models do not assume fractures within the pillars, and hence, do not account for stress redistribution due to fracture-driven load discharge; (ii) UCS tests were conducted on collected ex-situ rock samples, whose tolerated compressive stress values are not necessarily representative for the pillars.

Discussion, Recommendations, and Perspectives
Synthesizing the results of the field surveys of the past years and the different series of numerical models (Table 2), we have taken a closer look at specific aspects regarding the stability conditions of the Sabereebi Cave Monastery. Simulations with distinct objectives and scales shed light on the overall stability and failure potential of the cliff hosting the seven individual caves and sub-caves. The large-scale MC-Models revealed a classic landslide mechanism due to the overburden of the cliff coherent with the documented slope-parallel fracture sets across the cave walls. A point of criticism regarding this model series could be the assumed material behavior according to Mohr-Coulomb, which is one of the classic "first options" for general non-linear terrain models. A substantial flaw of Mohr-Coulomb Models is, however, that the natural nonlinear material response to applied stresses-such as, e.g., strain-softening or strain-hardening-is not approximated well by perfectly elasto-plastic material behavior. Better results could probably be obtained by using other material laws and criteria for rock mass failures like, for instance, the Hoek-Brown Failure Criterion (Hoek and Brown 1980). Table 4 Safety factors for pillars in caves 1, 2, 3, 4, and 5 (Fig. 16,  Fig. 17, Fig. 23 in the appendix) The structure named "p2b" is connected to the cave wall through a thin rock volume but evaluated as if it was a self-standing pillar-i.e., as a "semi-pillar" Relying, nevertheless, on the proven landslide mechanism ( Fig. 14a-b), countermeasures to instability within the more gently inclined cliff foot could be a terrace-styled combination of geotechnical installations such as anchoring and Nature-Based Solution (NBS) installations in the form of small rock and/or wood dams in the realms of landscape and cultural compatibility. A system of drainage conduits could prevent uncontrolled infiltration, in particular to the lowest lithological layer, which consists mainly of clay, making it susceptible to swelling phenomena . NBS could also be adapted for footpaths to ensure safe-yet inoffensive-access to the cliff. The small-scale MC-EL5-, MC-ELs-, and EL-Models focused on the stability and failure potential of architectural features around and within the caves and sub-caves, locating particularly vulnerable zones around pillars, thin walls, door frames, windows, niches, and arches (Figs. 16,17,18a,b,Fig. 23 in the appendix). Models show significant less failure potential in zones where a recent pillar failure ( Fig. 7a-b) and an entrance roof collapse (Fig. 8) took place; in contrast, we expect an imminent collapse of the rock spur between caves 6 and 7, including the pedestal of cave 6, as one of the next failure events. The safety factor analysis of five pillars and one semi-pillar ranked those in caves 2 (semi-pillar p2b), 4, and 5 as most likely to fail, whereas those in caves 1, 2 (pillar p2a), and 3 revealed themselves as slightly more stable. According to their priority, all of these architectural features could undergo appropriate supporting and reinforcement measurements. Examples could be strutting (Figs. 15d, 21), anchoring, as well as mechanical and/ or chemical enforcement such as nano-silica treatment; drillability conditions of the different lithologic layers should be evaluated via Pull-out Tests before installation.
Another obvious-yet not to be neglected-aspect of preservation should be the strict prohibition of vandalism in the form of scarring by visitors (Fig. 15c), as it weakens the rock structure and ultimately gives way to erosion in various forms.
The large-and small-scale MC-EL5-mSRM-Models point towards a probable strong dependency on stability and weathering. The model series already relying on uniform weathering horizons-i.e., assuming that material strengths decrease heterogeneously throughout the rock volume-suggest that zones of high failure potential focus around the caves and sub-caves within the sand-and siltstone layers (Fig. 20a-b). Better results could be obtained with numerical models employing weathering gradients (Fig. 22) which account for open-air exposure in-and outside the caves, free face area, rock type-dependent degradation patterns, eventual water infiltration paths, and the frequency and intensity of heavy precipitation. In this context, the aspect of chemical grain bond within the rock mass and, in a second stage, mineral leaching could be addressed (e.g., Ciantia et al. 2018). Furthermore, data from a recently initiated kinematic assessment of fractures, including on-site extensometers ( Fig. 15b; Frodella et al. 2021), could be included in improved numerical models, as fracture sets are reported not only to be slope-parallel but also to slope-perpendicular and randomly orientated (NACHPG 2019). In this context, the transformation of the continuum into a discontinuum model can be considered.

Conclusions
Similar to other artificial cave complexes in Georgia, the Sabereebi Cave Monastery, some 60 km southeast of the capital Tbilisi poses the challenge of preservation to geologists, engineers, and archaeologists. The cliff into which these Georgian-Orthodox caverns, chapels, and churches  Schematic representation of the influence of weathering, water infiltration, and erosion on the probable retrogressive failures above and beneath the cave entrances. Results from finite element method (FEM) modeling could be improved by including a weathering gradient rather than assuming heterogeneous weathering conditions throughout the rock volume were carved consists of a weak sedimentary rock sequence comprising claystone on the very bottom, overlain by exceptionally soft sand-and siltstone, marine conglomerates, and a thin layer of topsoil-all of which bear a considerable failure potential to erosion.
In this study, we present-in the first part-a strategy to process point cloud data from drone photogrammetry as well as from laser scanners acquired in-and outside the caves into high-resolution CAD objects that can be used for numerical modeling of critical failure zones ranging from macro-to micro-scale. In the second part, we explore a series of static elasto-plastic finite element stability models. Therefore, the study serves as a comprehensive 3D stability assessment of the Sabereebi Cave Monastery on the one hand; on the other hand, the established procedure should serve as a pilot scheme, which could be adapted to different sites in the future, combining non-invasive and relatively cost-efficient assessment methods, data processing and hazard estimation.
We discuss results from four distinct model series featuring different levels of detail, each of which focuses on specific geomechanical scenarios such as classic landsliding due to overburden, deformation of architectural features as a result of stress concentration, material response to weathering, and pillar failure due to vertical load. Generally, we conclude that particularly vulnerable zones are located around pillars, thin walls, door frames, windows, niches, and arches, and appropriate recommendations for preservation and countermeasures against further destruction are discussed alongside perspectives for future work. In addition to the above-mentioned model refinements concerning weathering gradients and the data incorporation of ongoing surveys, also dynamic models could be of particular interest as the region is characterized by strong seismicity.

Appendix
See Tables 5 and 6 and Fig. 23.   Funding The project was conducted within the framework of a "Multidisciplinary Survey and Monitoring of the Gareja Rock Cut Complex, the Monument of National Value" funded by the National Agency for Cultural Heritage Preservation of Georgia and managed by the Ilia State University, Georgia, together with the Italian Institute for Environmental Protection and Research, Italy, and the support of the University of Milano-Bicocca, Italy. Data used in this work was jointly acquired and, accordingly, belongs to all three institutions.

Conflicts of interest
The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.