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

1. Introduction

The construction of predictive cartographic models is an attractive and complex task. It consists of a territorial compilation of information from various fields of knowledge to, through their combination, build an explanatory model of the spatial distribution of a given phenomenon (Tomlin, 1990, 2017; De Smith et al., 2018; Bolstad & Manson, 2022). Geographic Information Systems (GIS) include a set of essential tools to assist in the construction of these models, as they allow to gather all the necessary information in a spatial database quickly, apply automatic/semi-automatically analytical procedures, and make the results available in maps, graphs, or vector/raster spatial data files (Church, 2002; Longley et al., 2011; Huang, 2018; Gao, 2021).

In the field of karst geomorphology (Jennings, 1985; White, 1988; Cunha, 1993; Ford & Williams, 2007; Palmer, 2007; Dimuccio, 2017; De Waele & Gutiérrez, 2022; and references herein), the underground karstification features (caves) assumes relevance when related to the activities of the human being. The human occupation of karst regions must be carried out through territorial planning and management that considers the inherent weaknesses of this type of natural environment, with particular attention to the aspects related to water resources (availability and vulnerability), natural hazards (e.g., collapses of caves roof, karst flash floods, landslides) and extraction of mineral resources (e.g., stone) (White, 1988; Crispim, 1995; Rodrigues, 1998). However, the karst environment not only has constraints on human activities but also has a peculiar (sometimes spectacular) landscape to enjoy (Brilha, 2005; Cunha et al., 2014, 2018; Aleksandar, 2019; Khalaf, 2022; Valjavec et al., 2022). This underground landscape becomes even more interesting for geo-tourism purposes (e.g., speleo-tourism) (Rodrigues et al., 2017; Torres et al., 2018; Ballesteros, 2019; Chiarini, 2022), but also in the scientific and didactic-pedagogical scope (Hobléa, 2004; Leslie, 2011; Abdellah et al., 2022).

In the most general multi-criteria decision-making (MCDM) analysis (Velasquez & Hester, 2013; Zlaugotne et al., 2020; Vakilipour et al., 2021; among others), several geostatistical techniques were previously used for defining and mapping the underground karstification potential of specific territories. Among this, standout the work of Taheri et al. (2015), where the AHP was used to map the sinkhole occurrence susceptibility in a karst region, as well as they used a magnitude-frequency analysis of the inventoried sinkholes to validate the model built, although some limitations due to the lack of some relevant geological data. In Moradi et al. (2016), the AHP was used as an additional tool in the hydrogeological research of a karst massif, which results in identifying the areas according to the karstification potential assessment. These last authors also used fuzzy logic with the same objective and, in the end, they considered that the fuzzy logic model has best prediction accuracy than the AHP model (for their specific study area); however, both methods allow the revision of the weights of parameters (= factors) to be used in another region. Zare et al. (2019) evaluate the water resources through the potential recharge in karst units using the APLIS (altitude, slope, lithology, infiltration, and soils) weighting model, then applying the AHP and TOPSIS (a technique for order of preference by similarity to ideal solution) to modify the weight of the APLIS model. These last authors produced three final maps (APLIS, AHP and TOPSIS) that were matched and compared with the location of karst springs and fractures - the TOPSIS was selected as the best method. In Portugal, a preliminary tentative about the interactive prospection of endokarst was made in the Arrábida Chain (Vargas et al., 2003), where the authors evaluate the parameters that affect the development of karst erosion at depth through a model based on interactive combinatory analysis of the dependent variables (lithology, fractures, intersections between fractures, slope and terra rossa deposits) taking into account their relative contributions to the development of endokarst. The implementation of the model was based on an interactive software application which permits as output a set of maps that indicate the likelihood of the occurrence of endokarst at any given point. Comparative analysis of the results and the location of the structures of two known caves provide an opportunity to understand the parameters that most contribute to the formation of endokarst.

Considering what was said in the previous paragraphs, the main objective of the present work is to identify, analyse and weigh the conditioning factors of underground karstification that affects Middle 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, 1949). Like some other studies, we used MCDM analysis to evaluate the local underground karstification potential (= endokarst potential), with the following conditioning factors: lithology (with the faciological and stratonomic characteristics - lithostratigraphic units), geological structure (which includes fractures and geometry of carbonate rock layers), topography and soil cover. First, we tried to understand the influence of each of these factors on underground karstification, and with that information, we established a hierarchy of their relative importance. With very detailed geological knowledge at the local scale and many known cave entrances, the present study may allow the construction of a more reliable model for a spatial representation of the underground karstification potential. Using the GIS potentialities, a prototype of the endokarst potential map for the study area is produced with the primary objective of helping and better directing the speleological prosecutions that, as a rule, should precede the environmental impact studies and the diversified tasks related to territorial planning and management in karst areas.

2. Study area

The so-called Santo António plateau, with a spatial extension of approximately 54 km², corresponds to a small area on the northwest edge of the most emblematic and suggestive Estremadura Limestone Massif (Central Portugal – Figure 1) (Rodrigues, 2020). The choice of this study area was based on the fact that it is affected (as the whole massif) by a well-developed karst (Fleury, 1925; Martins, 1949; Thomas, 1985; Ferreira et al., 1988; Crispim, 1992); as well as considering the remarkable work done by local speleological teams, the Natural Park of Serras de Aire e Candeeiros (PNSAC) and the Institute for the Conservation of Nature and Forests, I.P. (ICNF) in the inventory, exploration, and scientific dissemination of karstic caves (Canais and Fernandes, 1999; ICNF, 2020; PNSAC, 2020).

Geologically, the Estremadura Limestone Massif correspond to an extensive and thick outcrop of Jurassic carbonate units that identify the central sector of the Mesozoic Lusitanian Basin (Rocha and Soares, 1984; Kulleberg et al., 2013).

From a morpho-structural point of view, this massif comprises a set of limestone mountain ranges and plateaus, measuring about 36 km from north to south and having a maximum width of 23 km (Daveau, 2004). Martins (1949) describes it as a shape of an iron spear pointed to SW, where the coincidence between lithology and hypsometry gives it individuality. It is a large block, slightly folded, raised due to large fault systems.

The Santo António Plateau represents a flat surface to the SE with an approximately triangular shape (Rodrigues, 1998, 2020) and waves with a tendency towards a syncline arrangement with a large radius of curvature (Martins, 1949). To the north and northeast, it is bounded by the two cliffs “Costa de Minde” and “Costa de Alvados”; in the west by the cliff “Costa da Mendiga” (coincident with major faults); in the south, it is bounded by the Alcobertas depression and in the southeast by the Arrifes fault and its escarpment. The northern sector of the plateau is higher, mainly due to the tectonic uplift along the main regional faults (Minde, Alvados and Mendiga). Internally, this plateau is quite fractured, with a predominance of accidents in the NW-SE direction; igneous rocks inject some of these structures. The highest altitudes are recorded near the top of the tectonic accidents of the Alvados, Minde and Mendiga (Rodrigues, 1998, 2020).

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

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

Santo António plateau is a well-studied area in terms of Geology, Geomorphology, and Speleology (Fleury, 1925; Martins, 1949; Thomas, 1985; Crispim, 1992, 1995; Rodrigues, 1998; Azerêdo, 1993, 1998; Canais & Fernandes, 1999; Carvalho, 2013). This fact facilitated the construction and legitimation of the model. We highlight the importance of characterising the Middle Jurassic carbonated facies by Azerêdo (1993), the caves inventory job and the previous evaluation of lithological units in terms of susceptibility of karstification carried out by Crispim (1992, 1995). In the study area 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) (Martins, 1949; Ferreira et al., 1988; Rodrigues, 1998).

Figure 2. A) Stephead Valley; B) Mira-Minde Polje partially flooded in the winter of 2022.

Figure 2. A) Stephead Valley; B) Mira-Minde Polje partially flooded in the winter of 2022.

Martins (1949) mentions that, in the Estremadura Limestone Massif, the underground forms are more evolved than the surface ones and, therefore, must be older. He even talks about juvenile karst in consideration of superficial modelling. Through speleological exploration, it was found that the caves are connected to complex networks of underground conduits and extensive galleries. Exploration descriptions often indicate that some caves are part of a much more extensive underground network (Thomas, 1985; Canais & Fernandes, 1999), which clearly shows the high potential for the existence of a well-developed and mature endokarst in the region. In the words of Fleury (1925), we give an account of this developed endokarst: “This massif is a true sponge, increasingly corroded as one descends in-depth […]”.

