UAV-based doline mapping in Brazilian karst: A cave heritage protection reconnaissance

1.1 Study area

The study site description is a compilation of previous works [2,18,19,20,21] focusing on the fluvio-karst structures in the APANRV. The studied site covers an area of about 49.47 km² within the APANRV, situated in the states of Goiás and Bahia (Figure 1a). Sedimentary rocks, especially siliciclastic from the Urucuia Group, along with alluvial residues from the Serra Geral de Goiás dismantling were found abundantly in the region. Older carbonate rocks of the Bambuí Group (Neoproterozoic) occur in a subjacent way, such as the Lagoa do Jacaré formation, that is characterized by the intercalation of thin pelitic layers and the presence of cherts. Carbonates outcrop only in the middle and lower parts of the relief, where the unconsolidated siliciclastic materials have been removed. The study area can be classified into covered karst [22] and is in transition to being exposed to the surface by the actions of fluvial system. This transition is triggered by the high hydraulic gradients created by the lowering in regional base level. This is well illustrated by the presence of a deep trough of the Vermelho River, which resulted in the formation of a canyon of 80 m high walls [21]. This new condition seems to have dynamized the local karst system through caves interconnection, emptying previously filled galleries (paragenesis), and formation of deep veined canyons. There are also captures of surface drainage into the underground environment and the intensification of erosive processes near the gullies, especially the dolines (Figure 2). One of the largest caves (Tarimba cave > 10 km in horizontal projection and 69 m deep) in Brazil falls within this context, with its headwaters inducing erosion processes and the formation of dolines [23].

Figure 1

Location of the study area. (a) Location of the overflight area (UAV) within the APANRV. (b) Subareas of UAV aerial survey, control, and verification points, trails traveled and dolines identified in the field.

Figure 2

The major doline typologies found in the study area. (a) suffosion-type doline, (b) suffosion-type doline under forest, (c) collapse doline, (d) bottom of a cockpit-type depression with access to a cave, (e) buried doline, and (f) cover collapse doline.

The entire watershed of Serragem stream, which divides the study area approximately in half, is considered in the analysis. The Extrema stream delimits the area to the southeast, together with the river Vermelho stretching to the south; therefore, only a part of its watershed falls in the study area. Around 8 km² area is under the unconsolidated materials domain associated with the Urucuia Group, where dolines are not expected to occur significantly. This sector corresponds to the highest altimetric elevations in the region, in the north-northeast area portions. It also corresponds to the most preserved vegetation areas of the Cerrado type (savanna). In the lower portions, where there is greater availability of water, gallery forests and a few areas of Cerrado occur, which are being replaced by pastures and crops associated with rudimentary practices. The climate of the region is tropical, with dry (April to September) and rainy (October to March) seasons presenting a precipitation index of ∼1,260 mm/year [20].

2.1 UAV surveying and image processing

The Phantom 4 Pro DJI UAV vehicle was used for image acquisition (Table 1 – characteristics) over the two subareas (Area-1 and Area-2) for better execution of flight plans (Figure 1b). In parallel, a strategy was adopted for collecting data with a Hiper X dual-frequency GNSS receiver (RTK) at the points to be used for geometric correction. The UAV-based aerial images and control points were processed in Agisoft Photoscan software. The control points optimized the image alignment, resulting in a dense point cloud used in the digital surface model (DSM), digital terrain model (DTM), and orthomosaic photo generations (Table 1).

2.2 Field survey for validation

The field validation of the identified dolines was done by field visits to various prominent features to ensure better representation and spatial distribution of different features. As recommended in the literature, the coordinates were taken at the base of the dolines, in the deepest point [6,24]. After the UAV generated objects (dolines and possible dolines) for analysis, a GNSS navigation receiver was used to position true dolines (validated in the field) and characterize possible commission errors (Figure 1b). One hundred dolines were located and characterized (typology) in the field at a minimum density value of 2.02 dolines/km².

2.3 Automatic doline detection on DTM

