Assessing Endokarst Potential in the Northern Sector of Santo Ant ó nio Plateau (Estremadura Limestone Massif, Central Portugal)

: Karst is a peculiar natural landscape arising from high rock solubility and well-developed underground solutional channel porosity. It is unique for its surface relief (exokarst) and subsurface drainage, including cave systems (endokarst). In Portugal, karst areas mainly consist of marginal or low-density territories with great fragility and vulnerability and great geo-environmental richness that merits better policies and practices regarding their geo-conservation. Endokarst potential assessments can provide decision-makers and local authorities insight into present and future territorial management and planning. In this context, the main objective of this study was to produce a cartographic model to identify areas with a greater probability of containing karstic caves—i.e., a greater endokarst potential—in the northern sector of the Santo Ant ó nio Plateau (Estremadura Limestone Massif, Central Portugal). Geological, topographic, hydrogeological, and land cover data were collected, processed, and integrated into a spatial database using a Geographic Information System. The locations of known cave entrances in the study area were also identiﬁed from local public institutions and speleological team records. Subsequently, four conditioning factors were extracted from the data: lithostratigraphic units, fracture density, relief energy, and land cover. Using a multi-criteria decision-making analysis, each previously chosen conditioning factor and its respective classes were weighted using an analytic hierarchy process. The locations of known cave entrances served to evaluate the cartographic model built, with results showing an agreement of 81.9%. This prototype of the endokarst potential map for the study area may be used for strategic and operational environmental planning (at least on a local scale) to assist decision-makers, competent authorities, and local speleological teams. Its application may promote a more accurate and thoughtful deﬁnition of areas to be investigated, substantially reducing the time and costs associated with ﬁeld prospecting.


Introduction
Karst landscapes are characterised by a near absence of permanent surface water courses due to a high degree of meteoric water infiltration, which represents a disruption, in terms of solution continuity, of the surface drainage network.Additionally, they are characterised by the presence of depressions near the surface (such as sinkholes), bare rock sculpted by the action of water dissolution (such as karrens and karren fields), karstic springs usually located on the edge of massifs, and an underground landscape.This last is denoted as endokarst, where caves occur and sometimes form part of a complex network of conduits traversed by water in generally rapid flows [1,2].Several conditioning factors influence karstification, and its correct identification/processing can allow a systematic assessment of the endokarst potential.According to Dimuccio [3] ' [. . .]In terms of the conditions of susceptibility to karstification of a carbonate lithic massif, these result from the capacity of the water system to erode (erosivity) and from the resistance of the lithic body in being eroded (erodibility), both differentiating in space and time [. ..]'.
The construction of predictive cartographic models is an attractive yet complex task.Several authors consider that cartographic models comprise a compilation of different types of data combined to build an explanatory model of the spatial distribution of a given phenomenon [4][5][6][7].Geographic Information Systems (GIS) include a set of essential tools that assist in the construction of these models, allowing the quick gathering of all necessary information into a spatial database, applying automatic/semi-automatic analytical procedures, and providing results such as maps, graphs, or vector/raster spatial data files [8][9][10][11].
Multi-criteria decision-making (MCDM) analyses [32][33][34] use several geostatistical techniques to define and map the underground karstification potential of specific territories.Taheri et al. [35] employed an analytic hierarchy process (AHP) to map the sinkhole occurrence susceptibility in a karst region.They used a magnitude-frequency analysis of the inventoried sinkholes to verify their model, although they needed more relevant geological data.Moradi et al. [36] used an AHP as an additional tool to explore the hydrogeology of a karst massif and classify areas according to a karstification potential assessment.These authors also used fuzzy logic with the same objective, concluding that the fuzzy logic showed better prediction accuracy than the AHP for their specific study area.Nevertheless, both methods help revise the weight of used parameters (i.e., factors) in other regions.Zaree et al. [37] evaluated water resources through potential recharge in karst units using the altitude, slope, lithology, infiltration, and soil weighting model (APLIS).They then applied an AHP together with a technique for order of preference by similarity to the ideal solution (TOPSIS) to modify the weight of the APLIS model, thus producing three final maps.When these cartographic results were matched and compared with the location of karstic springs and fractures, TOPSIS was selected as the best method.
In Portugal, a preliminary tentative study to model endokarst potential was carried out in the Arrábida Chain [38], where researchers extensively evaluated the parameters that affect karst landscape development.They used a model based on an interactive combinatory analysis of some dependent variables (lithology, fractures, intersections between fractures, slope, and terra rossa deposits), considering their relative contributions to endokarst development.The implementation of the model was based on an interactive software application that provided a set of maps indicating the likelihood of endokarst occurrence at any given point.A comparative analysis of the results and the location of two known caves delivered an opportunity to understand the parameters that contribute most to the formation of endokarst.
Considering the availability of detailed geological knowledge at the local scale and many known cave entrances, the present study is expected to allow the construction of a model that spatially represents the underground karstification potential.Using GIS potentialities, a prototype of the endokarst potential map for the study area was produced with the primary objective of enhancing and better directing the speleological prospection that should precede environmental impact studies and the diversified tasks related to cave protection/conservation [39][40][41], scientific research/divulgation [42][43][44], territorial planning, and management in karst areas [45][46][47][48].The focus of the present work is to contribute to a more efficient management of the karst territory under study.Compared to other studies [33][34][35][36][37][38], this more practical objective led us to create a cartographic model procedure that could readily apply to other karstified carbonate massifs.Because it is easier to apply, it is easier to generalise its use, understand its weaknesses, and introduce all the necessary corrections to improve its robustness and efficiency.
After identifying the techniques used in other works with similar themes and objectives, as well as the strengths and weaknesses of these works, in the present study, we proceeded to identify, analyse, and weigh the conditioning factors of underground karstification affecting Jurassic carbonate units in the northern sector of the Santo António Plateau, located in the most emblematic and suggestive karstified massif of Portugal-the Estremadura Limestone Massif (sensu Martins [49]).