Figure 3. A) One entrance of “The Four Mouths Cave”; B) Temporary karst spring.

Figure 3. A) One entrance of “The Four Mouths Cave”; B) Temporary karst spring.

3. Material and methods

In the present work, some conditioning factors of underground karstification that affect carbonated units in the northern sector of the Santo António plateau are identified, analysed, and weight. We build a cartographic model based on the AHP (Saaty, 1977; Saaty, 1986; Saaty, 2009). This MCDM method was chosen because of its intrinsic ability to approximate human perception, thus also being more friendly for the decision-makers, as demonstrated by the excellent results obtained in many previous investigations (Ishizaka & Labib, 2009; Velasquez & Hester, 2013; Saaty & Ergu, 2015; Zlaugotne et al., 2020; Vakilipour et al., 2021).

Subsequently, the model’s predictive capacity is evaluated/verified through the location of the cave entrances known up to now. The degree of dependency (correlation) between these components was obtained by the analysis of the Receiver Operating Characteristic (ROC) curve, together with the associated metric corresponding to the area under the curve (AUC) (Sweets, 1988; Braga, 2001).

3.1. Data collection and pre-processing

The methodology adopted in the present work is described in the flowchart of Figure 4. We started by researching sources of information that allowed us to select an area of study adequate to the objectives (e.g., with known georeferenced endokarst details). Then we proceeded to an exhaustive collection of all the information that would serve to build the model: geological (including lithology, with the related faciological and stratonomic characteristics, as well as the strata geometry and tectonic structures), topographic, hydrogeological, and vegetation cover data. These factors are considered conditioning factors for the karstification of carbonate rocks. However, others influence the occurrence of caves and are associated with the dynamics of underground water circulation and the chemical properties of these waters. It is often referred to the great importance of chemically aggressive waters retained in confined aquifers to develop caves (speleogenesis; White, 1988; Klimchouk et al., 2000; Ford & Williams, 2007). The characteristics of the hydrochemical data and the underground water circulation (Crispim, 1995) did not allow us to produce adequate cartography for the modelling process. In addition, we verified that in several studies where it is intended to determine/mapping areas more susceptible to deep karstification, the used conditioning factors agree with our choices (Vargas et al., 2003; Paiva, 2014; Seif & Ebrahimi, 2014; Moradi et al., 2016; Nola & Bacellar, 2021). In some studies with slightly different objectives, factors associated with precipitation or temperature appear (in the sense that the temperature of meteoric waters influences the dissolution of carbonate rocks) (Seif & Ebrahimi, 2014; Moradi et al., 2016). Due to the small dimension of the study area, it would not make sense to use these factors because there would be no significant spatial differentiation.

The information collected that was integrated into the cartographic model (described in Table 1) is related to the following karstification conditioning factors: lithostratigraphic units, which permit to evaluate the susceptibility of karstification for each unit based on lithology, facies characteristics (granulometry, texture, and qualitative porosity), stratonomic character and strata geometry. The stratigraphic information and some aspects related to the strata geometry (strike and dip) come from the Geological Map of Portugal (Chart 27-A, Vila Nova de Ourém at 1:50 000 scale), as well as from the work of Crispim (1995). We also used Azerêdo (1993, 1998, 2007) information concerning the characterisation of the carbonated facies. At least we used data about the qualitative porosity of several sample rocks evaluated by Inês (2010) related to some facies identified/characterised by Azerêdo (1993). Fracture density was quantified in km of fracture per Km², and the data related to the faults (major faults, hidden fault, and probable fault) were from the Geological Map of Portugal (Chart 27-A Vila Nova de Ourém); joints and lineaments data from the work of Carvalho (2013). The calculation of the fracture density was performed using the Line Density tool from ArcGIS (ESRI^TM) for two groups of fractures: the “faults” group and the “joints and lineaments” group. Relief energy was determined as the difference between terrain elevation and local base level in each drainage sector (Salomon, 2000; Ford & Williams, 2007). The elevation data were obtained from the Digital Elevation Model (DEM) calculated with topographic data at a scale of 1:10 000. The limits of drainage sectors and the elevation of the local karst base level were obtained from Crispim (1995) using the position of the karst springs associated with each drainage sector. This data was digitalised from maps (~1:50 000 scale) that were part of Crispim´s work. Land cover categories were obtained from the National Land Cover Map of the year 2018 (1:25 000 scale). To employ the cartographic model, we use the most generic category of classes (level 1) because they are easily comparable and simply describe the type of land cover. We also used information about the location of known cave entrances that groups of speleologists have inventoried over the last five decades. Those caves allow us to evaluate the cartographic model by analysing spatial overlap between the endokarst potential and the location of the knows cave entrances.

Table 1. Geodatabase was constructed using GIS tools to obtain the chosen conditioning factors for elaborating the cartographic model to identify the endokarst potential in the study area. In bold are the derived raster thematic layers, used as input variables in spatial modelling, with a spatial resolution of 5-m pixels.

Classification	Data collection of the study area			Derived raster thematic layer
(spatial database)	GIS data type	Scale or resolution	Source or citation	(Conditioning factor) in a GIS
Geological data	Lithology (polygon-vector)A	1:50 000	IGM (1998)	Lithostratigraphic units
	Strata geometry (point-vector)	1:50 000	Carvalho (2013)
	Tectonic structures (line-vector)	1:50 000	Carvalho (2013)	Fracture density
Topographic data	Contours (line-vector)B	1:10 000	PMM (2020)	Relief energy
	Points elevation (point-vector)B	1:10 000	PMM (2020)
Hydrogeological data	Drainage sectors (polygon-vector)	1:50 000	Crispim (1995)
	Karst springs (point-vector)	1:50 000	Crispim (1995)
Vegetation cover data	Land cover (polygon-vector)C	1:25 000	DGT (2018)	Land cover
Speleological data	Karstic cave entrances (point-vector)D	1:100 000	ICNF (2021)	Location of known cave entrances
AVector data digitalised on Geological Map of Portugal 27-A Vila Nova de Ourém (courtesy of Porto de Mós Municipality - PMM). IGM = Instituto Geológico e Mineiro.
BVector data from approved cartography was used to generate a digital elevation model (DEM, 5-m pixel).
CAvaliable by Direção-Geral do Território (DGT) at: https://www.dgterritorio.gov.pt/Carta-de-Uso-e-Ocupacao-do-Solo-para-2018 (accessed on 31 May 2021).
DAvaliable by local speleological teams through the Instituto da Conservação da Natureza e das Florestas, IP de Portugal (ICNF).

Table 1. Geodatabase was constructed using GIS tools to obtain the chosen conditioning factors for elaborating the cartographic model to identify the endokarst potential in the study area. In bold are the derived raster thematic layers, used as input variables in spatial modelling, with a spatial resolution of 5-m pixels.

Classification	Data collection of the study area			Derived raster thematic layer
(spatial database)	GIS data type	Scale or resolution	Source or citation	(Conditioning factor) in a GIS
Geological data	Lithology (polygon-vector)A	1:50 000	IGM (1998)	Lithostratigraphic units
	Strata geometry (point-vector)	1:50 000	Carvalho (2013)
	Tectonic structures (line-vector)	1:50 000	Carvalho (2013)	Fracture density
Topographic data	Contours (line-vector)B	1:10 000	PMM (2020)	Relief energy
	Points elevation (point-vector)B	1:10 000	PMM (2020)
Hydrogeological data	Drainage sectors (polygon-vector)	1:50 000	Crispim (1995)
	Karst springs (point-vector)	1:50 000	Crispim (1995)