The detection starts with the DTM filling using the “Fill” tool (ArcMap, Hydrology). No “Z” limit is specified, which may carry a total filling of the depressions or model artifacts (here called Z-0). The recommendation for applying Z-limits is discussed elsewhere [2]. Next the subtraction of the original DTM from the filled DTM was made, resulting in a new DTM of depressions. Similar fill-difference routines have been documented in the literature [7,25]. As the number of objects that resulted from this operation was very high, the depressions deeper than 1 m were eliminated. The noise due to isolated pixels in the raster/polygon fragmentation process (single-cell depressions) and the overlaps were also eliminated (base Z-1) (Table 2).

The result was confronted with the orthomosaic photo and shaded relief (hillshade, DTM) for visual evaluation of the method’s effectiveness in detecting authentic dolines. A large number of pseudo-objects (dolines or no-dolines) that appeared as typical dolines on the shaded map, were taken as tree canopies using photos (Figure 3).

Figure 3

FPs associated with vegetation. At the top, in the Shaded DTM, it is possible to observe what appear to be small collapse dolines, outlined in red in the automated process. In the center is the topographic profile of the blue dashed line in the image above (A and B). Finally, below is the same scene in the orthophoto-mosaic, with the objects in red over the treetops. UTM 23 S, SIRGAS 2000.

Filtering the minimum size of dolines, non-propitious lithologies, drainage, and roads [9] is ineffective in reducing false positives (FPs). The morphometric analysis of these FP polygons can provide information about area and perimeter. In addition, with the QGIS zonal statistics tool, data were collected on the variation in the pixels internal to each object created from the DTM (e.g., minimum elevation, mean, range, and standard deviation). That was done using the DSM base without the interference caused by the automatic classification of the dense point cloud.

The statistical pixel-by-object analysis of DSM data allowed the identification of some peculiar behaviors of the depressions over treetops as they tended to be deeper and larger in perimeter compared to the area. Testing with arbitrary values of amplitude-to-area (ATA) ratio and perimeter-to-area ratio PAR can prove effective in the elimination of these spurious objects. In general, dolines follow a proportional relationship between depth and size, more like a circular geometry than other irregularly shaped objects [26]. Such relationship (depth-to-diameter ratios) has also been utilized in many karst morphometry studies aiming to understand the doline origin and evolution through their characterization [6].

For the morphological analysis, the cut-off limits are chosen based on field data and depth area relationship, as recommended by Rahimi and Alexander [27]. This way, the objects of true dolines (Z-1 DTM) were identified, and the highest value of the depth and perimeter per area ratios were found. Pixel statistics internal to the objects were performed using DSM and DTM data. These operations allowed the elimination of considerable amounts of abruptly shaped objects compared to the previous filtering procedures which resulted in the delineation of about 11,592 objects (Table 2). It is noteworthy that the identification of fewer objects in the DTM (44 polygons) as being relative to true dolines in the field (100 features) is mainly due to the coalescence tendency in compound depressions [28]. Automatic detection by fill occurs up to the level of overflow, which may imply the joining of adjacent smooth depressions.

For the redundancy checking, a cluster analysis including four tests, K4x4-DSM, K4x4-DTM, K3x5-DSM, and K3x5-DTM, in QGIS attribute-based clustering tool using non-altimetric parameters from the pixel statistics were performed. The abbreviation K4x4 refers to the adoption of the four parameters, including area, perimeter, range, and standard deviation and the distribution into four mandatory classes. Similarly, K3x5 refers to three parameters (e.g., area, perimeter, and range) and their distribution into five mandatory classes. As a result, the objects identified with Z1-DTM processing are found to represent the format better, but the pixel database varies for DSM and DTM.

The K-means clustering used in the attribute-based clustering tool (QGIS) is a non-hierarchical clustering method that requires a previously defined number of groups (clusters). It minimizes the sum of the squares of the distances of each object from the centroid of the group it belongs to [29]. For the optimal number of clusters determination, the Factoextra package of the R statistical program was used (R Core Team 2019). The method chosen was the intra-cluster sum of squares (“within-cluster sums of squares”), with several clusters ranging from 1 to 10. Given the result, it was decided to continue the analysis with four and five clusters. The tests performed in the cluster analysis have been summarized in Table 3, where it is possible to observe the variation in the parameter/group arrangements chosen for comparison based on the validated process with 44 dolines observed in the field.