Study Area
The Santo António Plateau, with a spatial extension of approximately 54 km 2 , corresponds to a small area on the north-western edge of the Estremadura Limestone Massif in Central Portugal (Figure 1) [50].This study area was selected since it is affected, as is the entire massif, by a well-developed karstification [17,49,[51][52][53][54].Additionally, the inventory, exploration, and scientific dissemination of karstic caves carried out by local speleological teams [55], the Natural Park of Serras de Aire e Candeeiros, and the National Institute for the Conservation of Nature and Forests were also considered.However, even though this area has been extensively studied by geologists [54,[56][57][58][59] and geographers [49,50,53], allowing us to obtain important data, it was not possible to extend the study area further as the spatial range of these scientific investigations largely coincide with the Santo António Plateau.
Geologically, the Estremadura Limestone Massif corresponds to an extensive and thick outcrop of Jurassic carbonate units that characterise the central sector of the Mesozoic Lusitanian Basin [60,61].From a morpho-structural perspective, this massif comprises a set of limestone mountain ranges and plateaus, measuring approximately 36 km from north to south, with a maximum width of 23 km [62].Martins [49] described its shape as an iron spear, pointing to SW, where the coincidence between lithology and hypsometry makes it peculiar.This massif corresponds to a large block that is slightly folded and raised owing to large fault systems.In particular, the Santo António Plateau has a flat surface to the southeast with an approximately triangular shape [18] and folds toward a syncline arrangement with a large radius of curvature [49].To the north and northeast, it limits with the 'Costa de Minde' and 'Costa de Alvados' cliffs and to the west with the 'Costa da Mendiga' cliff (coincident with major faults).To the south, it limits with the Alcobertas depression, and in the southeast with the Arrifes fault and its escarpment.The northern sector of the plateau is higher, mainly due to tectonic uplift along the main regional faults (Minde, Alvados, and Mendiga).Internally, this plateau is quite fractured, with a predominance of accidents in the northwest-southeast direction and igneous rocks injecting some of these structures.The highest altitudes were recorded near the top of the Alvados, Minde, and Mendiga tectonic accidents [18,50].
The Santo António Plateau is a well-studied area in terms of its geology, geomorphology, and speleology [17,18,49,51,52,[54][55][56][57]59], which facilitated the construction and legitimisation of the model.Mainly, we relied on the Middle Jurassic carbonate facies characterised by Azerêdo [56], as well as the cave inventory job and previous assessment of lithostratigraphic units in terms of susceptibility to karstification carried out by Crispim [17,54].In this plateau exists a profusion of surface and subsurface karst features, such as karren fields, dolines, uvalas, poljes, stephead valleys (amphitheatres) (Figure 2), rock shelters, and caves (Figure 3) [18,49,53].The Santo António Plateau is a well-studied area in terms of its geology, geomorphology, and speleology [17,18,49,51,52,[54][55][56][57]59], which facilitated the construction and legitimisation of the model.Mainly, we relied on the Middle Jurassic carbonate facies characterised by Azerêdo [56], as well as the cave inventory job and previous assessment of lithostratigraphic units in terms of susceptibility to karstification carried out by Crispim [17,54].In this plateau exists a profusion of surface and subsurface karst features, such as karren fields, dolines, uvalas, poljes, stephead valleys (amphitheatres) (Figure 2), rock shelters, and caves (Figure 3) [18,49,53].Martins [49] mentioned that the underground forms in the Estremadura Limestone Massif are more evolved than those from the surface and must, therefore, be older.Speleological exploration revealed that the local caves are connected to complex networks of underground conduits and extensive galleries.Some of these caves have often been indicated to be part of a much more extensive underground network [52,55], clearly showing the likeliness of a well-developed and mature endokarst in the region.In the words of Fleury [51] '[. ..]This massif is a true sponge, increasingly corroded as one descends in-depth [. ..]'.Martins [49] mentioned that the underground forms in the Estremadura Limesto Massif are more evolved than those from the surface and must, therefore, be older.Spe ological exploration revealed that the local caves are connected to complex networks underground conduits and extensive galleries.Some of these caves have often been in cated to be part of a much more extensive underground network [52,55], clearly showi the likeliness of a well-developed and mature endokarst in the region.In the words Fleury

Material and Methods
In the present study, specific conditioning factors of underground karstification affecting carbonate units in the northern sector of the Santo António Plateau were identified,

Material and Methods
In the present study, specific conditioning factors of underground karstification affecting carbonate units in the northern sector of the Santo António Plateau were identified, analysed, and weighted (Table 1).We built a cartographic model based on AHP [63][64][65].This MCDM method was chosen because of its intrinsic ability to approximate human perception, making it more friendly for decision-makers, as demonstrated by numerous prior investigations' results [32][33][34]66,67].
The model's predictive capacity was evaluated and verified based on known cave entrance locations.The degree of dependency (correlation) between these components was obtained by analysing the receiver operating characteristic (ROC) curve together with the associated metric corresponding to the area under the curve (AUC) [68,69].

Data Collection and Pre-Processing
The methodology adopted in this study is shown in Figure 4. We began by researching the sources of information that allowed us to select an area of study adequate for our objectives (e.g., with known georeferenced endokarst details).We then exhaustively collected all the information to build the model: geological (including lithology, with related faciological/stratonomic characteristics, strata geometry, and tectonic structures), topographic, hydrogeological, and land cover data.Four conditioning factors were extracted from these collected data: lithostratigraphic units, fracture density, relief energy, and land cover (Table 1).
included in the modelling process.However, we verified that several studies that intended to determine/map areas more susceptible to deep karstification selected the same conditioning factors as us [36,38,[71][72][73][74].In some studies with slightly different objectives, factors associated with precipitation or temperature were considered, as the meteoric water temperature influences carbonate rocks' dissolution [36,72].Considering the small dimensions of the study area (Figure 1), using these additional factors would be unwarranted because there would be no significant spatial differentiation [43,72].Regarding the chosen conditioning factors (Table 1), the lithostratigraphic units permit the evaluation of susceptibility to karstification of each unit based on lithology, facies identification (granulometry, texture, and qualitative porosity), stratonomic characteristics, and strata geometry.The stratigraphic information and some aspects related to the strata geometry (strike and dip) were obtained from the Geological Map of Portugal However, it is essential to highlight that other factors associated with the dynamics of underground water circulation and the chemical properties of these waters may influence the occurrence of karstic caves.For instance, chemically aggressive waters retained in confined aquifers can be crucial in developing caves (speleogenesis) [1,2,70].For the study area, the few available hydrochemical and underground water circulation data [17] did not allow us to construct an adequate cartography of these conditioning factors to be included in the modelling process.However, we verified that several studies that intended to determine/map areas more susceptible to deep karstification selected the same conditioning factors as us [36,38,[71][72][73][74].In some studies with slightly different objectives, factors associated with precipitation or temperature were considered, as the meteoric water temperature influences carbonate rocks' dissolution [36,72].Considering the small dimensions of the study area (Figure 1), using these additional factors would be unwarranted because there would be no significant spatial differentiation [43,72].
Regarding the chosen conditioning factors (Table 1), the lithostratigraphic units permit the evaluation of susceptibility to karstification of each unit based on lithology, facies identification (granulometry, texture, and qualitative porosity), stratonomic characteristics, and strata geometry.The stratigraphic information and some aspects related to the strata geometry (strike and dip) were obtained from the Geological Map of Portugal (Chart 27-A, Vila Nova de Ourém, at a 1:50,000 scale) and Crispim's work [17].Additionally, we obtained information concerning the characterisation of the carbonate facies from Azerêdo [56,57,75].Data on the qualitative porosity of several sample rocks related to facies identified and characterised by Azerêdo [56] were obtained from Inês's work [58].
Data related to the faults (major, hidden, and probable faults) were obtained from the Geological Map of Portugal (Chart 27-A Vila Nova de Ourém), whereas joints and lineaments were obtained from Carvalho's work [59].The calculation of fracture density, quantified as km of fracture per km 2 , was performed using the Line Density tool in ArcGIS (ESRI TM ) for two groups of fractures: the 'faults' and 'joints and lineaments' groups.
Relief energy was determined as the difference between the terrain elevation and local base level in each drainage sector [2,76].Elevation data were obtained from a digital elevation model constructed using topographic data at a 1:10,000 scale.The limits of the drainage sectors and the elevation of the local base level were obtained from Crispim [17] using the positions of the karstic springs associated with each drainage sector.These data were digitalised from maps (~1:50,000 scale) included in Crispim's work [17].
Land cover was obtained from the National Land Cover Map of 2018 (1:25,000 scale), and the most generic category of classes (Level 1) was used since they are easily comparable in describing the land cover types.
We also used information on the locations of known cave entrances inventoried by groups of speleologists over the last five decades.These caves allowed us to evaluate the performance of the cartographic model by analysing the spatial overlap between the obtained endokarst potential and the location of the known cave entrances.
To integrate cartographic information (with various scales and resolutions) into the modelling process, we built a set of raster thematic layers used as input variables with a homogeneous spatial resolution of 5 m pixels (see Table 1).

Model-Building Strategy
After collecting all the information, choosing the best analysis method was essential.As we had several influencing factors to consider with both quantitative and qualitative data, it was necessary to use a multi-criteria analysis [35][36][37][38][71][72][73][74]77,78].We integrated multi-criteria analysis with a cartographic predictive model and subdivided the study area into training (1/3 of the study area) and testing (2/3 of the study area) subsets of data.The small size of the study area and the existence of two distinct areas-one with a typical karst environment and many inventoried karstic caves and another with few karst characteristics and no inventoried karstic caves-conditioned the delimitation of the training and testing areas.The training area was delimited to have a high concentration of caves and simultaneously cover two existing environments (karst and non-karst).Splitting the data into training and testing areas gives an impression of the goodness-of-fit of the model and its ability to predict new data; therefore, it indicates its generality and transferability [79,80].The training area must represent a sample of the study area, including criteria-related data and data used for the assessment (inventoried caves) (Figure 4).
In the AHP, pairwise comparisons make it easier and more accurate to assign a preference between only two alternatives (factors/criteria) considering a given objective.The AHP uses a ratio scale from 1 to 9 associated with qualitative everyday life appreciation, where the operator makes a judgment, resulting in a quotient a/b that represents the dominance of a criterion over another (Table 2) [66,67,81].These pairwise comparisons are stored in a reciprocal matrix, where a normalised eigenvector (scale 0-1) represents the weight (Wi) of each factor for the objective (eigenvalue method) (Table 3).A consistency index is applied for priorities to make sense, which cannot exceed 10% [63,81].The AHP uses a logical procedure based on justifiable axioms, making it a robust method [67,81].It is scalable, with a hierarchical structure that can be easily adjusted to fit many-sized problems and is not data-intensive.Compared to some well-known MCDM (e.g., PROMETHEE, MAUT, TOPSIS, ELECTRE), with AHP, the decision-maker takes the lead in establishing preferences and assigns weights [32][33][34].

Criteria (Related to the Chosen
Conditioning Factors) However, in AHP, structuring occurs when criteria have many sub-criteria, and the decision-maker tends to give more weight than when they are less detailed.In our study, we implemented three criteria with five sub-criteria and one criterion (land cover) with seven classes (see Figure 5), minimising this problem.There are also some issues related to judgment scales.The linear 1-9 scale can represent a problem because of the lack of sensitivity in situations where weights are unequally dispersed; a high concentration in certain weights coexists with a high dispersion of others.Several alternative scales have been proposed to solve this problem, but choosing one is tricky because AHP is mainly related to subjective issues [66,82].We used a linear scale to consider the relatively homogenous nature of our judgment scales.However, the software application used, based on Microsoft ® Excel ® 2016 and developed by Goepel [83], allows employing other scales (e.g., logarithmic, square root, balanced-n).Rank reversal is one of the biggest criticisms of the AHP [32]; for example, if a criterion alternative is added or another is removed, the rank order of preferences can change.To avoid rank reversal (detectable in the aggregation of individual judgments), several authors have proposed the calculation of normalised weights using the geometric means (Formula (1)) instead of the eigenvalue method because the former indicates the central tendency and provides the same ranking [66].The Goepel's application implements the geometric mean method.Implementing the AHP allowed us to consider our chosen conditioning factors (criteria) with a pairwise comparison for all factor classes (sub-criteria) and between the factors for a common objective to assess potential underground karstification.The karstification susceptibility assessment of lithostratigraphic units, in agreement with several authors [2,70,84]; and references included therein], was based on the following assumptions: (1) we considered the purest carbonate rocks, with the highest percentage of calcite, to be more susceptible to karstification (e.g., limestones are more susceptible of karstification than dolostones); (2) we considered the influence of the insoluble constituents' present, e.g., clay minerals decreased the susceptibility to karstification and quartz grains increased rock porosity and the susceptibility to karstification; (3) as texture and granulometry are very dependent on the type of carbonate facies, we evaluated them together, and micritic and matrix-supported textures were considered more susceptible to karstification than sparitic and grain-supported ones; (4) the stratonomic data allowed an analysis of the lithostratigraphic units as a whole, highlighting other characteristics that enhance karstification but are not perceptible in the simple individual assessment of lithologies and facies-there was greater susceptibility to karstification in carbonate successions with relatively homogeneous beds in comparison with the more heterogeneous successions, which, together with the soluble strata, interstratify other relatively more insoluble layers (such as quartz sandstones, marls, and more or less carbonaceous mudstones); and (5) in terms of the carbonate bed geometry, we considered that the less inclined beds caused a slower circulation of water in depth, thus providing a higher rate of dissolution; even any stratification joints, more or less thick, could act as a barrier to vertical water circulation, thereby enabling a relatively more horizontal circulation that takes advantage of the stratification planes.

Designation
Around tectonic structures, a greater spatial density of the fractures, related to the frequency with which faults, joints and simple lineaments occur, and the presence of intersection zones between them, clearly indicate areas with greater capacity for water infiltration in-depth and the possibility to guide the development of caves [2,70,90].We considered that a greater or lesser propensity for water infiltration was associated with the fracture type.Slightly different criteria were used to calculate the density between the 'faults' and 'joints + lineaments' groups in a GIS: we applied a search radius of 200 m to the 'faults' group, assuming a greater concentration of water infiltration for the 'joints + lineaments' group, and we applied a survey radius of 600 m, taking a more diffuse water concentration.This exercise resulted in two raster maps used for the definitive calculation of fracture density using Formula (2) and the ArcGIS raster calculation tool.
Fracture density = 1 2 [(diaclases and lineaments) + (faults Fracture density (km/km 2 ) was ranked into five qualitative classes (from Very Low to Very High) using the ArcGIS natural-breaks classification method (Figure 5).
The potential energy available for the development of karstic caves has a component that represents the difference between the altitude of the terrain and that of the karstic spring(s), which is associated with the local base level in a specific groundwater circulation sector (i.e., drainage sector) of the carbonate massif [2,76].This indicator is relevant for assessing the importance of topography in underground karstification and is referred to as relief energy hereafter.Figure 5 shows the relief energy ranked into five classes obtained through the ArcGIS natural-breaks classification method.We assumed that the altimetric difference was a driving force (potential energy) in terms of topography.This force influences the hydrodynamic capacity of the water running through the interior of the rock mass, with an evident influence on the development of karstification.Thus, the relief energy represents the role of topography as a driving force in karstification [2,36,76].
The land cover type also influences karstification.Meteoric waters are acidified in areas covered by soil and vegetation, and their percolation is slower, enhancing carbonate rocks' dissolution.We considered that dissolution was enhanced in forest areas, pastures, and agricultural fields.In bush and uncovered areas, karstification dynamics are closer to those of bare karst, as the influence of these land cover types on chemical properties is reduced (although some influence exists in bush areas) [36,87,91].In artificialised territories, soil sealing does not enhance karstification; however, this land cover type can only influence the development of current karstification and not the endokarst forms developed over geological time without anthropogenic action (Figure 5).
The known cave entrances allowed us to verify the predictive capacity of the cartographic model built by correlating the endokarst potential results with the referred entrances' location and analysing the ROC curves [68,69,92,93].
In the first phase of AHP implementation, we weighted the classes of the four karstification factors using a pairwise comparison based on the assessment of karst territory experts-speleologists, geomorphologists, and geologists-after which we made pairwise comparisons between the factors.The quality of these comparisons was denoted by a consistency ratio, which was lower than 0.10 for all cases in this study, indicating satisfactory results.The relative weight of each factor (Wi) indicates the order of priority for each class associated with a factor and between factors (Tables 5-9); two classes may have the same weight when they are considered to have the same magnitude of importance, as is the case with the relief energy classes 'Very high' and 'High'.Expectedly, for lithostratigraphic units, fracture density, and relief energy, 'Very high' and 'High' classes amounted for more than 60% of the importance to underground karstification.For the land cover, the 'Forests' category represented ~30%, and the 'Pastures' category 22% (Tables 5-8).When comparing factors, the lithostratigraphic units had almost 50% importance, with fracture density 31%, relief energy 13%, and land cover residual importance 8% (Table 9).The geographic information layers (in raster format) corresponding to each karstification factor (Figure 5) were considered to model the endokarst potential, and after reclassification, their classes coincided with the normalised Wi.Using a GIS, we calculated the endokarst potential with the following formula: Endokarst potential = Wi Classes Litostratigraphic units × Wi Litostratigraphic units The adequacy of the model was verified through the spatial confrontation of the results and the known cave entrance.The predictive capacity of the model was verified based on the analysis of the ROC curves and the AUC, elaborated using the results of the binomial endokarst potential and the location of the cave entrances in both the training and testing subareas, as well as in the entire study area.In each ROC curve, the endokarst potential is represented in decreasing order on the abscissa axis.At the same time, a cumulative distribution function of the inventoried caves (with the data split into 20 classes) is represented on the ordinate axis.The variables have no dependence relationship if the 'curve' translates into a diagonal.A quick inflexion of the curve in the positive direction to that diagonal indicates that the model performs well; if the curve inflects in the negative direction, the model performs poorly and is considered unacceptable [68,69,92,93].The AUC value allows a quantitative assessment of the model's predictive capacity.The AUC values go from 0% to 100% (or 0 to 1), with values closer to 100% indicating a greater predictive capacity.The value of 50% graphically coincides with the diagonal line representing a casual predictive capacity.Values below 50% express a worse predictive capacity than a random model, and the corresponding model cannot be considered acceptable [68].To assess the predictive capacity of the cartographic model built, we used the values adopted in studies related to natural hazards (e.g., susceptibility to landslides), particularly those by Guzzeti et al. [93] and Oliveira [69].In these investigations, AUC values between 75% and 80% indicated an acceptable model, values between 80% and 90% stated a very good model, and values greater than 90% indicated an excellent model.
In general terms, we attempted to improve the spatial database whenever possible to construct the endokarst potential cartographic model in the study area.The components for assessing the susceptibility to karstification of the lithostratigraphic units were fine-tuned after more accurate literature research.Moreover, acquiring more data on the drainage sectors improved the relief energy calculation.This approach led us to obtain promising results and reduce the subjectivity inherent to the AHP.

Results and Discussion
From the analysis of the cartographic model developed for the endokarst potential of the study area, overlapped by the location of the entrances of the known caves (Figure 6), we verified that the weights attributed through the AHP are consistent since most (73%) of these entrances fall within the classes of "Very High" or "High" endokarst potential (Figure 7).Conversely, no cave entrances are in the areas classified as having 'Very Low' endokarst potential.In contrast to this promising result, a considerable area northwest of São Bento and other sectors distributed throughout the territory under investigation were also classified as 'Very High' or 'High' despite not having known cave entrances.If the present assessment is valid, we should consider that there may also be karstic caves in these sectors.We also observed that areas with higher elevation had greater endokarst potential in these areas, which may partly explain the distribution of known cave entrances.The used speleological data only referred to the locations of cave entrances that may belong to complex underground conduit systems.Thus, a more adequate analysis of the results could be achieved using topographic data that allowed mapping the development of these systems and comparing them with the produced model.Furthermore, the inventory of cave entrance locations concerns current knowledge (at least until February 2021); this distribution may be modified by discovering new caves, thus conditioning the model's spatial validity.Table 10 presents data on the number of caves in each class of the chosen conditioning factors and their density per km².Against expectations, we observed that the 'Moderate' and 'High' fracture density classes had more cave entrances than the 'Very High' class.Nevertheless, we verified a higher density of cave entrances in the 'High' and 'Very High' classes.When we analysed the spatial relationship between relief energy and the location of known cave entrances, we observed a more significant number of these in areas with greater potential energy relief, such as those in the 'High' and 'Very High' classes.It should also be noted that the highest density of cave entrances was higher for the 'Very High' class, leading us to conclude that there is a greater concentration of caves in sectors with higher relief energy.Finally, when considering land cover, we obtained different values than those expected.In the 'Bushes' class, a more significant number and greater density of cave entrances were observed than those in the 'Forests', 'Pastures', and 'Agriculture' classes.As there were apparent inconsistencies between the weighted classes for the fac- We observed a correlation between the lithostratigraphic units that were most susceptible to karstification and the cave entrances distribution, namely, in the micritic limestones of Serra de Aire (J 2 SA ) and the limestones of Chão das Pias (J 2 CP ), which represent 38.2% of the study area and contain approximately 63% of the known cave entrances.The marly limestones and marls of Fórnea (J 1−2 Fo ), marls and marly limestones of Zambujal (J 2 ZA ), and bioclastic limestones of Codaçal (J 2  Co ), which occupy 31.26% of the area, contain approximately 35% of the cave entrances.By analysing these data, we confirmed the propensity for karstification of units belonging to the Middle Jurassic, which were also classified as the most susceptible to karstification (Table 4).In addition to the results obtained from the modelling process, we compared our endokarst potential map with data from another study with a similar theme in the same area.As a result of this comparative exercise, we verified that the units categorised as 'Karstified Units' by Crispim [17,54] overlap areas in the higher endokarst potential classes identified in the present study.Additionally, the decrease in endokarst potential corresponded to a decrease in these 'Karstified Units'.Areas with 'Very High' endokarst potential overlap approximately 98% of the total area of 'Karstified Units', while 'High' endokarst potential areas overlap 'Karstified Units' at ~92%.Furthermore, 'Moderate' endokarst potential areas overlap 'Karstified Units' and 'Weakly Karstified Units' at ~93%, and 'Low' endokarst potential areas overlap 'Moderately Karstified Units' and 'Weakly Karstified Units' at ~70%.Finally, the 'Very low' endokarst potential class overlaps the Table 10 presents data on the number of caves in each class of the chosen conditioning factors and their density per km 2 .Against expectations, we observed that the 'Moderate' and 'High' fracture density classes had more cave entrances than the 'Very High' class.Nevertheless, we verified a higher density of cave entrances in the 'High' and 'Very High' classes.When we analysed the spatial relationship between relief energy and the location of known cave entrances, we observed a more significant number of these in areas with greater potential energy relief, such as those in the 'High' and 'Very High' classes.It should also be noted that the highest density of cave entrances was higher for the 'Very High' class, leading us to conclude that there is a greater concentration of caves in sectors with higher relief energy.Finally, when considering land cover, we obtained different values than those expected.In the 'Bushes' class, a more significant number and greater density of cave entrances were observed than those in the 'Forests', 'Pastures', and 'Agriculture' classes.As there were apparent inconsistencies between the weighted classes for the factors analysed and the expected cave entrance's locations, these results must be relativised owing to the low importance attributed to land cover and the good results obtained in the model assessment.
In addition to the results obtained from the modelling process, we compared our endokarst potential map with data from another study with a similar theme in the same area.As a result of this comparative exercise, we verified that the units categorised as 'Karstified Units' by Crispim [17,54] overlap areas in the higher endokarst potential classes identified in the present study.Additionally, the decrease in endokarst potential corresponded to a decrease in these 'Karstified Units'.Areas with 'Very High' endokarst potential overlap approximately 98% of the total area of 'Karstified Units', while 'High' endokarst potential areas overlap 'Karstified Units' at ~92%.Furthermore, 'Moderate' endokarst potential areas overlap 'Karstified Units' and 'Weakly Karstified Units' at ~93%, and 'Low' endokarst potential areas overlap 'Moderately Karstified Units' and 'Weakly Karstified Units' at ~70%.Finally, the 'Very low' endokarst potential class overlaps the units of 'Moderately Karstified Units', 'Weakly Karstified Units', 'Non-karstified Units', and 'Non-karstifiable Units' at ~52%.
As previously mentioned, the assessment/verification of the cartographic model was based on ROC curve analysis (Figure 8).Approximately 13% of the area classified with 'Very High' endokarst potential (87% to 100% in the abscissa axis) explains 70% of the area with cave entrances.In contrast, approximately 32% of the area with lower endokarst potential ('Low' and 'Very Low') explains only about 5% of the area with cave entrances.For the entire study area, the AUC value was 81.9%, thus considering the model as acceptable or with a good predictive capacity, regardless of whether other more robust and sophisticated analytical methodologies (such as artificial intelligence) may be used.A comparison of different methodologies is beyond the scope of this study.units of 'Moderately Karstified Units', 'Weakly Karstified Units', 'Non-karstified Units', and 'Non-karstifiable Units' at ~52%.As previously mentioned, the assessment/verification of the cartographic model was based on ROC curve analysis (Figure 8).Approximately 13% of the area classified with 'Very High' endokarst potential (87% to 100% in the abscissa axis) explains 70% of the area with cave entrances.In contrast, approximately 32% of the area with lower endokarst potential ('Low' and 'Very Low') explains only about 5% of the area with cave entrances.For the entire study area, the AUC value was 81.9%, thus considering the model as acceptable or with a good predictive capacity, regardless of whether other more robust and sophisticated analytical methodologies (such as artificial intelligence) may be used.A comparison of different methodologies is beyond the scope of this study.

Conclusions
The endokarst potential of a carbonate massif depends on several conditioning factors that contribute differently to the underground karstification process and, as a rule, almost always work interdependently.Details associated with the carbonate strata (i.e., faciological/stratonomic and geometry characteristics) can determine a greater or lesser endokarst potential.For instance, relatively purer, finer-textured carbonate rocks in thicker and fewer dip strata enhance underground karstification at the local scale.An increase in the entry of meteoric water into the carbonate massif is intrinsically linked to

Conclusions
The endokarst potential of a carbonate massif depends on several conditioning factors that contribute differently to the underground karstification process and, as a rule, almost always work interdependently.Details associated with the carbonate strata (i.e., faciological/stratonomic and geometry characteristics) can determine a greater or lesser endokarst potential.For instance, relatively purer, finer-textured carbonate rocks in thicker and fewer dip strata enhance underground karstification at the local scale.An increase in the entry of meteoric water into the carbonate massif is intrinsically linked to the existence of fissures, such as fractures (faults + joints) and bedding planes.The hydraulic gradient enhances the velocity of these waters due to the altimetric difference between the input (e.g., cave entrances) and output points (karstic springs), which the relief energy can describe.How meteoric waters enter the rocky massif and circulate is of particular importance.It is also essential to characterise the existing cover (detrital siliciclastic cover, soil, vegetation).Endokarst potential is associated with the greater or lesser development of subterranean forms (caves) that sometimes have a surface expression.
The geological, topographic, hydrogeological and land cover characteristics of the northern sector of the Santo António Plateau (Estremadura Limestone Massif, Central Portugal) proved decisive for assessing endokarst potential.At the same time, the location of the known cave entrances was vital for the assessment/verification of the proposed cartographic model.Applying the same modelling process in larger areas of the Estremadura Limestone Massif and other karstified massifs of Portugal is essential to legitimate its extrapolation.It is also clear that the obtained cartographic result highly depends on the quality/quantity of the data on which it is based.The cartographic modelling process used in this study represents a correct and easy-to-use methodology for defining sectors in a carbonate massif with environmental vulnerability associated with karst hydrogeology, with more significant potential for hazards, and that (eventually) can be used in speleo-tourism activities.Indeed, from a practical perspective, it is relatively easy to compile existing geological information in official geological maps (lithology, strata geometry, tectonic structures), elevation data, and land cover cartography to implement an AHP.We attempted to reduce the subjectivity inherent in AHP by applying improvements in treating the used spatial geodatabase.Hence, we consider disseminating the used methodology the best way to improve it, and comparing its results with those of other methods will be necessary.It should be noted that the results obtained in this study refer to a local scale.In a small study area, the spatial distribution of some conditioning factors can be quite different from that in large territories.Data representation at different scales associated with different spatial and thematic accuracies presents a challenge when integrating into a multi-criteria analysis.Knowledge of the endokarst itself is critical, but the available speleological information generally only represents the entrances of the caves and not their extension.
Finally, we highlight the importance of our results for obtaining knowledge on the endokarst and what it represents for anthropogenic activity.Human occupation of karst territories requires intelligent territory planning and management.Still, the weaknesses inherent to the karst environment are often not considered, especially those regarding water resources, natural hazards, and mineral resource extraction.Indeed, many examples in Portugal (but not only) have demonstrated that karst dynamics can put people and property at risk; furthermore, karst environments have a very peculiar landscape, both at the surface and underground, which can be used for geotourism, scientific and educational activities.
Informed Consent Statement: Not applicable.

Data Availability Statement:
Land cover data can be found at https://snig.dgterritorio.gov.pt/, the National Registry of Geographic Data, where open data can be downloaded (Direção-Geral do Território ® ).Other datasets used for the present research will be available on request.

Figure 1 .
Figure 1.Location map of the study area indicating terrain elevation, local administrative boundaries, karstic springs, and the entrance of known caves.

Figure 1 .
Figure 1.Location map of the study area indicating terrain elevation, local administrative boundaries, karstic springs, and the entrance of known caves.

Figure 2 .
Figure 2. Some examples of the most emblematic surface karst features in the study area.(A) St head Valley (amphitheatres); (B) Mira-Minde Polje partially flooded in the winter of 2022.

Figure 2 . 27 Figure 3 .
Figure 2. Some examples of the most emblematic surface karst features in the study area.(A) Stephead Valley (amphitheatres); (B) Mira-Minde Polje partially flooded in the winter of 2022.Sustainability 2023, 15, x FOR PEER REVIEW 6 of 27

Figure 3 .
Figure 3.Some examples of the subsurface karst features in the study area.(A) One entrance of 'The Four Mouths Cave'; (B) Temporary karstic spring.

Figure 4 .
Figure 4. Methodological flowchart adopted in this work to construct and evaluate the endokarst potential cartographic model of the study area.

Figure 4 .
Figure 4. Methodological flowchart adopted in this work to construct and evaluate the endokarst potential cartographic model of the study area.

Sustainability 2023 , 27 Figure 5 .
Figure 5. Chosen karstification conditioning factors for the northern sector of the Santo António Plateau.Figure 5. Chosen karstification conditioning factors for the northern sector of the Santo António Plateau.

Figure 5 .
Figure 5. Chosen karstification conditioning factors for the northern sector of the Santo António Plateau.Figure 5. Chosen karstification conditioning factors for the northern sector of the Santo António Plateau.

Sustainability 2023 , 27 Figure 6 .
Figure 6.Endokarst potential map and location of the known cave entrance on the northern sector of the Santo António Plateau.

Figure 6 .
Figure 6.Endokarst potential map and location of the known cave entrance on the northern sector of the Santo António Plateau.

Figure 7 .Table 10 .
Figure 7. Percentage of inventoried caves corresponding to each karstification potential class in the study area.Table 10.Number of caves in each class of the chosen condition factors and their density per km².

Figure 7 .
Figure 7. Percentage of inventoried caves corresponding to each karstification potential class in the study area.

Figure 8 .
Figure 8. ROC curve for the study area (more details in the text).

Figure 8 .
Figure 8. ROC curve for the study area (more details in the text).

Table 1 .
Geodatabase constructed using GIS tools to obtain the chosen conditioning factors for the cartographic model to identify the endokarst potential in the study area.

Table 4 .
Assessment of susceptibility to karstification of the lithostratigraphic units.Carbonate facies characterisation by Azerêdo

Table 8 .
Pairwise comparison for land cover types.

Table 10 .
Number of caves in each class of the chosen condition factors and their density per km 2 .