Vegetation cover data	Land cover (polygon-vector)C	1:25 000	DGT (2018)	Land cover
Speleological data	Karstic cave entrances (point-vector)D	1:100 000	ICNF (2021)	Location of known cave entrances
AVector data digitalised on Geological Map of Portugal 27-A Vila Nova de Ourém (courtesy of Porto de Mós Municipality - PMM). IGM = Instituto Geológico e Mineiro.
BVector data from approved cartography was used to generate a digital elevation model (DEM, 5-m pixel).
CAvaliable by Direção-Geral do Território (DGT) at: https://www.dgterritorio.gov.pt/Carta-de-Uso-e-Ocupacao-do-Solo-para-2018 (accessed on 31 May 2021).
DAvaliable by local speleological teams through the Instituto da Conservação da Natureza e das Florestas, IP de Portugal (ICNF).

3.2. Model-building strategy

With all the information collected, choosing the best method of analysis was essential. As we had several factors influencing karstification to consider, with quantitative and qualitative data, it was necessary to use a multicriteria analysis (e.g., Figueiredo, 2001; Vargas et al., 2003; Paiva, 2014; Ramos et al., 2014; Seif & Ebrahimi, 2014; Taheri et al., 2015; Moradi et al., 2016; Zaree et al.,2019; Nola & Bacellar, 2021). We idealized the use of multicriteria analysis integrated in a cartographic predictive model, therefore we subdivided the study area in a training subset of data (1/3 of the study areand in a testing subset of data (2/3 of the study area). Splitting the data into training and test areas gives the impression of the goodness-of-fit of the model and the ability to prednto new data, therefore, its generality and transferability (Fielding & Bell, 1997; Bennet & Ibrahim, 2014). The training area must represent a sample of the study area: criteria-related data and data used for the evaluation (inventoried caves) (Figure 4).

In AHP, using pairwise comparisons, it’s easier and more accurate to express one’s opinion between only two alternatives (factors/criteria) considering an objective. AHP uses a ratio scale from 1 to 9 associated with qualitative everyday life appreciation where the operator makes his judgement, which results in a quotient a/b that represents the dominance of a criterion over another (Ishizaka & Labib, 2009; Saaty & Ergu, 2015). Those pairwise comparisons are stored in a reciprocal matrix where a normalised eigenvector (scale 0-1) represents the weight of each factor for the objective (eigenvalue method). A consistency index is applied for priorities to make sense, which cannot exceed 10% (Saaty, 1977). The AHP uses a logical procedure based on justifiable axioms, making it a robust method (Saaty, 1977; Saaty & Ergu, 2015). It is scalable with a hierarchy structure that can easily adjust to fit many sized problems and is not data intensive. Comparatively, to some well-known MCDM (e.g., PROMETHEE, MAUT, TOPSIS, ELECTRE), the decision-maker takes the lead in making preferences, clearly assigned with weights (Velasquez & Hester, 2013; Zlaugotne et al., 2020; Vakilipour et al., 2021). 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 have three criteria with five sub-criteria and one criterion (land cover) with seven classes (see Figure 5), which minimise this problem. There are also some issues related to judgement scales; the linear scale 1-9 can represent a problem because lack of sensitivity in situations in which weights are unequally dispersed, that is, much concentration in certain weights coexisting with a high dispersion of others. Several alternative scales are proposed to solve this problem but it is difficult to choose one because AHP is mainly related to subjective issues (Hämäläinen & Salo, 1997; Ishizaka & Labib, 2009). In the face of the nature of our judgment scales (which are relativity homogenous), we use the linear scale. However, the software application we used (based on Microsoft Excel^® and developed by Goepel, 2018) allows using other scales (e.g., logarithmic, square root, balanced-n). The problem of rank reversal is one of the biggest criticisms of AHP (Velasquez & Hester, 2013); for example, if a criterion alternative is added or another is removed, the rank order of preferences can change. To avoid this rank reversal (detectable in the aggregation of individual judgments), several authors proposed the calculation of normalised weights by the geometric mean instead the eigenvalue method because the first indicates the central tendency and provides the same ranking (Ishizaka & Labib, 2009). The Goepel Excel application, used by us, implements the geometric mean method.

Figure 4. Methodological flowchart adopted in this work to construct and evaluate the cartographic model.

Figure 4. Methodological flowchart adopted in this work to construct and evaluate the cartographic model.

Implementing AHP allowed us to consider our factors (criteria) in a pairwise comparison for all classes of the factors (sub-criteria) and between the factors for a common objective - to assess the potential underground karstification. The assessment of karstification susceptibility of lithostratigraphic units, in agreement with several authors (Klimchouk et al., 2000; Ford & Williams, 2007; Barton, 2013; and references included therein), was based on the following assumptions: (1) regarding the type of rock, we considered the purest carbonate rocks, with the highest percentage of calcite, to be more susceptible of karstification. (e.g., limestone are more susceptible of karstification than dolostones); (2) type of constituents, such us the presence of insoluble ones, like clay minerals that decrease the susceptibility of karstification or quartz grains that can increase rock porosity and the susceptibility of karstification; (3) in relation to texture and granulometry, these are very dependent on the type of carbonate facies, hence they were evaluated together, micritic and matrix-supported textures are more susceptible of karstification than sparitic and grain-supported ones; (4) the collection of information on the stratonomic allowed an analysis of the lithostratigraphic units as a whole, highlighting other characteristics that enhance karstification that are not perceptible in the simple assessment of lithologies and facies, there is greater susceptibility of karstification in carbonated successions with relatively homogeneous beds in comparison with the more heterogeneous successions in which, together with the soluble strata, they interstratify other relatively more insoluble layers (such as quartz sandstone, marls and more or less carbonaceous mudstones); (5) in terms of the geometry of the carbonated beds, it is considered that the less dip soluble beds cause a slower circulation of water in depth, thus providing a higher rate of dissolution (even any stratification joints, more or less thick, can act as a barrier to vertical water circulation, thus enabling a relatively more horizontal circulation taking advantage of the stratification planes). There was a need to disaggregate the lithological information available for the study area, to retain, above all, the data related to the characteristics that influence the karstification in a more obvious way and by the predefined theoretical assumptions, supported by the literature of the specialty (Rauchspeciality, 1970; Waltham, 1981; James & Choquette, 1984; Klimchouk et al., 2000; Salomon, 2000; Ford & Williams, 2007; among others). The analysis of lithostratigraphic units focused on assessing their susceptibility to karstification, using descriptions and associations related to lithology; type of facies; granulometry; texture; stratonomic characteristics; primary (apparent) porosity; and bedding geometry. These components were individually evaluated on a quantitative scale (0 to 1) for greater or lesser susceptibility to karstification. Their sum indicates the susceptibility to karstification of each lithostratigraphic unit in a 0 to 4 quantitative scale and correspondent qualitative assessment for five classes: 0-1 (Very low); 1-1,5 (Low); 1,5-2,3 (Moderate);2,3-3 (High); 3-4 (Very high) (Table 2).

Table 2. Evaluation of susceptibility of karstification of the lithostratigraphic units. Carbonate facies characterisation by AZERÊDO (1993, 1998). Apparent porosity by INÊS (2010) or evaluated based on the type of lithology, granulometry and texture.

DesignationA	Lithology	SCB	Facies	Granulometry, texture and qualitative (apparent) porosityC	SC	StratonomyD	SC	Bed geometryE	SC	Σ - (qualitative assessment)F
Alluvium (a); detrital unit and terra rossa of Estremadura Limestone Massif (A)	Siliciclastic deposits, sometimes with a marly component	0	Pelitic and sandy facies	Pelites and sands with a generally clast-supported texture. Excellent porosity (evaluated based on the type of lithology, granulometry and texture).	0,9	Not show an apparent organisation in sedimentary beds but rather a massive structure filling the valley bottoms (alluvium) and some depressions and crevices of a karst nature (siliciclastic deposits with terra rossa).	0	-	0	0,9 (Very low)
Beds of Alcobaça (J3AI)	Marls, sometimes siltstones, limestones, and sandstones	0,2	Mudstones, silty-sandy carbonate clays and silty-clay sandstones	Pelites to matrix-supported sands containing various fossiliferous associations. Poor porosity (evaluated based on the type of lithology, granulometry and texture).	0,2	From base to ~57 m (thick beds of silty-sandy clays and silty-clay sandstones); 57 to 102 m (beds of sandstones (~30 m), sometimes fine or coarse, micaceous and silt-sandy clays, forming finer intercalations (1 to 5 m); 102 to 215 m (micaceous sandstones - ~113 m), from fine to very fine, grey and micaceous silt-sandy clays, grey green, with calcareous concretions, sometimes ferruginous); 215 to 269 m (silt-sandy carbonated clays - ~54 m), with fine calcareous intercalations); >269 m (micritic limestones).	0,2	Thickness about 150 to 200 m; average slope 25º (min. 15º and max. 38º).	0,3	0,9 (Very low)
Cabaços and Montejunto beds (J3CM)	Clay limestones and marls	0,4	Mudstones/wackstones related with rare packstones and grainstones.	Clay micritic limestones, limestones associated with marl, microsparitic limestones, bioclastic pelimicritic limestones. Less than 10% grain > 2mm; presence of carbonated microcrystalline matrix; non-contact grains (matrix-supported texture). Poor porosity (evaluated based on the type of lithology, granulometry and texture).	0,4	Base with 3 m (limestone at the base and yellowish marls with ferruginous concretions); 6 m (very bioclastic pelmicritic clayey limestone); 0.40 m (clay-limestone); 0.50 m (clay micritic limestone); 20 m (intraclastic micritic limestone); 2 m (a monogenic conglomerate of limestone matrix, compact)? (microsparite limestone, with peloids and intraclasts); 30 m (very bioclastic pelmicritic limestone).	0,3	Thickness about 65 m; average slope 34º (min. 8º and max. 85º).	0,2	1,3 (Very low)
Limestones of Moleanos (J2MI)	Limestones	0,5	Lithofacies 2 - rudstones; grainstones; Oolitic/bioclastic/ oncolytics/lithoclastics packstones	More than 10% grain > 2mm; Absence of carbonated crystalline matrix; grains in contact (grain-supported texture). Porosity evaluated by INÊS (2010): from poor to excellent.	0,7	Base 20-30 m (alternations of well-calibrated oolitic limestones and coarser calciclastics, sometimes with erosive base levels); 35-40 m (more compact, pelbiomicritic limestones); top; >100 m (succession becomes more clastic, with massive levels of calciclastic limestones).	0,6	Thickness is about 150 meters, probably 180-200 meters; average slope between 20 to 25º.	0,6	2,4 (High)
Micritic limestone of Serra de Aire (J2SA)	Limestones	0,8	Lithofacies 6 – Mudstones and wackstones oncolíticos with fenestrae and laminations (idem). Lithofacies 7 - Floastones, wackstones and mudstones with algal/oncoid nodules and rusting.	Dolomitic levels; micritic and dolomicritic limestones (compact or laminar). Less than 10% grain > 2mm; presence of carbonated microcrystalline matrix; non-contact grains; more than 10% grain > 2mm and non-contact grain. Poor for Lithofacies 6 and Poor to Reasonable for Lithofacies 7.	0,8	Base 50 m (cyclical sequences of micritic dolomitic limestones and limestones, fenestrated decimetric layers; laminar dolomicrites); <150 m (compact micritic, fenestrated or oncolitic limestones, with ferruginous tinges, in 40-50 cm to metric layers); >150 m (decreased fenestrated and oncosparitic limestones, becoming fossiliferous micritic limestones; biomicritic or pelimicritic limestones).	0,9	Thickness from 350 to 400 m; average slope 14º (min. 2º and max. 70º).	0,8	3,3 (Very high)
Bioclastic limestone of Codaçal (J2Co)	Limestones	0,7	Lithofacies 1 – Grainstones oolitic and bio-intraclastic with oblique bedding (idem). Lithofacies 2 - Rudstones, grainstones and bioclastic/oncolytic/lithoclastics packstones.	Biolclastic and oobioclastic and sporadically dolomitized limestones. More than 10% grain > 2mm; Absence of carbonated crystalline matrix; grains in contact. Poor to fair for Lithofaces 1 and poor to excellent for Lithofaces 2.	0,8	Base 8-10 m (well calibrated fine oolitic limestones with small-scale oblique lamination/stratification (0.5-2cm and 16-28º), millimetric hardground 5 m from the base); 5 m (fragmented and/or bioperforated micritic coatings at the base of laminae of oblique stratified bundles of oolitic limestones (10-13 cm and 18-20 cm); 10 m (bioclastic and oolitic limestones with intraclasts, oblique stratification with 12º and 22º); 5-6 m (bioclastic and oolitic limestones with average thickness of oblique bedding bundles: 30-40 cm to 70-100 cm).	0,5	Average thickness is about 50 to 60 m, with a tendency to increase to 70-80 m; average slope 8º (min. 5º and max. 10º).	0,9	2,9 (High)
Limestones of Chão das Pias (J2CP)	Slightly clayey or marly limestone, limestone, dolomitic limestone	0,7	Lithofacies 9a – Mudstones, wackstones and bioclastic packstones (compact limestone) (ibidem).	Less than 10% grain > 2mm; Presence of carbonated microcrystalline matrix; non-contact grains; Absence of carbonated microcrystalline matrix and grains in contact. Poor porosity evaluated.	0,8	First 15 m (slightly clayey or marly limestone in decimeter benches with siliceous nodules); At 40 m from the top the nodules become larger (botryoidal). Succession characterized by the alternation of micritic and calciclastic limestones.	0,7	Thickness about 50-60 m, reaching, however, values >80 m? Average slope 9º (min. 5º and max. 15º).	0,9	3,1 (Very high)
Marls and marly limestones of Zambujal (J2ZA)	Marls, marly limestones, clayey limestones, limestones.	0,4	Lithofacies 9a – Mudstones, wackstones and bioclastic packstones (compact limestones). Lithofacies 9b – Mudstones, wackstones and bioclastic packstones (limestones, marl-clay limestones, and marls).	Less than 10% grain > 2mm; Presence of carbonated microcrystalline matrix. Non-contact grains; absence of carbonated microcrystalline matrix and grains in contact. Poor porosity.	0,8	Rhythmic alternation of marls, marly limestones, and clayey limestones, in almost always thin layers. The succession becomes increasingly thick (decimetric to metric layers) and calcareous from the bottom to the top until the marly levels disappear. It appears significantly fractured.	0,7	Thickness about 220-250 m; average slope 14º (min. 4º and max. 34º).	0,8	2,7 (High)
Marl limestones and marls of Fórnea (J1-2Fo)	Marls and marly limestones	0,4	Grumose; wackstones; biomicrites to biosparites/grainstone; packstones to grainstones.	Less than 10% grain > 2mm/More than 10% grain > 2mm; Presence of carbonated microcrystalline matrix; non-contact grains/absence of carbonated microcrystalline matrix and grains in contact. Poor to reasonable porosity (evaluation by type of lithology and texture).	0,7	Succession dominated by thin to medium layers, centimeter to decimeter, sometimes without rhythmic organization. At 80 m from the top occurrence of biostromal bodies with metric thickness. The upper 50 m are dominated by micritic limestones.	0,7	Maximum thickness with about 220-250 m; average slope 33º (min. 19º and max. 58º).	0,5	2,3 (Moderate)
Beds of Coimbra (J1Co)	Dolomites	0,5	Wackstones to grainstones.	Less than 10% grain > 2mm/More than 10% grain > 2mm; Presence of carbonated microcrystalline matrix; non-contact grains/absence of carbonated microcrystalline matrix and grains in contact. Poor to reasonable porosity (evaluation by type of lithology and texture).	0,7	Cross-bedding and dolomites with parallel or wavy lamination, interstratified with pellets.	0,4	Thickness about 60 m (beds with vertical or slightly inverted slopes).	0,1	1,7 (Moderate)
Platelet dolomites (J1pi)	Dolomitic limestones	0,7	Mudstone	Micritic dolomitic limestone. Poor porosity (evaluated based on the type of lithology).	0,2	Layers with centimeter to decimeter thickness.	0,1	Thickness about 30-40 m.	0,2	1,2 (Low)
Marls of Dagorda (J1Da)	Sandy loams, gypsum and saliferous clays, Intercalations of dolomitic limestones	0,1	Sandy, pelitic and Mudstones facies for the carbonate ones.	Poor porosity (evaluated based on the type of lithology).	0,1	250-320 m (dolomitic member, essentially dolomitic or margo-dolomitic with red and/or greyish pelites and evaporites); 60-850 m (saliferous/dolomitic member, predominantly dolomitic and/or limestone and marl rich in evaporites - anhydrite and halite); 290-800 m (saliferous member, characterised by an accentuated domain of halite, sometimes interstratified with dolomitic marls and/or marly pelites and anhydrite); From about 1000 m (occurrence of “Dolomites in platelets”); 1600-1200 m (evaporitic salts with intercalations of evaporitic syngenetic dolomite and gypsum); ~3000-1600 m (thick saliferous series with frequent clayey intercalations, also saliferous, containing anhydrite inclusions);	0,1	Formation subjacent Jurassic limestones with significant thickness, >3000 m, according to sounding “São Mamede 1”.	0,4	0,7 (Very low)
Eruptive rocks	Dolerite	0	-	Poor porosity (evaluated based on the type of lithology).	0	Associated with fractures or discontinuities.	0,8	-	0	0,8 (Very low)

^A Information collected from MANUPPELLA et al. (2000) and KULLEBERG et al. (2006, 2013). For Middle Jurassic units (Aalenian to Batonian), the facies described were proposed by AZERÊDO (1993, 1998). ^B Susceptibility to karsification, weighting 0 to 1. ^C Granulometry and texture according to the classification of DUNHAM (1962) and EMBRY & KLOVAN (1971) in SALOMON, 2000. Evaluation for primary porosity based on the classification of qualitative porosity adopted by INÊS (2010), from the work of AHR (2008): poor (0-5 %); reasonable (5-10%); good (10-15%); good/excellent (15-20%) and excellent (> 20%). For lithologies that were not evaluated by INES (2010), an evaluation was made based on the type of lithology and/or texture of the carbonate facies. ^D Information collected from MANUPPELLA et al. (2000). ^E Information collected in CRISPIM, 1995; Geological Map of Portugal 27-A Vila Nova de Ourém, 1998; MANUPPELLA et al., 2000. ^F Values from 0 to 4, with the following weighting: 0-1 (Very low); 1-1.5 (Low); 1.5-2.3 (Moderate); 2.3-3 (High); 3-4 (Very high).

Table 2. Evaluation of susceptibility of karstification of the lithostratigraphic units. Carbonate facies characterisation by AZERÊDO (1993, 1998). Apparent porosity by INÊS (2010) or evaluated based on the type of lithology, granulometry and texture.

DesignationA	Lithology	SCB	Facies	Granulometry, texture and qualitative (apparent) porosityC	SC	StratonomyD	SC	Bed geometryE	SC	Σ - (qualitative assessment)F
Alluvium (a); detrital unit and terra rossa of Estremadura Limestone Massif (A)	Siliciclastic deposits, sometimes with a marly component	0	Pelitic and sandy facies	Pelites and sands with a generally clast-supported texture. Excellent porosity (evaluated based on the type of lithology, granulometry and texture).	0,9	Not show an apparent organisation in sedimentary beds but rather a massive structure filling the valley bottoms (alluvium) and some depressions and crevices of a karst nature (siliciclastic deposits with terra rossa).	0	-	0	0,9 (Very low)
Beds of Alcobaça (J3AI)	Marls, sometimes siltstones, limestones, and sandstones	0,2	Mudstones, silty-sandy carbonate clays and silty-clay sandstones	Pelites to matrix-supported sands containing various fossiliferous associations. Poor porosity (evaluated based on the type of lithology, granulometry and texture).	0,2	From base to ~57 m (thick beds of silty-sandy clays and silty-clay sandstones); 57 to 102 m (beds of sandstones (~30 m), sometimes fine or coarse, micaceous and silt-sandy clays, forming finer intercalations (1 to 5 m); 102 to 215 m (micaceous sandstones - ~113 m), from fine to very fine, grey and micaceous silt-sandy clays, grey green, with calcareous concretions, sometimes ferruginous); 215 to 269 m (silt-sandy carbonated clays - ~54 m), with fine calcareous intercalations); >269 m (micritic limestones).	0,2	Thickness about 150 to 200 m; average slope 25º (min. 15º and max. 38º).	0,3	0,9 (Very low)
Cabaços and Montejunto beds (J3CM)	Clay limestones and marls	0,4	Mudstones/wackstones related with rare packstones and grainstones.	Clay micritic limestones, limestones associated with marl, microsparitic limestones, bioclastic pelimicritic limestones. Less than 10% grain > 2mm; presence of carbonated microcrystalline matrix; non-contact grains (matrix-supported texture). Poor porosity (evaluated based on the type of lithology, granulometry and texture).	0,4	Base with 3 m (limestone at the base and yellowish marls with ferruginous concretions); 6 m (very bioclastic pelmicritic clayey limestone); 0.40 m (clay-limestone); 0.50 m (clay micritic limestone); 20 m (intraclastic micritic limestone); 2 m (a monogenic conglomerate of limestone matrix, compact)? (microsparite limestone, with peloids and intraclasts); 30 m (very bioclastic pelmicritic limestone).	0,3	Thickness about 65 m; average slope 34º (min. 8º and max. 85º).	0,2	1,3 (Very low)
Limestones of Moleanos (J2MI)	Limestones	0,5	Lithofacies 2 - rudstones; grainstones; Oolitic/bioclastic/ oncolytics/lithoclastics packstones	More than 10% grain > 2mm; Absence of carbonated crystalline matrix; grains in contact (grain-supported texture). Porosity evaluated by INÊS (2010): from poor to excellent.	0,7	Base 20-30 m (alternations of well-calibrated oolitic limestones and coarser calciclastics, sometimes with erosive base levels); 35-40 m (more compact, pelbiomicritic limestones); top; >100 m (succession becomes more clastic, with massive levels of calciclastic limestones).	0,6	Thickness is about 150 meters, probably 180-200 meters; average slope between 20 to 25º.	0,6	2,4 (High)
Micritic limestone of Serra de Aire (J2SA)	Limestones	0,8	Lithofacies 6 – Mudstones and wackstones oncolíticos with fenestrae and laminations (idem). Lithofacies 7 - Floastones, wackstones and mudstones with algal/oncoid nodules and rusting.	Dolomitic levels; micritic and dolomicritic limestones (compact or laminar). Less than 10% grain > 2mm; presence of carbonated microcrystalline matrix; non-contact grains; more than 10% grain > 2mm and non-contact grain. Poor for Lithofacies 6 and Poor to Reasonable for Lithofacies 7.	0,8	Base 50 m (cyclical sequences of micritic dolomitic limestones and limestones, fenestrated decimetric layers; laminar dolomicrites); <150 m (compact micritic, fenestrated or oncolitic limestones, with ferruginous tinges, in 40-50 cm to metric layers); >150 m (decreased fenestrated and oncosparitic limestones, becoming fossiliferous micritic limestones; biomicritic or pelimicritic limestones).	0,9	Thickness from 350 to 400 m; average slope 14º (min. 2º and max. 70º).	0,8	3,3 (Very high)
Bioclastic limestone of Codaçal (J2Co)	Limestones	0,7	Lithofacies 1 – Grainstones oolitic and bio-intraclastic with oblique bedding (idem). Lithofacies 2 - Rudstones, grainstones and bioclastic/oncolytic/lithoclastics packstones.	Biolclastic and oobioclastic and sporadically dolomitized limestones. More than 10% grain > 2mm; Absence of carbonated crystalline matrix; grains in contact. Poor to fair for Lithofaces 1 and poor to excellent for Lithofaces 2.	0,8	Base 8-10 m (well calibrated fine oolitic limestones with small-scale oblique lamination/stratification (0.5-2cm and 16-28º), millimetric hardground 5 m from the base); 5 m (fragmented and/or bioperforated micritic coatings at the base of laminae of oblique stratified bundles of oolitic limestones (10-13 cm and 18-20 cm); 10 m (bioclastic and oolitic limestones with intraclasts, oblique stratification with 12º and 22º); 5-6 m (bioclastic and oolitic limestones with average thickness of oblique bedding bundles: 30-40 cm to 70-100 cm).	0,5	Average thickness is about 50 to 60 m, with a tendency to increase to 70-80 m; average slope 8º (min. 5º and max. 10º).	0,9	2,9 (High)
Limestones of Chão das Pias (J2CP)	Slightly clayey or marly limestone, limestone, dolomitic limestone	0,7	Lithofacies 9a – Mudstones, wackstones and bioclastic packstones (compact limestone) (ibidem).	Less than 10% grain > 2mm; Presence of carbonated microcrystalline matrix; non-contact grains; Absence of carbonated microcrystalline matrix and grains in contact. Poor porosity evaluated.	0,8	First 15 m (slightly clayey or marly limestone in decimeter benches with siliceous nodules); At 40 m from the top the nodules become larger (botryoidal). Succession characterized by the alternation of micritic and calciclastic limestones.	0,7	Thickness about 50-60 m, reaching, however, values >80 m? Average slope 9º (min. 5º and max. 15º).	0,9	3,1 (Very high)
Marls and marly limestones of Zambujal (J2ZA)	Marls, marly limestones, clayey limestones, limestones.	0,4	Lithofacies 9a – Mudstones, wackstones and bioclastic packstones (compact limestones). Lithofacies 9b – Mudstones, wackstones and bioclastic packstones (limestones, marl-clay limestones, and marls).	Less than 10% grain > 2mm; Presence of carbonated microcrystalline matrix. Non-contact grains; absence of carbonated microcrystalline matrix and grains in contact. Poor porosity.	0,8	Rhythmic alternation of marls, marly limestones, and clayey limestones, in almost always thin layers. The succession becomes increasingly thick (decimetric to metric layers) and calcareous from the bottom to the top until the marly levels disappear. It appears significantly fractured.	0,7	Thickness about 220-250 m; average slope 14º (min. 4º and max. 34º).	0,8	2,7 (High)
Marl limestones and marls of Fórnea (J1-2Fo)	Marls and marly limestones	0,4	Grumose; wackstones; biomicrites to biosparites/grainstone; packstones to grainstones.	Less than 10% grain > 2mm/More than 10% grain > 2mm; Presence of carbonated microcrystalline matrix; non-contact grains/absence of carbonated microcrystalline matrix and grains in contact. Poor to reasonable porosity (evaluation by type of lithology and texture).	0,7	Succession dominated by thin to medium layers, centimeter to decimeter, sometimes without rhythmic organization. At 80 m from the top occurrence of biostromal bodies with metric thickness. The upper 50 m are dominated by micritic limestones.	0,7	Maximum thickness with about 220-250 m; average slope 33º (min. 19º and max. 58º).	0,5	2,3 (Moderate)
Beds of Coimbra (J1Co)	Dolomites	0,5	Wackstones to grainstones.	Less than 10% grain > 2mm/More than 10% grain > 2mm; Presence of carbonated microcrystalline matrix; non-contact grains/absence of carbonated microcrystalline matrix and grains in contact. Poor to reasonable porosity (evaluation by type of lithology and texture).	0,7	Cross-bedding and dolomites with parallel or wavy lamination, interstratified with pellets.	0,4	Thickness about 60 m (beds with vertical or slightly inverted slopes).	0,1	1,7 (Moderate)
Platelet dolomites (J1pi)	Dolomitic limestones	0,7	Mudstone	Micritic dolomitic limestone. Poor porosity (evaluated based on the type of lithology).	0,2	Layers with centimeter to decimeter thickness.	0,1	Thickness about 30-40 m.	0,2	1,2 (Low)
Marls of Dagorda (J1Da)	Sandy loams, gypsum and saliferous clays, Intercalations of dolomitic limestones	0,1	Sandy, pelitic and Mudstones facies for the carbonate ones.	Poor porosity (evaluated based on the type of lithology).	0,1	250-320 m (dolomitic member, essentially dolomitic or margo-dolomitic with red and/or greyish pelites and evaporites); 60-850 m (saliferous/dolomitic member, predominantly dolomitic and/or limestone and marl rich in evaporites - anhydrite and halite); 290-800 m (saliferous member, characterised by an accentuated domain of halite, sometimes interstratified with dolomitic marls and/or marly pelites and anhydrite); From about 1000 m (occurrence of “Dolomites in platelets”); 1600-1200 m (evaporitic salts with intercalations of evaporitic syngenetic dolomite and gypsum); ~3000-1600 m (thick saliferous series with frequent clayey intercalations, also saliferous, containing anhydrite inclusions);	0,1	Formation subjacent Jurassic limestones with significant thickness, >3000 m, according to sounding “São Mamede 1”.	0,4	0,7 (Very low)
Eruptive rocks	Dolerite	0	-	Poor porosity (evaluated based on the type of lithology).	0	Associated with fractures or discontinuities.	0,8	-	0	0,8 (Very low)

^A Information collected from MANUPPELLA et al. (2000) and KULLEBERG et al. (2006, 2013). For Middle Jurassic units (Aalenian to Batonian), the facies described were proposed by AZERÊDO (1993, 1998). ^B Susceptibility to karsification, weighting 0 to 1. ^C Granulometry and texture according to the classification of DUNHAM (1962) and EMBRY & KLOVAN (1971) in SALOMON, 2000. Evaluation for primary porosity based on the classification of qualitative porosity adopted by INÊS (2010), from the work of AHR (2008): poor (0-5 %); reasonable (5-10%); good (10-15%); good/excellent (15-20%) and excellent (> 20%). For lithologies that were not evaluated by INES (2010), an evaluation was made based on the type of lithology and/or texture of the carbonate facies. ^D Information collected from MANUPPELLA et al. (2000). ^E Information collected in CRISPIM, 1995; Geological Map of Portugal 27-A Vila Nova de Ourém, 1998; MANUPPELLA et al., 2000. ^F Values from 0 to 4, with the following weighting: 0-1 (Very low); 1-1.5 (Low); 1.5-2.3 (Moderate); 2.3-3 (High); 3-4 (Very high).

Related to fractures, a greater spatial density, 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 infiltration of water in depth and the possibility to guide the development of caves (Cunha, 1988; Klimchouk et al., 2000; Ford & Williams, 2007). The fracture density calculation was performed using the ArcGIS Line Density tool (ESRI^TM), with the density quantified in km of fracture per km². We considered that the greater or lesser propensity for water infiltration is associated with the fracture type; we used slightly different criteria to calculate the density between the groups: “faults” and “joints + lineaments”. In the ArcGIS Line Density tool, we applied a search radius of 200 m to the “faults” group, assuming a greater concentration of water infiltration; in the group “joints + lineaments”, we applied a survey radius of 600 m, assuming a more diffuse concentration of water. The result of this exercise was two raster maps used for the definitive calculation of the fracture density through the formula - using the ArcGIS raster calculate tool.

$$Fracturedensity=\frac{1}{2}[\left(diaclasesandlineaments\right)+(faults\mathrm{x}3\left)\right]$$

(1)

For this factor, we classify the results into five qualitative classes (from Very low to Very High) that represent every class of fracture density (km/km²) classified in ArcGIS through the natural break method (Figure 5).

The energy available for the development of karst caves has a component that represents the difference between the altitude of the terrain and the altitude of the karst spring (or karst springs), associated with the local base level in a specific groundwater circulation sector (= drainage sector) of a carbonated massif (Salomon, 2000; Ford & Williams, 2007). This indicator is a relevant sign for assessing the importance of topography in underground karstification. It is mentioned in the present work as relief energy. Figure 5 represents the relief energy in five classes obtained through natural breaks, corresponding to five qualitative classes. We assume that the altimetric difference is a driving force (potential energy) proportional to a specific area’s altimetric difference. This force influences the hydrodynamic capacity of the waters that run through the interior of the rock mass, with clear influences on the development of karstification. The concept of relief energy is the role of topography as a driving force in karstification (Salomon, 2000; Ford & Williams, 2007; Moradi et al., 2016).

The type of land cover also influences karstification. In areas covered by soil and vegetation, meteoric waters are acidified, and their percolation is slower, which enhances the dissolution of carbonate rocks. We consider that the dissolution is enhanced in forest areas, in areas occupied by pastures and agricultural fields. In bush areas and uncovered areas, the karstification dynamic is closer to that of bare karst, as the influence of these soil occupations on chemical properties is reduced (although there is some influence in bush areas) (James & Choquette, 1984; Trudgill, 1985; Moradi et al., 2016). In artificialised territories, soil sealing does not enhance karstification; however, this last land cover can only influence the development of current karstification and not in the endokarst forms developed over geological time without any influence of anthropic action. We consider the land cover categories from the original National Land Cover Map (Figure 5).

The known entrances of caves will allow us to implement verification of the predictive capacity of the cartographic model built by correlating the results of the endokarst potential with the location of the referred entrances and by the analyse of ROC curves (Bi & Bennett, 2003; Ghung & Fabbri, 2003; Guzzetti et al., 2005; Oliveira, 2012).

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

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

In the first phase of the AHP implementation, we weighted the classes of the four karstification factors through the pairwise comparison; after that, we made a pairwise comparison between those factors. The quality of these comparisons is described through the consistency ratio (CR), which, in all cases of this work, has a value lower than 0,10 - indicative of satisfactory results. The relative weight of each factor (Wi) gives us the order of priority for each class associated with a factor and between factors (Table 3 to Table 7), there may be cases in which two classes have the same weight because they are considered to have the same importance, such as the relief energy classes, “Very high” and “High”. As expected in factors classes comparison, for lithostratigraphic units, fracture density and relief energy, the “Very high” and “High” classes have more than 60% importance to underground karstification; in land cover, the category “Forests” represents ~30% and the category “Pastures” 22% (Table 3 to Table 6). In comparing factors, the lithostratigraphic units have almost 50% importance, the fracture density 31%, relief energy 13% and land cover residual importance (8%) (Table 7).

Table 3. Pairwise comparison for susceptibility of karstification classes in lithostratigraphic units.

Lithostratigraphic units (SK*)	Very high	High	Moderate	Low	Very low	$Wi$
Very high	1	2	5	7	9	0,46
High	1/2	1	4	6	8	0,32
Moderate	1/5	1/4	1	3	6	0,13
Low	1/7	1/6	1/3	1	3	0,06
Very low	1/9	1/8	1/6	1/3	1	0,03
CR = 0,055 * susceptibility of karstification

Table 3. Pairwise comparison for susceptibility of karstification classes in lithostratigraphic units.

Lithostratigraphic units (SK*)	Very high	High	Moderate	Low	Very low	$Wi$
Very high	1	2	5	7	9	0,46
High	1/2	1	4	6	8	0,32
Moderate	1/5	1/4	1	3	6	0,13
Low	1/7	1/6	1/3	1	3	0,06
Very low	1/9	1/8	1/6	1/3	1	0,03
CR = 0,055 * susceptibility of karstification

Table 4. Pairwise comparison for fracture density classes.

Fracture density	Very high	High	Moderate	Low	Very low	$Wi$
Very high	1	2	4	6	9	0,47
High	1/2	1	3	4	7	0,30
Moderate	1/4	1/3	1	2	3	0,12
Low	1/6	1/4	1/2	1	1	0,07
Very low	1/9	1/7	1/3	1	1	0,05
CR = 0,013

Table 4. Pairwise comparison for fracture density classes.

Fracture density	Very high	High	Moderate	Low	Very low	$Wi$
Very high	1	2	4	6	9	0,47
High	1/2	1	3	4	7	0,30
Moderate	1/4	1/3	1	2	3	0,12
Low	1/6	1/4	1/2	1	1	0,07
Very low	1/9	1/7	1/3	1	1	0,05
CR = 0,013

Table 5. Pairwise comparison for relief energy classes.

Relief energy	Very high	High	Moderate	Low	Very low	$Wi$
Very high	1	1	2	3	4	0,32
High	1	1	2	3	4	0,32
Moderate	1/2	1/2	1	2	3	0,19
Low	1/3	1/3	1/2	1	1	0,09
Very low	1/4	1/4	1/3	1	1	0,08
CR = 0,007

Table 5. Pairwise comparison for relief energy classes.

Relief energy	Very high	High	Moderate	Low	Very low	$Wi$
Very high	1	1	2	3	4	0,32
High	1	1	2	3	4	0,32
Moderate	1/2	1/2	1	2	3	0,19
Low	1/3	1/3	1/2	1	1	0,09
Very low	1/4	1/4	1/3	1	1	0,08
CR = 0,007

Table 6. Pairwise comparison for land cover categories.

Land Cover	Forests	Pastures	Agriculture	Bushes	Open spaces or with little vegetation	Artificialised territories	$Wi$
Forests	1	2	3	3	3	8	0,35
Pastures	1/2	1	1	2	3	8	0,22
Agriculture	1/3	1	1	2	2	6	0,18
Bushes	1/3	1/2	1/2	1	1	7	0,11
Open spaces or with little vegetation	1/3	1/3	1/2	1	1	7	0,11
Artificialised territories	1/8	1/8	1/6	1/7	1/7	1	0,03
CR = 0,035

Table 6. Pairwise comparison for land cover categories.

Land Cover	Forests	Pastures	Agriculture	Bushes	Open spaces or with little vegetation	Artificialised territories	$Wi$
Forests	1	2	3	3	3	8	0,35
Pastures	1/2	1	1	2	3	8	0,22
Agriculture	1/3	1	1	2	2	6	0,18
Bushes	1/3	1/2	1/2	1	1	7	0,11
Open spaces or with little vegetation	1/3	1/3	1/2	1	1	7	0,11
Artificialised territories	1/8	1/8	1/6	1/7	1/7	1	0,03
CR = 0,035

Table 7. Pairwise comparison for karstification conditioning factors.

Karstification factors	Lithostratigraphic units	Fracture density	Relief energy	Land cover	$Wi$
Lithostratigraphic units	1	2	4	5	0,49
Fracture density	1/2	1	3	4	0,31
Relief energy	1/4	1/3	1	2	0,13
Land cover	1/5	1/4	1/2	1	0,08
CR = 0,018

Table 7. Pairwise comparison for karstification conditioning factors.

Karstification factors	Lithostratigraphic units	Fracture density	Relief energy	Land cover	$Wi$
Lithostratigraphic units	1	2	4	5	0,49
Fracture density	1/2	1	3	4	0,31
Relief energy	1/4	1/3	1	2	0,13
Land cover	1/5	1/4	1/2	1	0,08
CR = 0,018

The geographic information layers (in raster format), corresponding to each karstification factor (Figure 5), are considered to model the endokarst potential; after reclassification, their classes coincide with the normalised Wi. In a GIS, we calculated the endokarst potential by applying the following formula:

$$Endokarstpotential=\left({Wi}_{ClassesLitostratigraphicunits}\times {Wi}_{Litostratigraphicunits}\right)+\left({Wi}_{ClassesFracturedensity}\times {Wi}_{Fracturedensity}\right)+\left({Wi}_{ClassesReliefenergy}\times {Wi}_{Reliefenergy}\right)+\left({Wi}_{ClassesLandcover}\times {Wi}_{Landcover}\right)$$

(2)

The adequacy of the model was verified through the spatial confrontation of the results with the known natural caves. The verification of the predictive capacity of the model was based on the analysis of the ROC curves and the AUC (area under the curve), elaborated through the results of the binomial endokarst potential - location of the entrance of the caves for both the training and test subareas, as well as for the entire study area. In the ROC curve, the endokarst potential was represented in decreasing order on the abscissa axis and a cumulative distribution function of the inventoried cavities (with the data split into 20 classes) 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 has a better performance; if the curve inflects in the negative direction, the model has a poor performance, being also considered unacceptable (Bi & Bennett, 2003; Ghung & Fabbri, 2003; Guzzetti et al., 2005; Oliveira, 2012). The AUC value allows the quantitative assessment of the model’s predictive capacity. The AUC values are between 0% and 100% (or between 0 and 1), and the closer to 100%, the greater the predictive capacity of the model. The value of 50% graphically coincides with a diagonal line representing a casual predictive capacity, as mentioned above. Values below 50% represent a worse predictive capacity than random, and, in this case, the respective model cannot be considered acceptable (Bi & Bennett, 2003). For values that allow the assessment of the predictive capacity of the cartographic model, we opted for the values adopted in works related to natural hazards (e.g., susceptibility to landslides), namely by Guzzeti et al. (2005) and Oliveira (2012). Values between 75% and 80% correspond to an acceptable model, values between 80% and 90% indicate a very good model and values greater than 90% represent an excellent model.

Following the flowchart in Figure 4, we work iteratively on the predictive model until we obtain acceptable AUC values. Between iterations, we change some aspects related to the factors (e.g., calculation method, increment of quality of base data) and how we evaluate the weights in the AHP pairwise comparison. This type of approach has led us to obtain promising results.

4. 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). No cave entrances in areas classified as “Very Low” endokarstic potential exist. Although the promising result, we must consider that in a considerable area to the NW of São Bento and other sectors distributed throughout the study area, the endokarst potential is also classified as “Very High” or “High”, and here, there are no known cave entrances. Considering the present assessment as valid, we can say there may also be caves in those sectors. We also observe that areas with higher elevation have greater endokarst potential, existing in these areas some differences that may explain the distribution of known cave entrances. On the other hand, cave data only refer to the location of cave entrances that may belong to complex underground tunnel systems. With topographic data that allowed mapping the development of these systems and comparing them with the model produced, we could perform a more adequate analysis of the results achieved. Furthermore, the inventory of the location of the cave entrances concerns the current knowledge (at least until February 2021), and this distribution can be modified with the discovery of new caves and thus condition the validity of the model.