A confusion matrix was used to evaluate the results of the tests (morphological and clustering), as documented in the literature [7,10,27,30,31]. This matrix type compares the data collected in the field with the data predicted by the processing performed to separate false objects from dolines. The true positives are the points where there is compatibility between real dolines and dolines predicted with the model. Similarly, the true negative (TN) points present high compatibility with no doline points. FPs also occur, where the model points to erroneous objects inferred as no dolines from the field surveying results. Finally, true dolines can be classified as no dolines by the model, called false negatives (FNs).

The best cluster analysis results served as the data source for the confusion matrix (Table 3). The highest occurrence of true dolines in one of the classes was adopted as a true positive, considering that the cluster analysis was able to discern that set as belonging to the same type of objects. The occurrences of no dolines verified in the field that fell into the class with the highest result of true dolines were considered a FP. A high value of true positive (TP) with a low FP tends to be considered a better model for the overall filtering of the objects. However, the confusion matrix provides for considering the negative parameters as well. TNs considered those no doline objects (verified in the field) that fell into all classes other than true dolines (which had more correct objects). Finally, FNs are those true dolines that were classified in classes other than those with the highest recurrence of dolines, i.e., they were controversially classified as no dolines.

To better balance the analyses in different samples (44 dolines vs 403 non-dolines), the overall accuracy, the TP and TN rates, and the average accuracy were calculated. In the cluster analyses, more significant results were observed for the K3x5-DSM database. It presented rates consistently above 80%, an average accuracy of almost 89%, and relatively few objects remaining for analysis. Despite the fewer objects of the K4x4-DTM model, it presented much lower indices.

For the morphological analyses, the best results were with the perimeter-to-area ratio (PAR), achieving an average accuracy of over 96%. This shows that at least for characterizing spurious objects, the product of this ratio can be quite effective (93%), since the data from true dolines were considered 100% correct. When comparing the performance of this ratio with the K3x5-DSM model, the latter is more effective in TN rate (96% vs 93%) and overall accuracy (94.6% vs 93.7%). This means that the better result for average accuracy of the PAR only occurred due to the unbalance in adopting all objects considered to be dolines (TPs) representative of this phenomenon. Nevertheless, considering only the tests of morphological analysis, it also resulted in a much lower number of objects for visual analysis. In addition to the PAR, having obtained better rates than the others (amplitude/area-DTM and amplitude/area-DSM).

The objects generated from the best clusters and morphological analyses were gathered in the same base for future visual inspection (Table 4). Of the 933 remaining objects, 41% were obtained from the K3x5-DSM base, 25% from the PAR, and 34% convergence between the two analyses performed. Additionally, the confirmed FPs were eliminated.

2.4 Automatic detection of dolines (DSM)

Considering that dolines were identified in the field, but did not generate polygons in the DTM procedures, the difference by filling in the DSM was also performed. This measure aimed to detect milder depressions in the open field or under the forest, whose vegetation eventually followed the relief. The automatic classification from point cloud to DTM could misrepresent important information.

Similarly, vast amounts of fill difference objects were created in the DSM. The adoption of a depth filter (e.g., Z-1) did not prove to be meaningful. Considering that the shape of the spurious objects was more irregular compared to real dolines (Figure 4), and based on field data, a maximum PAR (0.4) was chosen, which resulted in a significant cleanup (Table 5). This new base, considered ancillary, was kept independent of the objects generated from the DTM. Despite the locational correspondence of many objects between these bases, due to the significant difference in format, it was considered that they should not be brought together for visual analysis.

Figure 4

Shape difference between doline (blue) and spurious objects (red) with DSM data. The spurious objects have a much higher product of PAR than the dolines observed in the field. UTM 23S, SIRGAS 2000.

2.5 Manual/visual procedures

Objects were classified in the visual analysis between no dolines and possible dolines, with support from the orthomosaic photo, topographic profiles, and shaded relief (DTM and DSM). The first group consists of objects associated with remnant tree canopies, depressions linked to drainages, edges of sloping surfaces between pastures and forests (especially in the DSM base), and anthropogenic depressions associated with dams and roads. Similar object-cleaning procedures have been extensively adopted in literature [7,8,9,10,32].

