1. Introduction

2. Study area

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.

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.

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

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.

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.

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.

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.

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.

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.

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.

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

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.

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.

Figure 8. ROC curve for the study area.

Figure 8. ROC curve for the study area.

5. Conclusions

References

1. Abdellah, A.B., Jamila, R., Ahmed, B., Elhassan, L., Samir, N. (2022). The Karst Landscapes of Beni Mellal Atlas (Central Morocco): Identification for Promoting Geoconservation and Tourism. Quaestiones Geographicae, Vol.41(3), pp.87-109.
2. Aleksandar, A., Nemanja,T., Slobodan, M. (2019). Karst geoheritage and geotourism potential in the Pek River lower basin (eastern Serbia). Geographica Pannonica, Vol. 23 (1), pp. 32-46.
3. Azerêdo, A.C. (1993). Jurássico Médio do Maciço Calcário Estremenho (Bacia Lusitânica): análise de fácies, micropaleontologia, paleogeografia. Tese de Doutoramento, Departamento de Geologia da Faculdade de Ciências da Universidade de Lisboa, 366 p.
4. Azerêdo, A.C. (1998). Geometry and facies dynamics of Middle Jurassic carbonate ramp sandbodies. West-Central Portugal, University of Birmingham, 34 p.
5. Azerêdo, A.C. (2007). Formalização da litostratigrafia do Jurássico Inferior e Médio do Maciço Calcário Estremenho (Bacia Lusitânica). Comunicações Geológicas, t. 94, pp. 29-51.
6. Ballesteros, D., Fernández-Martínez, E., Carcavilla, L., Jiménez-Sánchez, M. (2019). Karst Cave Geoheritage in Protected Areas: Characterisation and Proposals of Management of Deep Caves in the Picos de Europa National Park (Spain). Geoheritage 11, pp. 1919–1939.
7. Barton, H.A. (2013). Biospeleogenesis. Treatise on Geomorphology, J. Shroder (Editor in Chief), A. Frumkin (Editor), Academic Press, San Diego, CA Vol. 6, Karst Geomorphology, pp. 38-56.
8. Bennett, B., Ibrahim, A.M. (2014). The Assessment of Machine Learning Model Performance for Predicting Alluvial Deposits Distribution. Procedia Computer Science 36, pp. 637-642.
9. Bi, J., Bennett, K.P. (2003). Regression Error Characteristic Curves. Proceedings of the Twentieth International Conference on Machine Learning (ICML-2003), Washington DC, pp. 43-50.
10. Bolstad, P. & Manson, S. (2022). Chapter 13: Spatial Models and Modeling. In GIS fundamentals: A first text on geographic information systems, 7th ed., Eider Press, pp. 561-582.
11. Brilha, J. (2005). Património geológico e geoconservação a conservação da natureza na sua componente geológica. Palimage, Viseu, 190 p.
12. Canais, F., Fernandes, J. (1999). Lapas e Algares da Serra de Santo António. SubTerra - Grupo de Espeleologia 187 p.
13. Carvalho, J.M. (2013). Tectónica e caraterização da fraturação do Maciço Calcário Estremenho, Bacia Lusitânica. Contributo para a prospeção de rochas ornamentais e ordenamento da atividade extrativa. Tese de Doutoramento, Departamento de Geologia da Faculdade de Ciências da Universidade de Lisboa, 443 p.
14. Chiarini, V., Duckeck, J., De Waele, J. A. (2022). Global Perspective on Sustainable Show Cave Tourism. Geoheritage 14 (82), 27 p.
15. Church, R. L. (2002). Geographical information systems and location science. Computers and Operations Research 29 (6), pp. 541–562.
16. Crispim, J.A. (1992). Características cársicas das rochas carbonatadas da região central do Maciço Calcário Estremenho. Algar, Boletim da Sociedade Portuguesa de Espeleologia 3, Lisboa, pp. 9-18.
17. Crispim, J.A. (1995). Dinâmica cársica e implicações ambientais nas depressões de Alvados e Minde. Tese de Doutoramento em Geologia, Faculdade de Ciências da Universidade de Lisboa, 388 p.
18. Cunha, L. (1988). As Serras Calcárias de Condeixa-Sicó-Alvaiázere. Estudo de Geomorfologia. Tese de Doutoramento, Faculdade de Letras da Universidade de Coimbra, 325 p.
19. Cunha, L. (1993). Geomorfologia estrutural e cársica. O Quaternário em Portugal, Balanço e Perspectivas, Colibrí, Lisboa, pp. 63-74.
20. Cunha, L., Dimuccio, A., L., Aubry, T. (2014). Património cársico no Centro de Portugal. O caso do Maciço de Sicó. Património Geomorfológico e Geoconservação. Experiências em Portugal e no Brasil, Atas/proceedings do I Encontro Luso-Brasileiro de Património Geomorfológico e Geoconservação, pp. 27-34.
21. Cunha, L., Dimuccio, A., L., Paiva, I. (2018). Geopatrimónio cársico e desenvolvimento local no Maciço de Sicó. Desafios para afirmar a Lusofonia na Geografia Física e Ambiente, II Encontro Luso-Afro-Americano de Geografia Física e Ambiente, pp.1227-1233.
22. Daveau, S. (2004). O Relevo de Portugal. Grandes Unidades Regionais. In: Feio, M. & Daveau, S. (Organizadores), Coimbra, Associação Portuguesa de Geomorfólogos, pp. 61-73.
23. De Smith, M.J., Goodchild, M.F., & Longley, P.A. (2018). Geospatial analysis and model building. In Geospatial Analysis: A Comprehensive Guide to Principles, Techniques and Software Tools, 6th ed., 618 p.
24. Ferreira, A.B., Rodrigues, M.L., Zêzere, J.L. (1988). Problemas da Evolução Geomorfológica do Maciço Calcário Estremenho. Finisterra, XXIII, Vol. 23 (45), Lisboa, pp. 5-28.
25. Fielding, A.H., Bell, J.F. (1997). A review of methods for the assessment of prediction errors in conservation presence/absence models. Environmental Conservation, Vol. 24 (1), pp. 38-49.
26. Figueiredo, R.F. (2001). Modelação cartográfica em ambiente SIG para apoio à decisão: aplicação ao estudo da afectação potencial de usos do solo no sector Norte do Maciço Marginal de Coimbra. Dissertação de mestrado em Geografia Física e Estudos Ambientais, Faculdade de Letras da Universidade de Coimbra, 204 p.
27. Fleury, E. (1925). Portugal Subterrâneo, Ensaio de espeleologia portuguesa. Coleção Natura, Lisboa, 45 p.
28. Ford, D., Williams, P. (2007). Karst Geomorphology and Hydrology. John Wiley & Sons, Ltd, 562 p.
29. Gao, J. (2021). Fundamentals of Spatial Analysis and Modelling, 1st ed., CRC Press, 368 p.
30. Goepel, K. (2018). AHP Excel Template with multiple Inputs. Version AHPcalc-v2018-09-15.
31. Ghung, C.F. Fabbri, A.G. (2003). Validation of Spatial Prediction Models for Landslide Hazard Mapping. Natural Hazards 30, pp. 451-472.
32. Guzzetti, F., Reichenbach, P., Cardinali, M., Galli, M., Ardizzone, F. (2005). Probabilistic landslide hazard assessment at the basin scale. Geomorphology 72, pp. 272–299.
33. Hämäläinen, R.P., Salo, A.A. (1997). The Issue is Understanding the Weights. Journal of Multi-Criteria Decision Analysis, Vol. 6, pp. 340-343.
34. Hobléa, F. (2004). Expertise et médiation scientifique pour la gestion durable de l’environnement karstique. Actes de la réunion annuelle de la Société Suisse de Géomorphologie (SSGm), Gestion durable de l’environnement karstique, pp. 25–42.
35. Huang, B. (2018). Comprehensive geographic information systems (Huang, B., Editor in Chief). Vol. 1, GIS methods and techniques. Vol. Editors: Thomas J.C., Ming-Hsiang T., Elsevier, 439 p.
36. Inês, N. (2010). Paleotipologias diagenéticas em formações carbonatadas do Jurássico Médio da Bacia Lusitânica: análise multiscala e aplicações à caraterização de reservatórios. Dissertação de Mestrado na Área de especialização: Estratigrafia, Sedimentalogia e Palentologia, Departamento de Geologia da Faculdade de Ciências da Universidade de Lisboa, 230 p.
37. Ishizaka, A., Labib, A. (2009). Analytic Hierarchy Process and Expert Choice: Benefits and limitations. Operational Reasearch Society, Vol. 22, 4, pp. 201-220.
38. James, N.P., Choquette, P.W. (1984). Diagenesis 9, Limestones - The Meteoric Diagenetic Environment. Carbonate Sedimentology and Petrology, Vol. 4, pp. 45-78.
39. Khalaf, E.E.D.A.H. (2022). Karst Heritage as a Tourist Attraction: a Case Study in the White Desert National Park, Western Desert, Egypt. Geoheritage 14 (94) (2022), 30 p.
40. Klimchouk, A., Ford, D., Palmer, A., Dreybrodt, W. (2000). Spleogenesis, Evolution of Karst Aquifers. National Speleological Society January 2000 Edition, 527 p.
41. Kullberg, J.C. Rocha, R.B., Soares, A.F., Rey, J., Terrinha, P., Axerêdo, A.C., Callapez, P., Duarte, L.V., Kullberg, M.C., Martins, L., Miranda, R., Alves, C., Mata, J., Madeira, J., Mateus, O., Moreira, M., Nogueira, C.R. (2013). A Bacia Lusitaniana: Estratigrafia, Paleogeografia e Tectónica. In: Dias, R., Araújo, A., Terrinha, P. & Kullberg, C. (Editors), Portugal, Vol. II, Geologia Meso-cenozóica de Portugal. Escolar Editora, pp. 195-347.
42. Leslie, A.N. (2011). Informal Karst Education in the United States and Internationally. Ph.D. Thesis, Department of Geography, University of South Florida, 491 p.
43. Longley, P.A., Goodchild, M. F., Maguire, D.J., Rhind, D. W. (2011). Geographic Information Systems & Science. Third Edition, John Wiley & Sons, 539 p.
44. Martins, A.F. (1949). Maciço Calcário Estremenho. Contribuição para um Estudo de Geografia Física. Tese de Doutoramento. Faculdade de Letras da Universidade de Coimbra, 248 p.
45. Moradi, S., Kalantari, N., Charchi, A. (2016). Karstification Potencial Mapping in Northeast of Khuzestan Province, Iran, using Logic and Analytical Hierarchy Process (AHP) techniques. Geopersia 6, pp. 265-282.
46. Nola, I.T.S., Bacellar, L.A.P. (2021). Multi-criteria analysis for mapping susceptibility to iron formation caves development in the Gandarela mountain range (MG), southeast Brazil. International Journal of Speleology 50 (2), Tampa, FL (USA), pp. 137-187.
47. Oliveira, S.M.C. (2012). Incidência Espacial e Temporal da Instabilidade Geomorfológica na Bacia do Rio Grande da Pipa (Arruda dos Vinhos). Tese de doutoramento em Geografia Física, Universidade de Lisboa, 496 p.
48. Paiva, I.M.R. (2014). Hidrossistema Cársico de Dregracias-Sicó. Estudo do funcionamento hidrodinâmico a partir das suas respostas naturais. Tese de doutoramento em Geografia, Faculdade de Letras da Universidade de Coimbra, 445 p.
49. Ramos, A., Cunha, L., Cunha, P. (2014). Application de la Méthode de l’Analyse Multicritère Hiérarchique à l’étude des glissements de terrain dans la région littorale du centre du Portugal: Figueira da Foz – Nazaré. Geo-Eco-Trop, Liège 38 (1), pp. 33-44.
50. Rauch, H.W., White, W.B. (1970). Lithologic controls on the development of solution porosity in carbonate aquifers. Water Resources Research 6 (4), pp. 1175-92.
51. Rocha, R.B., Soares, A.F. (1984). Algumas reflexões sobre a sedimentação jurássica na orla meso-cenozóica ocidental de Portugal. Mem. Notícias Universidade Coimbra 97, pp. 133-142.
52. Rodrigues, L.M.E. (1998). Evolução geomorfológica quaternária e dinâmica actual. Aplicações ao ordenamento do território. Exemplos no Maciço Calcário Estremenho. Tese de Doutoramento em Geografia Física, Universidade de Lisboa, 868 p.
53. Rodrigues, M.L. (2020). The Limestone Massif of Estremadura. In: Vieira, G., Zêzere, J., Mora, C. (Editors) Landscapes and Landforms of Portugal. World Geomorphological Landscapes. Springer, pp. 229-250.
54. Rodrigues, H.A., Rodrigues, M.L., Lobo, H. S. (2017). Espeleoturismo em Portugal: Panorama Geral do Uso Turístico das Cavidades Naturais. Revista Rosa dos Ventos - Turismo e Hospitalidade 9 (I), pp. 92-106.
55. Saaty, T.L. (1977). A scaling method for priorities in hierarchical structures. Journal of Mathematical Psychology 15 (3), pp. 234-281.
56. Saaty, T.L. (1986). Axiomatic foundation of the analytic hierarchy process. Management Science 32 (7), pp. 841-855.
57. Saaty, T.L. (2009). Theory and Applications of the Analytic Network Process: Decision Making with Benefits, Opportunities, Costs, and Risks. RWS Publications, Pittsburg, PA USA, 352 p.
58. Saaty, T.L., Ergu, D. (2015). When is a Decision-Making Method Trustworthy ? Criteira for Evaluating Multi-Criteria Decision-Making Methods. International Journal of Information Technology & Decision Making, Vol. 14, pp. 1-17.
59. Salomon, J. N. (2000). Précis de Krastologie. Presses Universitaries de Bourdeaux, Pessac, Bourdeaux, 249 p.
60. Seif, A., Ebrahimi, B. (2014). Using GIS to evaluate degree of karstification according to some important factors in carbonate rocks in Iran. Carbonate Evaporites 29, pp. 107-126.
61. Taheri, K., Gutiérrez, F., Mohseni, H., Raeisi, E., Taheri, M. (2015). Sinkhole susceptibility mapping using the analytical hierarchy process (AHP) and magnitude-frequency relationships: A case study in Hamadan province, Iran. Geomorphology vol. 234, pp. 64-79.
62. Thomas, C. (1985). Cova da Velha. Espeleo Divulgação4, Núcleo de Espeleologia da Universidade de Aveiro, pp. 30-35.
63. Tomlin, C.D. (1990). Geographic Information Systems and Cartographic Modeling. Englewood Cliffs N.J., Prentice-Hall (now Pearson Education), p. 572.
64. Tomlin, C.D. (2017). Cartographic Modeling in International Encyclopedia of Geography. American Cancer Society, pp. 1-6.
65. Torres, E., Spoladore, A. and Cunha, L. (2018). Espeleoturismo: nota preliminar sobre as realidades paranaense e portuguesa. Paisagem e Território, Londrina, UEL, pp. 83-88.
66. Trudgill, S.T. (1985). Field observations of limestone weathering and erosion in the Malham District, North Yorkshire. Filed Studies 6, pp. 201-236.
67. Vakilipour, S., Sadeghi-Niaraki, A., Ghodousi, M., Choi, S. (2021). Comparison between Multi-Criteria Decision-Making Methods and Evaluating the Quality of Life at Different Spatial Levels. Sustainability 13 (7), 4067, 36 p.
68. Valjavec, M.B., Pejnović, N.D., Draženović, M., Čonč, Š., Horvat, K.P. (2022). The Transboundary Approach to Landscape Geointerpretation: Challenges in Interpretive Planning and Geoconservation. Geoheritage 14 (116), 13 p.
69. Vargas, H.S., Almeida, J.A., Kullberg, J.C., Brito, M.G. (2003). Prospecção interactiva de endocarso através de uma aplicação SIG. Finisterra XXXVIII (76), pp. 51-64.Velasquez, M., Hester, P.T. (2013). Na analysis of multi-criteria decision making methods. International Journal of Operations Research 10 (2), pp. 56-66.
70. Waltham, A.C. (1981). Origin and development of limestone caves. Progress in Physical Geography: Earth and Environment Vol. 5, pp. 242-256.
71. White, W.B. (1988). Geomorphology and hidrology of karst terrains. O.U.P. New York, 464 p.
72. Zaree, M., Javadi, S., Neshat, A. (2019). Potential detection of water resources in karst formations using APLIS model and modification with AHP and TOPSIS. Journal of Earth System Science 128 (76), 12 p.
73. Zlaugotne, B., Zihare, L., Balode, L., Kalnbalkite, A., Khabdullin, A., Blumberga, D. (2020). Multi-Criteria Decision Analysis Methods Comparison. Environmental and Climate Technologies, Vol. 24 (1), pp. 454-471.