We noticed a correlation between the geologic units that are most susceptible to karstification with the distribution of caves, namely in the Micritic Limestones of Serra de Aire (J²_SA) and the Limestones of Chão das Pias (J²_CP), which represent 38.2% of the study area and contain about 63% of the caves. The Marly limestones and marls of Fórnea (J^1-2_Fo), Marls and marls of Zambujal (J²_ZA), Bioclastic limestones of Codaçal (J²_Co), which occupy 31.26% of the area, contain about 35% of the caves. Analysing these data, we conclude that the propensity for karstification of units belonging to the Middle Jurassic is confirmed, which are also classified as most susceptible to karstification (Table 2).

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

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

Regarding the remaining factors (Table 8), we present data with the number of caves for each class and their density per Km². Contrary to what would be expected, we observed that the “Moderate” and “High” fracturing density classes with a more significant number of caves instead of the “Very High” class. Even so, we verified a higher density of caves in the “High” and “Very High” classes. When we analysed the spatial relationship between the relief energy and the location of the known cave entrances, we verified a more significant number of these in areas with greater potential energy relief, such as in the “High” and “Very High” classes. It should also be noted that the highest density of caves is higher in the “Very high” class, which leads us to conclude that there is a greater concentration of caves in sectors with higher relief energy. Finally, in the land cover factor, we verified different values than expected; it is in the “Bushes” class where a more significant number and greater density of caves are observed, instead of “Forests”, “Pastures”, and “Agriculture”. These are apparent inconsistencies between the weighted classes for the factors now analysed and the (expectable) location of the caves; we must relativise these results because of the low importance attributed to land cover and the goods results obtained in the evaluation of the model.