In addition to the objects already confirmed as dolines (44), possible karst depressions were classified. Objects in which the visual elements indicate a chance of approximately 90% or more of being true was called “probable,” while objects with probabilities above 75% were classified as “suspected.” The presence of centripetal drainage, vegetation change, morphological aspects, sinkholes, and carbonate rocks, among others, have contributed to the greater conviction. There are, however, inconclusive objects (50% chance), called shallows, amid dense vegetation or woods that may be non-apparent drainages, deep clearings, or even sinkholes, as many have seen in the field (Figure 5). Figure 6 summarizes the steps of this research.

Figure 5

Visual classification of objects. (a) Objects associated with drainages, generally elongated in which it is possible to observe the channel; (b) objects associated with vegetation edges at different altimetry, common at the pasture/woodland junction; and (c) objects called shallows, in which it is not possible to visualize the deepest point and unequivocally decide for a doline or not. In this case, the shallow is 10 m deep. UTM 23S, SIRGAS 2000.

Figure 6

Schematic presentation of methods.

3.1 Field data

One hundred dolines identified in the field were duly characterized based on their typologies (simple or compound type), land use (e.g., such as forest, Cerrado, pasture, and exposed soil classes), the existence of water (stream, lake), and cave-ground connections. The genetic typologies most commonly used to classify dolines are dissolution, collapse, suffosion, cover collapse, and burial [2,33,34]. Additionally, the cockpit type can also be included here, such as those remodeled by surface runoff channels. Considering these processes can occur in nature in a combined and simultaneous way, this can be classified into predominant and secondary processes (Table 6).

The suffosion and cockpit are the predominant typologies, responsible on average for 37.8 and 31.7% of the total identified features in the region, respectively. Both are predominantly conditioned by the presence of cover materials (abundant in the region) and surface shaping processes (slope or fluvial systems). Rock collapse is another significant typology that accounts for approximately 20% of the features. The burial and cover collapse processes are more localized without the presence of any typical dissolution dolines.

At several places, these typologies can also occur in conjugation, such as cockpit type, followed by collapses and several suffosion dolines (Figure 7). An abrupt drop in elevation (15–20 m) together with lowering of the base level causes capturing of the drainage, collapses, and subsidence at the downstream side. This is frequent in gallery forests and paleochannels, showing that small streams can be captured and a surface aspect of active drainage under the forest can be maintained.

Figure 7

Typical sequence of dolines found in the area, often associated with caves. Dolines occur where siliciclastic deposits or soils over carbonates are thin.

The analysis of land use in the sample of dolines also indicated some main behaviors, among the land use types such as forest, Cerrado, pasture, and exposed soil classes (Table 7). Generally, forest predominates 62% of the time in the dolines, influencing 85% of features. This leads to a major challenge for the UAV and visible spectrum camera (RGB) application in doline detection in the study area. The pasture reached about 33%, presenting a predominant land use type and 59% as secondary. There were few dolines in Cerrado or exposed soil, which indicates the transition of the native Cerrado to the pasture land, especially, over the open areas, without the occurrence of abrupt gradients or exposure of the epikarst (outcropping rocks). In the forest areas, the deeper dolines occur (collapses, for example) and have been preserved because of unfavorable conditions for farming.

About 32% of the dolines had water inside (perennial streams or lakes) and resulted from the strong underground capture, as most of these are collapses (karst windows) or perennial drainage cockpits. On the other hand, 54% of the dolines in the area serve as cave-ground connection points, demonstrating their significance as an indirect indicator of the presence of the speleological feature. Slightly more than half is embedded in larger, compound depressions (51%), which is consistent with the recurring typology of cockpits, generally polygonal depressions containing lower-order dolines.

3.2 Remote sensing

After the analysis of the two databases, first the one from the DTM processing (K3x5-DSM and PAR) with 933 objects and then the one obtained from the DSM (perimeter-to-area), with 730 objects, the following results were obtained (Table 8). In general, most of the analyzed objects were categorized as negative, in an approximate and coincident relation of 70% for both bases. In the DTM and DSM-driven objects, most of these negatives are attributed to depressions over treetops or next to drainages and related to vegetation edges, respectively. This difference between the two bases may be related to the ground-level correction effects during the creation of the DTM. In the DSM base, no canopy-related depressions were observed. On the other hand, there were consequences of the presence of forest masses on slopes acting as virtual barriers, with the use of hydrological restitution tools in the fill-difference routine. Apparently, this could lead to the impression that the use of one base or the other would only imply the substitution of problems: either false depressions relative to tree crowns would occur, as in the DTM base, or spurious depressions close to the vegetation edge would occur, as in the DSM base.

The analysis of the positive results, which matters when considering the goal of detecting dolines, revealed that both analyzed bases resulted in similar relative amounts of suspected or probable objects, ranging from 7 to 8% in the first case and around 3% in the second. Polygons associated with dolines identified in the field also orbited close to 5% of the totals in both bases. However, contrary to expectations, the number of shallows was greater for the DTM base than for the DSM base. Since such depressions generally occur in dense forests, it was believed that the interpolation leveling processing in creating the DTM would eliminate such depressed areas more generally. This may be due to the existence of pronounced gradients within these forested areas, deciding to keep the shallows as positive for future field validation.

The joint analysis of the two bases allowed the identification of the object detection overlaps and, above all, the additional objects brought by the inspection of the accessory DSM base (Table 9). As can be seen, the quantity of positive objects added by the latter is significant, including the creation of three objects related to field dolines that had not been detected in the DTM base. The addition of more than 43% (122) in total positives by the DSM base reveals that the joint analysis of the DSM proved to be very productive.

The 47 polygons connected to dolines are associated with 78 of the 100 points verified in the field. This is due to the fill-difference technique sometimes encompassing lower-order dolines in composite depressions in their surroundings. Nevertheless, considering the limitations of the sensor used (RGB), we can observe a 78% success rate. On the other hand, dolines that were not detected by any database (22) were identified in the field. To ensure a better representation of the distribution of true dolines in the area along with the possible dolines, it was decided to transform the polygons into points and replace those related to dolines (47) with those collected in the area (100). Thus, the original 403 objects (Table 9) generated 456 points among shallows, suspected, probable, and dolines, as seen in the map shown in Figure 8a. This measure prevented a compound depression, with many internal dolines, from being represented by just one point in the density analysis (Figure 8b).

Figure 8

Depression distribution in the study area. (a) Distribution of the objects and their respective classification, in which it is possible to notice a greater occurrence near the drainage incisions. (b) Density of all possible objects (Kernel density, ArcGIS) and their relationship to large cave systems. UTM 23S, SIRGAS 2000.

The analysis of the presented maps (Figure 8) shows that practically no dolines and possible objects occur in the highest parts of the relief, where sediments associated with the Urucuia Group predominate. There is a higher density near the drainage incisions, especially accentuated from downstream to upstream towards the northeast. In the west-southwest block, lower sectors, there is a greater dispersion of occurrences, very much linked to the existence of shallows. This condition seems to be in line with the idea of a greater doline-forming process out at the points where the carbonates start to be exposed (Figure 8). This precisely indicates the capture of surface drainages for subterranean systems that develop in the area in a simultaneous and often not much interdependent (underlying karst).

When considering only the field dolines and the suspected and probable objects (249 in all), without including the shallows (inconclusive objects), we have a previous density of 5.03 dolines/km². Checking the shallows and possible dolines internal to compound depressions may further increase this value, indicating that this estimate is compatible or undersized with the present karstic area.

A density map considering this analysis, without shallows, was constructed (Figure 9a). A second map (Figure 9b) regarding only shallows and possible objects, both not yet confirmed (no dolines), was constructed and was used as a comparison. The density maps produced (Figures 8b and 9) show slightly different results because they considered different combinations of objects given as positive (Table 9). The first map (presented in Figure 8b) shows all positive points (456) and therefore presents maximum density values close to 46 depressions/km². Because of the greater number of shallows, there is a spread of lower densities over the area, especially in the southwest. The map that presents only the suspected (75% chance), probable (90%), and doline (100%) objects shows maximum concentrations of just over 38 dolines/km² (Figure 9a) and represent what the total area was calculated as an average density of 5.03 dolines/km². Finally, the last map (Figure 9b) shows the positive features without the inclusion of the field dolines. In comparison with the previous one, this reinforces the points for future checks, especially those of higher densification. This comparison shows areas where there is a greater doline forming process, or strong evidence of its occurrence vs areas where prospecting for dolines should be prioritized. The maps serve as indirect indications of the possibility of the occurrence of caves and greater karstification areas.