Figure 7. Percentage of inventoried caves by karstification potential class in the study area.

Figure 7. Percentage of inventoried caves by karstification potential class in the study area.

Table 8. Natural inventoried caves by underground karstification factor classes.

Factors (criteria)	Classes	Area (km2)	N.º caves	N. º caves/Km2
Fracture density	Very low	21,99	5	< 1
	Low	14,67	6	< 1
	Moderate	9,55	11	1,15
	High	5,24	19	3,63
	Very high	2,55	7	2,75
Relief energy	Very low	11,66	0	0
	Low	9,89	2	< 1
	Moderate	8,15	1	< 1
	High	11,81	9	< 1
	Very high	12,49	36	2,88
Land cover	Forests	7,22	2	< 1
	Pastures	7,81	4	< 1
	Agriculture	14,71	1	< 1
	Bushes	20,54	40	1,95
	Open spaces or with little vegetation	0,08	0	0
	Artificialised territories	3,64	1	< 1

Table 8. Natural inventoried caves by underground karstification factor classes.

Factors (criteria)	Classes	Area (km2)	N.º caves	N. º caves/Km2
Fracture density	Very low	21,99	5	< 1
	Low	14,67	6	< 1
	Moderate	9,55	11	1,15
	High	5,24	19	3,63
	Very high	2,55	7	2,75