Figure 9

Different combinations of densities by objects. (a) The density generated without the shallows is observed. Therefore, only the field dolines and the most probable and suspicious objects, accounted for an overall average density of 5.03 dolines/km². (b) Field dolines were removed from the analysis, keeping the shallows, suspects, and probable objects. It, therefore, represents the potential for field prospecting of dolines and caves. UTM 23S, SIRGAS 2000.

4.1 Covered vs exposed karst

The particular doline typologies found in covered and exposed karst systems of the world and Brazil indicate the significance of tropical and humid climates in doline genesis. A pattern significantly influenced by the cover materials and drainage insertion, which is characteristic of the tropical environment, was observed in the typology of the dolines mapped in the field. A relation of doline geometries with hydrological features of the area is established. The region of the present work is characterized by the predominance of mobile (alluvium-glacial) covers over soluble rocks, influenced by the intense rainfall events. In the semi-arid northeast of Brazil, very different situations were observed, and doline typologies of dissolution and collapse were exposed [35]. Additionally, other areas present large polygonal-type features in the subtropical humid forest zone with annual rainfall up to 1,860 mm [36]. The predominance of suffosion and cockpit-type features can be linked with the dissolution dolines found in areas of exposed rock with little or no soil. This situation is quite distinct from the one found in a comparatively drier climate (576–698 mm), such as the Taurus Mountains in Turkey, where dissolution dolines predominate in exposed carbonate areas [37,38]. On the other hand, in an analysis of doline types in northeastern Italy, Calligaris et al. [39] identified the predominance of suffosion (29%), followed by collapse processes (rock 16%, cover 15%, and caprock 11%).

4.2 Sand erosion

The studied karst is a covered fluvio-karst type, and where the carbonate is found buried under a thin layer of alluvial sediments and the insertion of surface drainage occurs, processes of capture, subsidence and collapse through dolines occur more frequently. Few dolines in the studied areas preferentially occur in areas of contact between carbonates and siliciclastic sediments (Urucuia Group). The dolines may operate as catalysts in the emptying process of the sandy surfaces of mobile materials leading to a connection with the underground environment [1]. In this case, it is not about liquefaction processes of the cover material by flooding processes near channels or variation in the aquifer level, as some work on covered karst has pointed out [9,40]. It is the interconnection of surface erosion of the slopes accentuated by the gradient with subsurface emptying structures already well developed in the epikarst. The occurrence of cave systems and drainages driven by the lowering of the regional base level act coupled with the process of headward erosion, promoting the capture of surface sediments through suffosion and collapses. This explains the concentration of dolines near the flumes or paleochannels in the region (Figure 8), showing a concentrated and linear distribution in lowered zones and not dispersed on adjacent plateaus as observed in other works [9,40,41,42]. This can also be interpreted as prior underlying karstification, which, once exposed, allows for rapid capture of surface systems since the area fits as covered karst in transition to exposed [22].

Another interesting erosion and exposition of the already existent doline and cave-underground connectivity have been explained by Hussain et al. [18]. Some of the cover-collapse dolines were found filled with the sediments at a considerable degree of compaction. Underground cave streams during flooding days may lead to bottom erosion of filling material which may trigger the collapse of less cohesive material leading to geological hazards. Further details can be found in the study by Hussain et al. [18].

4.3 Area size and doline density