Relief energy	Very low	11,66	0	0
	Low	9,89	2	< 1
	Moderate	8,15	1	< 1
	High	11,81	9	< 1
	Very high	12,49	36	2,88
Land cover	Forests	7,22	2	< 1
	Pastures	7,81	4	< 1
	Agriculture	14,71	1	< 1
	Bushes	20,54	40	1,95
	Open spaces or with little vegetation	0,08	0	0
	Artificialised territories	3,64	1	< 1

In addition to the results obtained in evaluating the model, we compared the endokarst potential with data from another work with a similar theme in the same area - Crispim (1992, 1995). As a result of this comparative exercise, we generally verified that the units categorised in this work as “Karstified Formations” overlapping areas in the classes of higher endokarst potential and that the decrease in the endokarst potential corresponds to a decrease in those “Karstified Formations”. For “Very high” endokarst potential there is an overlapping with “Karstified Formations” at about 98%; for “High” endokarst potential the “Karstified Formations” overlaps in ~92%; the “Moderate” endokarst potential overlaps in ~93% the group of “Karstified Formations” and “Weakly karstified formations”; for “Low” endokarst potential the overlapping with the group of “Moderately karstified formations” and “Weakly karstified formations” it is ~70%; finally the “Very low” endokarst potential class overlaps ~52% the units of the group “Moderately karstified formations”, “Weakly karstified formations”, “Non-karstified formations”, “Non-karstificable formations”.