The adoption of smaller areas, however, can point to much higher values, such as the one carried out by Piló [43] in the Lagoa Santa karst region, reaching 18 dolines/km² using coarse-scale aerial photos (1:10,000). Results of LiDAR and UAV applications in small areas in Europe (e.g., Slovenia, Italy, Austria, and Hungary) show significant improvement in doline detection. The size of the study area and its impacts on increasing doline density were highlighted in the literature (Table 10). Bauer [28] studied a small (1.3 km²) isolated karst mountain in the Styrian Basin (Austria) using LiDAR data together with the dimensional filtering of the objects and found 109 dolines/km². Similarly, the work by Čeru et al. [42] on small portions of the Dobrava (0.76 km²), Poljšica (0.31 km²), and Podbrezje (0.38 km²) plateaus in the Ljubljana Basin (Slovenia) resulted in the highest doline density as presented in Table 10. These results over smaller karstified areas are compatible with what has been found in the present work, as doline density exceeded 45 dolines/km² in areas near Tarimba Cave (Figure 8b). On the other hand, Plan and Decker [44] analyzed an area of 59 km² in Austria equivalent to that of this work (50 km²) and found a high density of dolines (122 dolines/km²) delineated using an orthomosaic photo. Yechieli et al. [45] also found a high frequency of dolines (67 dolines/km²) associated with evaporites on the coast of the Dead Sea (Israel), in an area compatible with that of this study (60 km²). These discrepancies between high densities in studies that used different area sizes show the need for the adoption of area limits for density analysis, which should preferably consider the average frequencies rather than the maximum.

The average doline density (5 dolines/km²) in the studied area is a significant number for Brazilian karst [46] and it is five times more than the other karst in the region as of Lagoa Santa (Brazil), with one cave density of 1 doline/km² [47]. The density is well below the great (small area) majority of the karst regions studied around the world (Table 10). However, the density is compatible with that achieved in the region of the largest cave in the world (e.g., Mammoth Cave, Kentucky), where Wall et al. [8] performed detection with LiDAR data, achieving a density of 5.1 dolines/km². Equivalent results are observed in Ghizhou, China, with up to 5.06 dolines/km² [10], in part of Slovakia, with 6 dolines/km² [48], in Georgia and Florida in the United States, with 4.39 dolines/km² and 6.07 dolines/km² [49,50], respectively.

4.4 Doline density and cave proximity

Large concentrations of depressions are found near the caves in the region (Figure 8), connecting the largest underground system in the area (Tarimba Cave), with the surface. Some sites, such as in the headwaters of the Tarimba Cave, a higher density is also found without any cave occurrence, as inferred from the partial mapping found in the literature. This density indicates the higher speleological potentiality of these sites, which stresses the need for further detailed high-resolution mapping.

4.5 Environmental vulnerability and geohazard

The underground vulnerability of karst is a function of doline density (i.e., the higher the doline density, the higher the vulnerability and vice versa). This requires the addition of an extra GIS layer in the already existent vulnerability assessment models, especially the environmentally sensitive areas such as APANVR. A study on the same area by Ferreira et al. [2] highlights the need for such inclusion, as documented in the literature. The well-known karst vulnerability assessment model, concentartion, overlying layers and precipation, uses the O factor, which ranked the areas based on the soil thickness and doline presence. The small index values are assigned for poor soil with the presence of less or no dolines. In comparison, small to moderate ranking is used for the regions of deep soil (1 m) and numerous dolines and moderate to high ranking for the soil is mainly developed to a depth of 2 m and fewer dolines and very high, very deep soils (>2 m) and no occurrence of dolines. Detailed information can be found elsewhere [51,52]. This model requires some modifications to assess APANRV’s vulnerability better considering doline density and cave openings. In addition to that, high-resolution dolines mapping using UAV can help in geologic hazard assessment in the karst regions. The karst landforms create favorable conditions in the formation of depression/dolines through local denudation of the cover layer (deposits). Because of disturbed physical and mechanical properties, these depressions are dangerous and require their identification with the aid of UAV results. This mechanism has been explained in Section 4.2. The revealed regularities make it possible to carry out predictive assessments and identify hazardous processes. This can work as a reconnaissance in the application of geophysical techniques for depth and volume estimation of these sediments in the doline [18]. In short, detailed and well-presented hydrogeological conditions of the karst are required as a reconnaissance for further thorough investigations. Only then the proper land use planning of the areas on rational grounds is possible.