As we already said, the evaluation/verification of the cartographic model produced was based on the ROC curve analysis. Lastly, following the flowchart in Figure 4, we applied the cartographic model iteratively until we obtained acceptable AUC values indicative of the model’s predictive capacity. Observing Figure 8, with the ROC curve of the study area, we notice that about 13% of the area classified with “Very high” endokarst potential (87% to 100%) explains 70% of the areas with caves and approximately 32 % of the area with lower endokarst potential (“Low” and “Very Low”) explains only about 5% of the areas with caves. The AUC values are for the training and test areas 78,3% and 86%, respectively. For the entire study area, the value is 81,9%, typifying the model as good or, in other words, with a good predictive capacity.

Figure 8. ROC curve for the study area.

Figure 8. ROC curve for the study area.

5. Conclusions

The endokarst potential of a carbonate massif depends on several conditioning factors that contribute differently to the underground karstification process and that, as a rule, almost always work interdependently. Details associated with the faciological and stratonomic characteristics and the geometry of the carbonated layers can be decisive for the 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 waters into the carbonate massif is intrinsically linked to the existence of fissures (fractures and bedding planes). The hydraulic gradient enhances the velocity of these waters, that is, by the altimetric difference between the input (e.g., cave entrances) and output points (karst springs) -being that this difference can be characterised through the relief energy. How the meteoric waters enter the rocky massif and how they travel, have particular importance. It is essential to characterise any existing coverings (e.g., detrital cover, soil, vegetation). The endokarst potential is associated with the greater or lesser development of subterranean forms (caves) that sometimes has an expression on the surface (e.g., cave entrances, karst springs).

The geological and geomorphological characteristics of the northern sector of the Santo António plateau proved to be decisive for the evaluation of the endokarst potential, as well as the knowledge of the location of the cave entrances for the evaluation/verification of the proposed cartographic model. Using the model in larger areas of the Estremadura Limestone Massif or other karstified massifs is essential for its legitimacy in generalising its application. It was also clear that the model is highly dependent on the quality/quantity of the information on which it is based, namely, knowledge of the endokarst.

We also highlight the importance of our results for the knowledge of the endokarst and what it represents for human activities. Human occupation of karst territories requires intelligent Territory Planning and Management; we see that the weaknesses inherent to the karst environment are often not considered, especially in water resources, natural hazards, or the extraction of mineral resources. Many examples show how karst dynamics can put people and property at risk. Furthermore, the karst environment has a peculiar (sometimes spectacular) landscape (at the surface and in deep) to enjoy; this landscape can be used for geotourism, scientific, and educational activities. The cartographic model built in this work can also represent an excellent methodology to define sectors in a carbonated massif with environmental vulnerability (associated with karst hydrogeology), sectors with more potential for hazards, or sectors that can be used in speleotourism activities.