An Integrated Geophysical and Unmanned Aerial Systems Surveys for Multi-Sensory, Multi-Scale and Multi-Resolution Cave Detection: The Gravaglione Site (Canale di Pirro Polje, Apulia)

Abstract

Gravaglione represents one of the main swallow holes of the Canale di Pirro, low Murge, Apulia region, Italy. Here, after an intense rainstorm, a huge volume of rainwater accumulates at the surface. The drainage dynamics suggest that the Gravaglione could be part of a large, and potentially unknown, karst system. To verify this hypothesis and to acquire useful information on the possible karst environment features, an integrated aerial and geophysical multiscale and multimethod approach was applied. In particular, aerial photogrammetry, ground penetrating radar measurements and electrical resistivity tomography surveys were hence conducted and integrated to potentially detect the caves, define the subsurface volume possibly affected by karst systems and to verify the existence of links between the surficial morphology and the subsoil structure. The results provided interesting insights regarding the presence of a complex karst system extending up to 200 m b.g.l. and with a marked 3D nature. Overall, the Gravaglione case study demonstrates the geophysical approach validity and poses the basis for the development of an expeditive and low-cost high-resolution strategy for detecting and characterizing karst caves.

1. Introduction

Karst environments present intrinsic anisotropy given by the peculiar physical characters of the subsoil, so that karst water may flow in large quantities in conduits and caves, in concentrated and high-energy modalities, whilst only few meters away from the main active karst network, the rock material might be characterized by lack, or minor development, of underground voids and water flow [1,2,3]. In general, the unsaturated zone in karst comprises the soil (if there is any), the epikarst [4], and the vadose zone; in the unsaturated zone the water flows essentially in vertical direction driven by the gravity force. Below it, the saturated zone hosts the freshwater body, where water flows generally in a horizontal direction driven by the hydraulic gradient. In addition to the spatial anisotropy, the karst subsoil changes constantly in time [1,2,3,4,5,6,7,8]. CO₂ in flowing waters dissolves the carbonate rocks, and enlarges the initial fractures into larger conduits and caves. This is, in general, a slow process, but can also change rapidly (often within hours), due to strong response to hydrological events such as storm rainfall and snowmelt [9,10,11]. Further, the underground voids network may change in time, in consequence of variations in the conduit morphology (collapses, sediments filling, etc.), so that the saturated zone may transform into unsaturated, and vice versa [12]. The variable nature of karst terrains requires specific adaptation in the investigating techniques, to properly characterize such a peculiar environment. Although it is certainly necessary to collect field data about mineralogy, lithology, stratigraphy, geological structures and geomorphology of the investigated area [12], monitoring actions at high temporal resolutions are required, too. Detailed study about both underground and superficial karst features can be considered the right way to characterize the karst environment and may help in the choice of the best management actions. The comprehension of hydraulic functioning of dolines and endorheic basins give a significant contribution in the understanding of hydrogeological groundwater dynamics, while investigating underground karst voids close to the surface might help in the process of hazard reduction linked to collapse sinkholes [13,14,15].

In this framework, geophysical methods can play a key role for their ability in characterizing subsoil structure from the surfaces up to hundreds of meters, which configures them as a useful tool for cavity detection and mapping. This possibility is even more valuable if we compare geophysical methods with traditional methods of cave research, such as speleological exploration and surface surveys. While these last can be time-consuming, dangerous, and limited in scope, geophysical methods can provide a non-invasive and cost-effective means of detecting and mapping caves.

Generally, given a specific investigation target, the higher is the contrast of physical properties between the target itself and the surrounding material, the more efficient will be the geophysical methods in its detection. As a result, they can be particularly effective for karst environments characterization since voids, conduits and caves usually have a physical response largely different from the rest of the environment in which they develop. The most common geophysical techniques for karst environment characterization are ground-penetrating radar (GPR), electrical resistivity tomography (ERT), gravity, and seismic methods [16,17,18,19].

Here, we present a multi-method approach for karst environment characterization based on the integration of information gathered from unmanned aerial systems (UAS), GPR and ERT surveys which allows for a detailed subsoil structure reconstruction from the surface up to 200 m depth. Each of the three chosen methodologies will address a specific target. Since some karst features are visible at the surface, UAS will be used to verify their presence, and to have insights on how the surficial morphology is controlled by the characteristics of the karst system. GPR will furnish detailed information on the first 10 m of depth. Having a sub-metric vertical resolution, it could provide useful information to locate possible surficial access points to the karst system. Finally, ERT will depict the electrical properties distribution of the subsoil up to a few hundreds of meters allowing a low-resolution imaging of the karst environment. According to [20], karst aquifers are in some cases the only water sources available and generally provide water approximately to the 10% of global population [21]; therefore, it is of great importance to know their distribution and features in karst zones. In karst terrains, ERT (electrical resistivity tomography) can be hence useful not only to discover underground structures but also to obtain information about water pathways and the role of sinkholes in discharging water during the strongest rainfall events.

2. Canale di Pirro Overview: Geology, Geomorphology, Hydrogeology

Apulia (Southern Italy) is a region characterized by wide outcrops of carbonate rocks [22]. In general, the green colors in Figure 1 represent the limestone bedrock (Jurassic–Cretaceous), covered by the more recent carbonate calcarenites (Eocene–Miocene) and marine terraced deposits (Pleistocene). These latter are present mostly in Salento and at the coastal sectors of Murge, in central Apulia. Faults acted in modeling the topography, resulting in a horst and graben landscape, mainly arranged in NW–SE direction [23,24,25], with subordinate orientation of the landforms along NE–SW direction [26] (Figure 1). Due to the soluble nature of such geological formations, the regional landscape essentially developed through karst processes. Apulia can be subdivided into three main karst districts: Gargano, Murge and Salento [27]; they have been subjected to different geodynamic history with consequent differences in evolution of the karst processes.

Murge is the main karst area in central Apulia, extending several tens of kilometers inland from the Adriatic coastline. It can be subdivided into High (or NW) Murge, the upper portion of the plateau elongated in a NW–SE direction, with maximum elevation of 680 m a.s.l. at Mt. Caccia, and Low (or SE) Murge, the portion of the karst plateau closer to the sea, with a maximum elevation of 500 m a.s.l. [27,28,29]. They are commonly characterized by the absence of superficial water bodies, due to the nature of carbonate karst rocks, with rapid infiltration into the subsoil through fractures in the rock mass, or at specific infiltration points (dolines, swallow holes). In the whole area, abundant Terra rossa and residual karst deposits are present to fill the bottom of karst valleys and depressions such as polje, dolines and endorheic basins, making possible agricultural practices [30,31,32]. At these sites, due to the impermeable character of these deposits, accumulation of water may occur during and after rainstorms, and is at the origin of the main events of flash floods [33,34,35,36].

The Canale di Pirro polje is located in the Low Murge [31,37]. It can be described as a WE-oriented tectono-karst valley with an overall length about 12 km, with flat bottom and tectonically controlled ridges at its N and S boundaries. It extends from the towns of Gioia del Colle and Putignano toward the east [27,37], until reaching the Murge escarpment, that is the main NW–SE tectonic line separating inland Murge from the Brindisi plain towards the Adriatic coastline [38].

The structural characterization of carbonate rock mass discontinuities in the polje (summarized in Figure 2 as contour charts in stereographic projections) indicate a NE–SW general trend.

At the polje and in its surrounding, the most abundant superficial karst features are dolines and endorheic basins [31]. Their orientation is generally controlled by the mainly structural setting, along the NW–SE direction, which is in agreement with the tectonic control of many karst landforms in Apulia. The most common shapes are elliptical and sub-circular; this is probably because dolines, characterized by more regular shapes, and minor size as well, are more numerous than endorheic basins [31].

Under the surface, many caves open within Canale di Pirro; even though typically of limited size, and often used as stables for livestock, the longest and deepest cave of the region (Grave Rotolo, [40,41,42]) opens at the polje bottom. The karst geomorphological data collected at the Canale di Pirro polje have been summarized in the map shown in Figure 3 and they are described in detail by Pisano and co-authors [31].

As typical of poljes in other karst regions of the world [43,44,45,46,47,48,49], the Canale di Pirro polje often presents flooded sectors, becoming temporary lakes after severe storms [46,50], due to the hydraulic functioning of the swallow holes at its bottom. Generally, swallow holes are clogged by Terra rossa deposits; the slow infiltration of water throughout the residual deposits promotes the opening of ways for the water, so that, after it remains at the surface for hours/days, later it can be drained in short times. The Canale di Pirro bottom hosts two main swallow holes: Rotolo cave (PU 355) and Gravaglione (PU 354). Another minor swallet, corresponding to Grotticella del Canale di Pirro (PU 1509), was destroyed years ago to create a vineyard on the southern side of the polje, and its precise location has been lost [40].

Gravaglione is an elongated doline and the most important adsorption point of the polje; after an intense rainstorm a huge volume of rainwater accumulates at the surface, flooding the surrounding fields (Figure 4). The water needs time to be absorbed by the swallow hole; in some cases, the area results drained for a few hours, other times it remains flooded for at least 48 h, as in the case of a rainstorm that occurred in December 2020 (Figure 4).

The other swallow hole, about 500 m downvalley of Gravaglione, is the access to Grave Rotolo: originally a small swallet in the field, after the explorations in 2012 it has become a 324 m deep karst system (the deepest in Apulia), reaching the groundwater at about 260 m from the ground [40]; the cave is still under explorations and many cavers and scuba-divers are being involved in this activity.

3. The UAS and Geophysical Integrated Approach

The geophysical characterization of the area surrounding the Gravaglione swallow hole was realized by combining and integrating three different survey techniques: UAS, GPR and ERT. In the characterization of karst environments, the application of these investigation techniques does not represent itself a novelty. All methods, in fact, have been successfully applied for cave detection or, more in general, to gather useful information of the surface, subsurface and on the dynamics of the processes therein taking place. Even though they can be standalone applied, these methods perform better when integrated. This is particularly true when, due to the environmental condition, the result of a single method application is not uniquely interpretable (e.g., Ref. [51] and references therein).

Based on this, UAS, GPR and ERT survey results will be firstly individually presented and then integrated to create a common interpretational framework. The surveys presentation order will reflect our intention of providing a characterization “from the surface to the deep” of the investigated area. Accordingly, UAS results will be firstly presented, followed by GPR and ERT. The different surveys were performed on a two-year span of time and, at the time of acquisition, they were not planned as part of an integrated survey. For this reason, even covering and investigating the same area (the Gravaglione), they do not exactly overlap in terms of total area covered and of profile traces (for the GPR and the ERT) (Figure 5). Most of the data were collected in December 2021 and in February 2023, in periods of relatively dried ground conditions. The tomographies that will be shown in paragraphs 5.3.1 and 5.3.2 were collected in July 2021. In all the cases, no relevant rainfall events occurred in the week prior the fieldwork operations. Assuming similar hydrological conditions both in December 2021 and February 2023, we do not expect relevant changes in the physical properties distribution within the subsoil. This guarantees the integrability of the survey results.

3.1. UAS Mapping

Surface characterization is nowadays routinely conducted by creating three-dimensional models from passive sensors, such as photographs acquired with UAS based on the technique of photogrammetry. Photogrammetry makes it possible to create three-dimensional products (e.g., point cloud, DSM—Digital Surface Model, and 3D model) of a real scene, using a technique called SfM (Structure from Motion [52,53,54,55,56]). In the field of close-range remote sensing (RS), it is a widely used technique because it is economical and guarantees excellent results. An optimal approach to photogrammetric surveying is to combine the drone survey with GNSS acquisition so that the data can be aligned with great precision, especially on the height (z) coordinate [57,58]. In addition, this method relies on precise mathematical and geometric rules (projection) to extract the coordinates of homologous pixels in the images. Each image, in fact, contains all of the information necessary to reproduce the geometry (point cloud and mesh) and color (texture) of the photographed object [58,59]. Photogrammetry aims to obtain metrically accurate, precise, and geometrically reliable information based on the principle of collinearity [60], correlating the coordinates in the sensor plane (2D) with those of the three-dimensional object (3D) [60,61].

3.2. GPR Investigations

Ground penetrating radar (GPR) is widely applied for geological issues thanks to the high resolution offered by the method despite the low depth of investigation, generally limited to a few meters. Indeed, GPR was successfully applied for characterizing and identifying geological faults [62,63,64] or for hydrogeological purposes [65,66,67,68]. Furthermore, the high sensitivity offered by GPR for the identification of voids and/or caves makes the method particularly efficient for detecting and identifying subsurface karst features [69,70,71,72,73,74,75]. The method is based on the introduction of electromagnetic (e-m) signals that in presence of variations in the three physical properties of dielectric permittivity, electrical conductivity and magnetic permeability are subject to phenomena of attenuation and reflection. The investigation depth varies from centimeters to a few to dozens of meters, depending on the antenna center frequency (2 GHz–10 MHz). For geological studies, the operating frequencies are generally lower than 200 MHz that provide the possibility to have sub-metric resolution at depths greater than 2–3 m. Considering the reduced depth obtainable with the GPR and the noise affecting the data when low frequency antennae, generally not shielded, are adopted [62,63,64], it is fundamental to use an integrated and compared approach which involves the use of other non-invasive (i.e., seismic analyses, electrical tomography) or invasive methodologies (borehole data) [76,77,78].

3.3. Geoelectrical Analyses

The geoelectrical method proved to be successful in subsoil characterization in a wide range of frameworks. When applied in the karst, its effectiveness is based on electrical conductivity contrast between the hosting material and the voids, allowing cave detection, reconstruction of their geometries and nature of the filling materials if present. As a result, several successful examples of electrical resistivity tomography (ERT) applications for revealing karst aquifer structure and cave detection are available in the literature (e.g., [18,79,80,81,82,83]). The interpretation of electrical resistivity can, however, be challenging especially if no external constrains are applied in the interpretation of the geoelectrical results. Cave, conduits, and other elements of karst systems may have very different electrical behaviors depending on the filling material. If air-filled, the karst system is generally imaged as a highly resistive feature. In case of water or clay-filled voids, it is imaged as a low resistive feature, with the additional difficulty in distinguishing if the low resistive areas are associated with hydraulic conductivity (water presence) or not. Additionally, other site-specific condition may decrease the magnitude of the resistivity contrast between the karst system and the hosting rocks (e.g., [84,85,86,87]).

4. Data Acquisition and Processing

4.1. UAS

Three-dimensional acquisition of the study area was carried out using a DJI Mavic 2 Pro drone, with a 1-inch Hasselblad 20 mega-pixel sensor. The flight mission was performed automatically using the UGCS software (v.4.6.520). The software allowed to create a flight plan by selecting the ‘Photogrammetry’ mode and setting the altitude and flight speed, camera options, capture mode, lateral and frontal overlap between shots [67]. In addition, thanks to the terrain awareness function of the UGCS, the flight was carried out using a DSM (Digital Surface Model), useful for the drone to fly at a constant altitude above the ground and to have a constant GSD (Ground Sample Distance) throughout the survey. The DSM used was provided by Tinitaly (http://tinitaly.pi.ingv.it/, accessed on 4 April 2023) [88,89].

The survey was conducted with these specifications: (i) altitude 70 m AGL (Above Ground Level); (ii) 4 m/s in safe shooting mode, i.e., the drone stops to take the photo in order to avoid the rolling shutter effect; (iii) nadiral camera (90°); and (iv) overlap between images: 75% frontal, 70% lateral. A high-accuracy GNSS (Trimble R2), linked with national GNSS networks, was used to acquire ground targets as ground control points (GCPs). Twenty targets were positioned and acquired throughout the area of interest. Of these, 15 were used as control points and five as verification points.

The drone flight captured 624 images. The photogrammetric process was carried out entirely within the Agisoft Metashape^® (Agisoft LLC, St. Petersburg, Russia) software (v.2.n) following these steps: (i) dataset creation; (ii) camera calibration; (iii) photo alignment and point cloud creation (110,829 tie points); (iv) identification and addition of coordinates to GCPs; (v) data correction and alignment adjustment (native correction in Agisoft); (vi) creation of a dense cloud (13,382,622 points); (vii) creation of DSM; and (viii) creation of orthophotos.

4.2. GPR

The Gravaglione area was studied performing several GPR profiles at the frequency of 40 MHz (Radarteam Sweden AB, Gränsvägen, Sweden) according to the scheme depicted in Figure 5. Data were acquired with the GSSI SIR-3000 GPR System and markers were placed every 5 m for assigning the right coordinates to each recorded trace.

The use of a low frequency antenna (40 MHz) was preferred to other possible low-frequency antennas (with a central frequency up to 200 MHz) to achieve a reasonable high investigation depth along the whole investigated transect lines. Limitation in the GPR signal penetration, in fact, must be expected considering the fieldwork activities periods (December and February) and the Gravaglione geological framework, where the Terra rossa cover may have a thickness that in some areas exceed 5 m.

In this case study, the observation time window of each radargram is set to 400 ns and discretized by 1024 time-samples. The acquisitions were performed with the use of markers placed every 5 m. Furthermore, on the obtained radargrams we performed the zero timing by setting $t_0=40 ns$. Then, the BKG removal procedure was applied to remove this constant signal over the entire time window and a DeWOW filter was used for enhancing events greater than 15 ns. A bandpass filter was employed to reduce the effects on signal characterized by frequencies far from the operating one. To further reduce the noise, three traces at a time were stacked. Then, we performed the topographic correction of the data and the two-way-time-depth conversion by fixing the relative dielectric permittivity value. In the Gravaglione geological framework, it is almost impossible to define a dielectric constant value which can be considered representative of the whole geological volume. Carbonates, Terra rossa and air (assuming that possible air-filled cave are present in the area) strongly differs, in fact, in terms of dielectric constant values. We hence assumed a relative dielectric constant of 10 ($\epsilon =10)$ by considering that the prevalent geological unit in the area is composed by carbonate rocks. Furthermore, this assumption is validated by the comparison of the GPR migrated data that have highlighted the most interesting anomalies at similar depths to those detected by ERT.

In this scenario, characterized by a v_em of 0.1 mns⁻¹, the expected vertical resolution provided by the antenna is 0.5 m (f = 40 MHz), approximately. Finally, to better visualize the encountered anomalies and make easier the comparison with ERT data, the envelope of each trace was calculated.

In detail, data processing performed with the support of the commercial software Reflex-W (Dr. Sandmeier, Karlsruhe, Germany), provided few steps to avoid artifacts according to the processing chain shown in Figure 6, which involves the following procedures:

• -. Data editing, for assigning the real coordinates to each trace of the acquired dataset;

• -. Zero timing, which defines the actual starting time, $t_0$, of the observation time-window;

• -. Frequency filter applying

  • DeWOW procedure, which subtracts at each collected waveforms its mean value along the time axis;

  • Background removal: this procedure helps removing or mitigating the signal contributions due to antenna coupling, air-material interface and (undesired) horizontal reflectors;

  • Butterworth filter, that allows for simultaneous performing of lowpass and highpass filters, in our case the frequency content between 20 and 70 MHz was considered.

• -. Stacking of the traces for lowering the data noise: a number of five traces was considered;

• -. Topographic correction to insert the real altitude of each acquired trace and data migration. This step was efficiently supported by the high resolution DSM recorded with UAV: for each profile a vertical section was extracted by the DSM using the software Global Mapper and the data were used for correcting the topography. As regards data migration, in addition to information from the literature, the ERT results have been also considered;

• -. Hilbert transform that provides an equivalent transformation between time domain and frequency domain by Fourier transform [90] for improving the detection of e-m anomalies. This step was adopted to make it easier to compare the GPR data with ERT data.

4.3. ERT

Geoelectrical survey was based on the acquisition of Dipole–Dipole (DD, both in direct and reverse configuration) and Wenner–Schlumberger (WS) data along each of the transect lines reported in Figure 5. DD dataset includes reciprocal measurements for data error analysis (e.g., [91]) and the maximum spacing between current and potential dipoles was set to “6*a” (where “a” is the dipole length). A Syscal pro 48 channels (Iris Instruments) connected to multielectrode cables was used for performing all surveys.

A total of 11 ERTs were acquired, as follows: ERTs T1–T9 by adopting a minimum electrode spacing of 5 m; ERT M1 and M2 by adopting a minimum electrode spacing of, respectively, 10 m and 20 m. ERT T1–T8, M1 and M2 were planned and positioned assuming a possible karst system with a main development corresponding to the Gravaglione location. ERT T1–T8 were realized aiming to define the possible lateral boundaries of the karst system, while ERT M1 and M2 were acquired to investigate its possible vertical extension. ERT T9 was acquired lately on, to verify some of the ERT M2 results. All the geoelectrical datasets were filtered eliminating anomalous apparent resistivity values and by using ResIPy [92].

Geoelectrical results will be presented following a three-level scheme.

The first level consists in the identification of the best geoelectrical array for cave detection. From experimental cavity studies and numerical simulations, Wenner–Schlumberger (WS) and Dipole–Dipole (DD) configurations proved to be the most effective. Doyoro et al. find out that DD, despite its low signal to noise ratio, provided the highest model resolution that shows relatively distinct anomaly geometries [93]. WS is also effective, even less than DD, but its results are more challenging when it is necessary to determine the correct anomaly geometries. Zhou et al. show that the use of a mixed array, obtained by combining and processing together apparent resistivities from different arrays, is the most effective for karst characterization, but implies a significant increase in data collection time and effort [94]. All of these results converge hence on indicating DD as the best configuration to use for mapping the subsoil in karst areas, but they do not consider the possibility that the geoelectrical arrays performances also depend on the environments in which they are applied and on the existing geophysical contrast between the cave and the surrounding environment. For an effective data analysis and interpretation, a comparison between geoelectrical models obtained along profile T1 by inverting apparent resistivities from different arrays (DD, WS and Mixed array) will be performed and the most effective configuration for Gravaglione area will be identified.

The second level is realized by integrating different scales of geoelectrical data to create a multi-scale and multi-resolution model of the Gravaglione area. The multi-scale survey is obtained by combining geoelectrical data with different resolution and investigation depth, acquired along a common transect. According to [95], this approach results in an increase in the target resolution also at depth, which is made possible by the greater number of receivers that sense the resistivity changes in a fixed cell, thus better constraining the inversion process. As this approach can be applied only to geoelectrical data, with different resolutions and investigation depths, acquired along a common transect, only ERT T1, M1 and M2 dataset were used.

The third level consists in the realization of a quasi-3D geoelectrical model of the area. A 2D geoelectrical image, in fact, cannot be sufficient for characterizing complex subsoils like those interested by karst phenomena, which are rarely 2D in nature [96]. This is especially true if a priori knowledge of a possible preferential development direction of the survey target is missing and the two 2D geoelectrical profiles are not correctly positioned. It is hence desirable, whenever possible, to prefer 3D ERT surveys to 2D ERT [97]. In the Gravaglione area, given the instrumental impossibility to acquire real 3D geoelectrical data with a sufficient spatial data coverage, the approach proposed by Cheng et al. was preferred, which consists of the use of a quasi-3D approach based on the acquisition of multiple 2D surveys combined in a 3D analysis [18].

5. Results

5.1. UAS

The acquisition and photogrammetric processing of the UAS data resulted in the output of (i) high resolution RGB orthophotos and (ii) DSM, both georeferenced (Figure 7).

The use of GCPs resulted in data with an estimated average error (due to ground validation points) of approximately 2 cm (East and North) and 4 cm (Altitude). The resolution of the DSM was approximately 14 cm/pixel, while that of the orthophoto was 1.8 cm/pixel. The orthophoto in the visible spectrum made it possible to appreciate from above the entire area covered with other tools. The DSM, on the other hand, showing features of geo-morphological interest, made it possible to appreciate the trend of features and micro-features in the study area. In particular, it was possible to observe the traces left by (i) depressions and concavities on the ground, possibly linked to cavities in the subsoil, and (ii) preferential pathways of the water passage. In particular, the latter data were then used for the integrated multi-sensor analysis.

5.2. GPR

GPR radargrams provided interesting information about the presence of potential caves in the subsoil and presence of fractured rocks. Assuming valid the adopted dielectric constant value of 10, the following images can provide a useful indication of the GPR anomalies’ depth. Figure 8 shows the results obtained with the two perpendicular profiles F2 and F6 (for location, see Figure 8a). The first transect allows for the identification of the presence of two anomalies placed at distances ranging between 25 and 40 (C1) and 110 and 170 (C2), respectively. The top of both anomalies is placed at depth of 3.5–4 m. The anomalies seem to be characterized by a time-dimension of 50 ns that, if we assume e-m propagation in the void (v_em = 0.3 mns⁻¹) would correspond to a depth of 5 m. The second radargram (F6) intercepts the anomaly C1 and, then, close to the Gravaglione detects two separated anomalies, C3 and C4, presumably linked to the presence of fractured or altered rocks or voids. In this case, the two anomalies are more superficial, and their top is at depth of 3 m, approximately.

Similar results are obtained for the radargrams F3 and F5 acquired at the eastern edge of the site (see Figure 9). The radargrams permitted to presumably identify the presence of zones characterized by fractured rocks, from a depth of two meters. In F3 two anomalies are mainly detectable, C5 and C6, located between 70 and 95 m, and between 105 and 115 m from the starting point of the acquisition, respectively. Given their proximity, it can be hypothesized that they might be linked. F5 shows no reflection in the first 100 m of acquisition: the only two anomalies are at the end of the radargrams, at distances from 100 to 125 for C2, which is the same anomaly already detected by the radargram F2 and between 150, and from 130 up to the end of the radargram for C7.

Figure 10 shows a 3D visualization of the radargrams acquired in the Gravaglione area. The strongest reflections are imputable to the presence of voids and fractured rocks that appear to develop near the swallow hole, as well as close to the south-eastern edge of the investigated area (see Figure 10a,b). The highest amplitude reflection areas, projected in plan (Figure 10c), allow for highlighting the zones characterized by a higher probability to encounter cavities or fractured rocks. All the results are addressed to identify the presence of a structure aligned in NW–SE direction.

5.3. ERT

5.3.1. Assessing the Most Effective Geoelectrical Array

To assess what is the most effective geoelectrical array for characterizing the Gravaglione area, we produced geoelectrical models by inverting DD and WS datasets independently, and then jointly by creating the mixed array. The analysis focused on the T1 profile, the closest to the Gravaglione swallow hole. For this survey, quality of both DD and WS data were good. The acquired DD dataset was composed of 906 measuring points among which, due to imperfect installation of one electrode, 80 measurements placed in the western portion of the profile were discarded. The remaining dataset (826 measuring points) was characterized by reciprocal errors lower than 1.5%. The imperfect installation of the electrode also affected the WS dataset. Once the bad data were removed, the WS dataset passed from 565 measuring points to 531, with stacking errors lower than 2%.

Figure 11 shows the comparison between geoelectrical models related to the separate inversion of WS and DD dataset (upper and middle panel, respectively) and of the mixed array (WS + DD) (lower panel). The three models are presented by using the same color scale for allowing an easier comparison.

All of the models image a surficial conductive layer (ρ < 100 Ω*m), whose thickness varies along the profile, reaching its maximum value in correspondence of the eastward portion of the section, and that can be associated to the topsoil cover mainly composed of Terra rossa. Below it, the resistivity increases in a not uniform way, imaging the carbonate rocks. Under the Gravaglione area, a high resistivity zone is present (marked with capital letter A in Figure 11), that could be associated to presence of a cave, or to a highly fractured area. Carbonate rocks may be generally characterized by higher resistivity values when compared to those retrieved by the Gravaglione geoelectrical surveys. We retain that the presence of Terra rossa, with a high content of iron oxides, could explain the observed resistivity values if we assume filling, even-partial, of the carbonate pore system.

As clearly visible, WS (see Figure 11a) is the less appropriate array configuration for describing the Gravaglione karst system, while both DD (Figure 11b) and Mixed array (Figure 11c) image what could be the electrical trace of a possible cave, between 10 and 30 m of depth. In particular, at least in the geological context of the Gravaglione, DD performs much better, and this justifies the choice to show, in the following paragraphs, only the DD results.

5.3.2. Multi-Scale and Multi-Resolution Characterization

To better image the deeper structure of the Gravaglione, other two geoelectrical surveys, M1 and M2 (see Figure 5 for location), were performed along the same trace of T1 with a minimum electrode spacing of 10 m (M1) and 20 m (M2), respectively. Data were acquired by using both DD and WS configurations but, considering the results previously discussed, only DD data will be presented.

Both M1 and M2 datasets were filtered eliminating all of the data characterized by a reciprocal error higher than 15%. This operation filtered out the 12% of the M1 dataset, and the 6% of the M2 dataset, respectively.

In Figure 12, a comparison between T1, M1 and M2 inversion results is presented. In the figure, the distances on the x-axis in the T1 and M1 panels (upper and central panels) were shifted in order to have the position x = 0 coincident with the starting point of the M2 survey (lower panel). The surficial part of the three models is basically similar and depicts the presence of a topsoil cover with variable thickness along the profiles. An exception is observable in the eastern section of the largest scale tomography (M2) where a highly resistive area is present (capital letter B in Figure 12). As regards to the resistive area (capital letter A in Figure 11 and Figure 12), previously associated to the karst system (or more in general, to a highly fractured area), it appears not to be confined in the first 30 m of depth as it was imaged in the small-scale tomography (T1). In M1 (central panel in Figure 12), it extends and deepens eastward down to a depth of 50 m. At greater depth and farther to the east, enhanced resistivity values are still visible, but they seem not higher enough to be indicative of the presence of a karst system.

Again, as for T1, the deeper portion of the model M1 is the less resolved area, and that with the lowest model sensitivity. By enlarging the investigation scale, M2 images a continuous U-shaped resistive body, with its western boundaries in correspondence of the Gravaglione area, where it includes the anomaly A; it develops eastward, reaching a maximum depth at about 100 m below the ground surface, and emerges then’ toward the surface, about 650 m along the profile. At the surface, in correspondence with the right branch of the U-shaped resistivity anomaly, no morphological depression of the ground is visible. Farther toward east, from 700 m to 800 m along the profile, a shallow resistive layer is evident (capital letter B).

To investigate with higher resolution the last two hundreds of meters along profile M2, an additional tomography, named T9 (see Figure 5 for location) and with electrode spacing of 5 m, was realized along the M2 trace, centered at about 700 m (dashed purple line in the M2 panel of Figure 11). No data were discarded due to filtering operation with a mean reciprocal error of the dataset lower than 1%.

The T9 geoelectrical model is shown in Figure 13, confirming the presence of a high resistive surficial area (capital letter B) in the east sector of the tomography, but not showing a marked vertical resistivity anomaly in correspondence of the right branch of the U-shaped anomaly visible in Figure 12.

Considering the coherence between the four geoelectrical models (T1, M1, M2 and T9) and the spatial overlapping of the surveys, an additional model was obtained by following the approach proposed by [94] and realizing a multi-scale survey.

We hence integrated T1, T9, M1 and M2 data. As a result, a dataset with a non-uniform data distribution (Figure 14, upper panel) was formed and inverted for obtaining the resistivity model in Figure 14 (lower panel). The model retains the same features of the four single models (Figure 11, Figure 12 and Figure 13) in a unique geoelectrical image, adding to the large scale M2 survey the surficial and intermediate-depth details, respectively coming from T1, T9 and M1. In this way, the U-shaped resistive body is better resolved and imaged. In particular, the integration of the higher resolution tomographies (T1 and T9) adds details in the surficial part of the models, allowing for a more precise interpretation of the geoelectrical results in terms of the possible presence of karst systems, and for a better definition of their characteristics.

The resistive area (ρ > 2000 Ω*m), whose shallower portion is located just below the depressed area corresponding to the Gravaglione swallow hole, can represent hence the electrical signature of a large karst system mainly developing toward East, and that reaches a maximum depth around 150 m b.s.l (below surface level). Moreover, outside this area, the geoelectrical model depicts a subsoil in which carbonate rocks are characterized by a not uniform resistivity value distribution, with generally not so high resistivity values due to the effect of the Terra rossa deposits. Considering that in Grave Rotolo the groundwater is at about 260 m b.g.s. (well below the maximum reached investigation depth), this spatial variability of the carbonate rocks resistivity could be likely linked to a different fracturing level and filling materials of the rock matrix.

5.3.3. Quasi-3D Approach for Karst Environment Characterization

What presented above testifies the successful use of the geoelectrical method in identifying the possible area interested by presence of a karst system. However, a 2D geoelectrical images cannot be sufficient for characterizing complex subsoils affected by karst phenomena.

To this purpose, the 3D inversion program R3t [98] was used, that was included in ResIPy, which is based on an Occam’s solution to the 3D inverse problem, and employed an unstructured tetrahedral finite element mesh for forward modeling and resistivity model parameterization. The 3D modeling was limited to the area covered by the ERT surveys T1-T8 for having a good spatial data coverage and to limit the amount of computational resources required. A 3D unstructured finite element mesh consisting of approximately 310,000 elements was therefore created to combine the eight ERT lines, for a total of 384 electrodes. To save computational resources, the topography of the area was considered flat. The combined ERT dataset consists of 7840 measurements (with reciprocal data for error analysis) in DD configuration.

Figure 15 presents some of the 3D modeling results. The uppermost slice is a DEM model extracted from the UAS survey whose extension was limited to the area covered by the 3D modeling. Resistivity slices at different depths are shown starting from 300 m a.s.l. down to 260 m a.s.l. The shallower part of the model is clearly affected by non-uniform data distribution (Figure 15). In the easternmost and westernmost sections of the model, where data density is lower, the distribution of the anomalies mainly follows the profile traces. The same does not occur in the central portion, where data coverage is high and, in general, in the deeper sections of the model. As for these latter, the resistivity distribution appears to be more regular for the combined effect of a less complex geological setting (in terms of presence of geological units) and a coarser model mesh. By comparing slices of the 3D model corresponding to the 2D section (comparison not shown here), the 3D modeling depicts similar resistivity distributions, but with a more compressed resistivity range showing a maximum observed resistivity just above 1000 Ω*m (it was about 3000 Ω*m in the 2D sections).

From the surface up to z = 300 m a.s.l., the more conductive zones can be associated with the presence of a topsoil cover (Terra rossa) whose thickness, as also highlighted by the 2D models, varies in a non-uniform manner in the investigated area. It is worth noting that topsoil thickness shows a good spatial correlation with the morphologically depressed area (Figure 15). The more depressed the surface (e.g., in correspondence of Gravaglione), the thicker the topsoil cover is.

Increasing the investigation depth, the depth slices are dominated by the presence of a resistive volume which extends from z = 300 m a.s.l. down to 260 m a.s.l. This volume may represent the 3D electrical signature of the karst system that was also imaged by the previously shown 2D surveys.

A more effective visualization of the possible extension of the karst system in the Gravaglione area, probably linked to that partly described in previous works few about half a km downvalley [99], may be conducted by extracting from the 3D model the isovolumes related to zones characterized by a resistivity equal or higher of a given value. This operation is performed in Figure 16, where two possible isovolumes were calculated. The dark red (Figure 16a) is the more conservative in terms of spatial extension, and only includes the 3D model portion with resistivity equal or higher to 1000 Ω*m. This isovolume should cover, in the limit of the ERT method and model resolution, the minimum extension of the karst system. The gray isovolume includes the areas with a resistivity equal or higher than 500 Ω*m (Figure 16b); it is more extended and may also include areas characterized by a higher fracture level, but not necessarily interested by presence of voids and caves. Note that the limit of this isovolume was arbitrarily chosen, as it only represents the resistivity value higher than that observed in other portion of the 3D models where the carbonate rocks seem to be not interested by the possible presence of any karst system, and that have hence a lower resistivity value.

6. Discussion

The combined use of different geophysical investigation techniques allowed a detailed characterization of the Gravaglione area from the surface down to 200 m b.g.l.

GPR and ERT results pointed out how the subsoil is not homogeneous in terms of physical properties imaged by the methods. This inhomogeneity has been interpreted as indicative of the presence of both highly fractured areas within the subsoil and/or voids and caves of different sizes. The best confirmation of these conclusions could derive from the GPR and ERT data integration, or from direct comparison. Unfortunately, in the case study presented here, this operation could not be systematically performed, due to the only partial overlapping of the surveys and to their different resolution and investigation depth. The GPR, in fact, has a high resolution but its depth of investigation is limited to the first 8–10 m of depth b.g.l.; on the other hand, ERT has a much lower resolution (at least in the used configuration), but it reaches far more relevant depth of investigation.

It derives that a direct comparison between the two methodology results could be limited only to the shallow portion of the ERT sections, with the possibility that the discontinuities imaged by the GPR do not have a counterpart in the low resolution geoelectrical image.

Figure 17 shows the comparison between the GPR line F2 (yellow arrow in Figure 17a) and the near ERT T1 (blue arrow in Figure 17a). In particular, the radargram (Figure 17b) is compared with the most superficial part of the ERT1. The comparison result allows us to enhance the significance of the data and to demonstrate the validity of the assumption made on the relative dielectric constant ($\epsilon =10)$ performed in Section 4.2. Indeed, it is possible to note the presence of some zones characterized by a strong attenuation in the radargrams (a1, a2 and a3) which correspond to the most conductive strata identified by the ERT1. Likewise, the stronger reflections in the radargrams (C1 and C2) correspond to the more resistive nuclei of the ERT. However, if the reflections recorded with GPR seem to suggest the presence of likely caves near the surface (depth < 5 m), the resistivity values recorded in the ERT for C1 and C2 allow us to hypothesize the presence of areas characterized by highly fractured carbonate. For this reason, the results show the importance of GPR and ERT data integration to improve the interpretation of the geophysical results, and to potentially distinguish between geological scenarios with presence of fractured rocks and/or caves.

Another interesting perspective for data analysis and integration between GPR and ERT is offered by the processing of the GPR data with the Hilbert transform (Figure 17c) which allows for better highlighting of the discontinuities in the subsoil by the analysis of instantaneous amplitude of the processed signal. Where abrupt variations in the energy distribution within the recorded traces occur, the instantaneous amplitude significantly increases; on the other hand, where the attenuation of the signal is higher, more conductive conditions are encountered in the subsoil. The powerful advantage of the Hilbert transform consists of its capability to reproduce a “simplified image” of the subsoil, more easily comparable with ERT results (if compared to classical processing of the GPR data).

A further alternative strategy to integrate GPR and ERT results is based on verifying if some spatial correlations exist between the occurrence of GPR anomalies (see Figure 10) and the ERT ones (isovolumes in Figure 16). From this point of view, the UAS survey has a twofold role: first, it provides a detailed “base” on which to project the other methodologies results; secondly, and most importantly, it allows for checking if the presence of geophysical anomalies affects, or has a signature, also on the surficial morphology.

The UAS DSM was, hence, re-plotted limiting the elevation range between 308.7 m and 309.6 m (Figure 17a) to highlight the presence of small morphological variations, possibly masked in Figure 6 by the use of an elevation range which includes the maximum and the minimum elevation of the area. From the DSM, the isoline related to the following elevation values were thus extracted (Figure 17b): 308.7 m, 309.0 m, 309.3 m, 309.6 m.

It is possible to see that both the GPR and ERT anomalies show a striking correlation with the ground morphological variations. GPR anomalies almost never fall in the more depressed areas, but are mostly aligned to small morphological steps, or placed across these latter (Figure 18a).

ERT derived isovolumes (the same presented in Figure 16 but with a coloring scheme, which reflects the depth of the isovolume areas) seem to drive the surficial morphology with their most surficial portions (warm colors) enveloped by the 309.3 m DEM contour line (Figure 18b).

Further, GPR and ERT anomalies are well correlated. The most elevated portions of the ERT isovolumes always occur where GRP anomalies are visible. Where a correlation is not visible, it is possible that (i) the ERT anomaly is deeper than the maximum GPR investigation depth, or (ii) the ERT survey does not reveal any presence of anomaly due to the small anomaly size compared to the ERT resolution. Furthermore, it must be considered that while the isovolumes are the results of a real 3D modeling, the GPR anomaly maps was derived by the 2D analysis, and present, therefore, only the anomalies located along the acquired profiles.

7. Conclusions and Future Perspectives

Karst environments are complex systems characterized by an intrinsic anisotropy given by the peculiar physical characteristics of the subsoil in which they develop. Their characterization is fundamental for understanding the hydrogeological groundwater dynamics, but it might also help in the process of hazard evaluation, linked to collapse sinkholes. It is hence clear the necessity to perform detailed studies about both the underground and superficial karst features, which may help to take the best management actions in terms of mitigation of the risk. To this aim, geophysical methods proved to be effective and valuable tools for karst environment characterization, provided that the framework conditions (in terms of geophysical contrast between the karst system and the environment in which it develops) are favorable, as in the case of the Gravaglione area. There, the geological and morphological conditions, combined with the knowledge of karst features in the area, reduced the degree of ambiguity associated with the interpretation of the geophysical data, thus improving the related outcomes.

As a result, the integrated UAS and geophysical approach adopted successfully imaged the Gravaglione area subsurface structure from the ground level up to 200 m b.g.l., revealing the possible presence of a wide karst system whose presence also affects the surficial morphology. The karst system has its most surficial part in correspondence of the Gravaglione swallow hole and, from there, it deepens and develops towards the East.

Most of the results were derived by the ERT survey due to its investigation depth compatible with the characteristics of the karst system. GPR and UAS, however, demonstrate their capability for characterizing the most superficial part of the investigated karst environment. The obtained results allow us in this way to constrain and improve the ERT results because both the methods are more cost- and time-effective when compared to the ERT, and offer a higher resolution. These are not secondary aspects when it is necessary to scan large areas and where one of the survey aims is to locate a possible explorable access to the karst system.

This could be the case, for example, of the eastern sector of the investigated area where GPR lines show strong shallow reflective areas, and where the only performed ERT highlights the presence of a surficial resistive area.

The investigation at Gravaglione and surrounding areas is hence an ongoing process that did not end with the results here presented. Additional, higher resolution investigations will be performed with the twofold purpose of determining the real extension of the Gravaglione karst system, and of investigating its possible connection with the Grave Rotolo system as well.

From a methodological standpoint, the obtained results demonstrate the validity of the three-levels investigation strategy adopted for defining karst environment features by using geoelectrical methods and GPR and UAS survey.

With regards to GPR, the adopted frequency (i.e., 40 MHz) has provided impressive results of the possible presence of voids and fractured rocks down to the depth of ten meters. The absence of vegetation, obstacles, and structures in the area of investigation has allowed the acquisition of clean data, with conditions of acquisition that proved to be ideal for the use of this no-shielded antenna. The possibility to work with this type of antenna, not coupled to the ground, provides the ability to investigate large areas with high accuracy, representing an effective tool for the detection and characterization of fractured rocks and caves in the Apulian carbonate terrains. Three-dimensional representations have supported the identification of the main anomalies present in the shallower layers of subsoil, allowing fundamental information about the possible locations of cave entrances, or at least about the area where it is more recommendable to locate drillings with the aim to identify the highly fractured rocks. The use of the Hilbert transform, further, highlights with more clarity the fractured areas, and its implementation can simplify the comparison between GPR and other low resolution geophysical methodologies (i.e., ERT, seismic, …).

With regards to the ERT investigation strategy, we can confirm:

• -. the better performance of the DD array compared to those of other geoelectrical arrays, in agreement with [91];

• -. the importance of adopting a multi-scale approach, resulting in an increase in the target resolution also at depth [93],

• -. the necessity of a 3D or quasi 3D geoelectrical modeling for a clear understanding of the real spatial development of a karst system, that is rarely 2D in nature.

We can also affirm that, even if each of the three levels has its own consistency and can be applied individually, we retain that the Gravaglione results prove how, in the whole, their combined use can greatly improve the quality of the retrievable information.

In detail, the multi-scale/multi-resolution approach, and the quasi-3D geoelectrical modeling, respectively, and the second and third level of the investigation strategy, take strong advantage from the outcomes of the first level. The identification of the most effective geoelectrical array allows, in fact, for minimizing both the acquisition phase time, avoiding the acquisition of different arrays of datasets, and data analysis and modeling time.

Author Contributions

Conceptualization: L.C., I.S.L., M.P. and G.R.; Methodology: N.A., L.C., I.S.L., M.P. and G.R.; Formal analysis, N.A., L.C., M.D.G., D.P. and G.R.; Investigation: N.A., L.C., M.D.G., D.P. and G.R.; Writing—original draft preparation, N.A., L.C., I.S.L. and G.R.; Writing—review and editing: L.C.; I.S.L., M.P. and G.R.; Visualization, N.A., L.C., M.D.G., D.P. and G.R.; Funding acquisition, M.P. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Geophysical and UAS data supporting the results can be found by contacting the corresponding author.

Acknowledgments

The authors greatly thank Giuseppe Consoli, owner of the fields surrounding the Gravaglione site, for his great availability during the surveys.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Geological map of Apulia Region with details of karst zones and study site location (mod. after [28]).

Figure 2. Contour diagrams for the structural characterization of primary and secondary discontinuities in the carbonate rock masses at the surface, highlighting the preferential direction of discontinuities in SW–NE direction, according to the anti-Apennine tectonic regional trend [39].

Figure 3. Karst geomorphological map of Canale di Pirro polje, realized in collaboration with CNR-IRPI researchers, showing the superficial morphological and karst features (mod. after [31]). Gravaglione and Rotolo swallow holes are highlighted in red.

Figure 4. The Gravaglione swallow hole after a rainstorm in December 2020. Red circles highlight the main adsorption areas and the water entrances to the underground (photo courtesy by Rodrigo Torres).

Figure 5. Geophysical survey in the area surrounding the Gravaglione swallow hole.

Figure 6. Data analysis chain adopted for the processing of the GPR data (a) raw radargram acquired at the Gravaglione site (b) processed radargram with topographic correction and migration (c,d) and with application of the Hilbert Transform (e).

Figure 7. The photogrammetric survey of the area: (a,b) colored point cloud and colored point cloud for relative height (bird’s eye view); (c) RGB orthophoto (1.8 cm/pixel); (d) DSM (Digital Surface Model) (14 cm/pixel).

Figure 8. Location of the radargrams F2 and F6 (a); pradargram F2 with indication (highlighted by white dashed rectangles) of the anomalies C1 and C2 (b); radargram F6 with indication (highlighted by white dashed rectangles) of the anomalies C1, C3 and C4 (c).

Figure 9. Location of the radargrams F3 and F5 (a) radargram F3 with indication (highlighted by white dashed rectangles) of the anomalies C5 and C6 (b) radargram F6 with indication (highlighted by white dashed rectangles) of the anomalies C2 and C7 (c).

Figure 10. 3D visualization of the 40 MHz radargrams acquired in the area, overlapped to the aerial image of the Gravaglione area (a) with focus of the most reflective anomalies characterizing the west sector (b) and location of the main GPR anomalies in plan (white gray-red circles) (c).

Figure 11. T1 inversion results: WS array, R.M.S. 1.10 (a); DD array, R.M.S. 1.18 (b); Mixed array, R.M.S. 0.82 (c).

Figure 12. Inversion results obtained by inverting T1 dataset (upper panel), M1 dataset (central panel, R.M.S. 1.43) and M2 dataset (lower panel, R.M.S. 1.39). In the M1 section, the continuous black line indicates the area also investigated by T1. In the M2 section, the continuous and dotted black lines indicate, respectively, the areas also covered by T1 and M1; the dashed purple line indicates the area covered by T9.

Figure 13. Inversion results obtained by inverting T9 dataset (R.M.S. 1.32).

Figure 14. Multiscale survey inversion result. Upper panel: experimental data distribution. Lower panel: Multi-scale geoelectrical model (R.M.S. 1.82).

Figure 15. Depth slices extracted from the 3D model every 10 m of depth. White lines on each slide indicate the traces of the profile T1–T8. The color scale and the resistivity range here adopted differs for practical reasons by those adopted in Figure 11, Figure 12, Figure 13 and Figure 14. In the uppermost part of the figure, the DSM of the modeled area, as derived by the UAS survey, is plotted.

Figure 16. Geoelectrical section and isovolumes from the 3D geoelectrical model. As indicated by the label, the sections were extracted in correspondence of the 2D tomographies T1, T2, T3 and T4. More specifically, shows the isovolume that includes model areas associated to resistivity values higher than 1000 Ω*m (a); shows the isovolume that includes model areas associated to resistivity values higher than 500 Ω*m (b); summarizes the previous panels by showing at the same time both the isovolumes (c).

Figure 17. Localization of ERT8 and F2 radargram in the Gravaglione area (a). F2 radargram individuates interesting electromagnetic anomalies (b) clearly highlighted processing the GPR data by the Hilbert Transform (c). The GPR anomalies are in agreement with the superficial ones detected by the ERT8 (d).

Figure 18. UAS derived DSM presented by using a limited elevation range (308.5 m, 309.0 m) with indication of the locations of GPR anomalies (white ovals) (a); superimposition of the >1000 Ω*m isovolume, GPR anomalies locations and DSM isolines related different elevations (b); superimposition of the >500 Ω*m isovolume, GPR anomalies locations and DSM isolines related different elevations (c). In (b,c) the isovolumes coloring scheme is reported next to panel (b). The coloring scheme reflects the elevation of the isovolumes.

References

1. Ford, D.; Williams, P. Karst Hydrogeology and Geomorphology; Wiley: Chichester, UK, 2007; [DOI: https://dx.doi.org/10.1002/9781118684986]

2. Parise, M.; Ravbar, N.; Živanovic, V.; Mikszewski, A.; Kresic, N.; Mádl-Szőnyi, J.; Kukuric, N. Hazards in Karst and Managing Water Resources Quality. Karst Aquifers—Characterization and Engineering; Stevanovic, Z. Professional Practice in Earth Sciences Springer: Berlin/Heidelberg, Germany, 2015; Chapter 17 pp. 601-687. [DOI: https://dx.doi.org/10.1007/978-3-319-12850-4_17]

3. Parise, M.; Gabrovsek, F.; Kaufmann, G.; Ravbar, N. Recent advances in karst research: From theory to fieldwork and applications. Advances in Karst Research: Theory, Fieldwork and Applications; Parise, M.; Gabrovsek, F.; Kaufmann, G.; Ravbar, N. Geological Society, London, Special Publications: London, UK, 2018; Volume 466, pp. 1-24. [DOI: https://dx.doi.org/10.1144/SP466.26]

4. Williams, P.W. The role of the epikarst in karst and cave hydrogeology: A review. Int. J. Speleol.; 2008; 37, pp. 1-10. [DOI: https://dx.doi.org/10.5038/1827-806X.37.1.1]

5. Bonacci, O. Karst Hydrogeology with Special Reference to the Dinaric Karst; Springer: New York, NY, USA, 1987.

6. Dreybrodt, W. Processes in Karst Systems—Physics, Chemistry and Geology; Springer: New York, NY, USA, 1988.

7. White, W. Geomorphology and Hydrology of Karst Terrains; Geological Magazine Oxford University Press: New York, NY, USA, 1988; Volume 126, pp. 203-204. [DOI: https://dx.doi.org/10.1017/S0016756800006415]

8. Klimchouk, A.; Ford, D.C.; Palmer, A.N.; Dreybrodt, W. Evolution of Karst Aquifers. The Formation of Epikarst and Its Role in Vadose Speleogenesis; National Speleological Society: Huntsville, AL, USA, 2000; pp. 91-99. [DOI: https://dx.doi.org/10.1016/S0022-1694(00)00341-3]

9. Abusaada, M.; Sauter, M. Studying the flow dynamics of a karst aquifer system with an equivalent porous medium model. Ground Water; 2013; 51, pp. 641-650. [DOI: https://dx.doi.org/10.1111/j.1745-6584.2012.01003.x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23039080]

10. Anaya, A.; Padilla, I.; Macchiavelli, R.; Vesper, D.; Meeker, J.; Alshawabkeh, A. Estimating Preferential Flow in Karstic Aquifers Using Statistical Mixed Models. Ground Water; 2014; 52, pp. 584-596. [DOI: https://dx.doi.org/10.1111/gwat.12084] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23802921]

11. Tamburini, A.; Menichetti, M. Groundwater Circulation in Fractured and Karstic Aquifers of the Umbria-Marche Apennine. Water; 2020; 12, 1039. [DOI: https://dx.doi.org/10.3390/w12041039]

12. Goldscheider, N.; Drew, D.; Worthington, S. Methods in Karst Hydrogeology. International Contributions to Hydrogeology; CRC Press: London, UK, 2007; Volume 26, 264. [DOI: https://dx.doi.org/10.1201/9781482266023]

13. Gutierrez, F.; Parise, M.; De Waele, J.; Jourde, H. A review on natural and human-induced geohazards and impacts in karst. Earth Sci. Rev.; 2014; 138, pp. 61-88. [DOI: https://dx.doi.org/10.1016/j.earscirev.2014.08.002]

14. Parise, M. Sinkholes. Encyclopedia of Caves; 3rd ed. White, W.B.; Culver, D.C.; Pipan, T. Academic Press: Cambridge, MA, USA, Elsevier: Amsterdam, The Netherlands, 2019; pp. 934-942. ISBN 978-0-12-814124-3

15. Parise, M. Sinkholes, Subsidence and Related Mass Movements. Treatise on Geomorphology; Shroder, J.J.F. Elsevier: Amsterdam, The Netherlands, Academic Press: Cambridge, MA, USA, 2022; Volume 5, pp. 200-220. [DOI: https://dx.doi.org/10.1016/B978-0-12-818234-5.00029-8] ISBN 9780128182345

16. Mirzanejad, M.; Tran Khiem, T.; McVay, M.; Horhota, D.; Wasman Scott, J. Sinkhole detection with 3D full seismic waveform tomography. Geophysics; 2020; 85, pp. 169-179. [DOI: https://dx.doi.org/10.1190/geo2019-0490.1]

17. Martínez-Moreno, F.J.; Galindo Zaldívar, J.; Pedrera, A.; Teixidó, T.; Ruano, P.; Peña, J.A.; González-Castillo, L.; Constán-Ruiz, A.; López-Chicano, M.; Martín-Rosales, W. Integrated geophysical methods for studying the karst system of Gruta de las Maravillas (Aracena, Southwest Spain). J. Appl. Geophys.; 2014; 107, pp. 149-162. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2014.05.021]

18. Cheng, Q.; Chen, X.; Tao, M.; Binley, A. Characterization of karst structures using quasi-3D electrical resistivity tomography. Environ. Earth Sci.; 2019; 78, 285. [DOI: https://dx.doi.org/10.1007/s12665-019-8284-2]

19. Pazzi, V.; Di Filippo, M.; Di Nezza, M.; Carlà, T.; Bardi, F.; Marini, F.; Fontanelli, K.; Intrieri, E.; Fanti, R. Integrated geophysical survey in a sinkhole-prone area: Microgravity, electrical resistivity tomographies, and seismic noise measurements to delimit its extension. Eng. Geol.; 2018; 234, pp. 282-293. [DOI: https://dx.doi.org/10.1016/j.enggeo.2018.07.016]

20. Stevanović, Z. Global Distribution and Use of Water from Karst Acquifers; Geological Society, London, Special Publications: London, UK, 2018; Volume 466, pp. 217-236. [DOI: https://dx.doi.org/10.1144/SP466.17]

21. Olarinoye, T.; Gleeson, T.; Marx, V.; Seeger, S.; Addingham, R.; Allocca, V.; Andreo, B.; Apaéstegui, J.; Apolit, C.; Arfib, B. et al. Global karst springs hydrograph dataset for research and management of the world’s fastest flowing groundwater. Sci. Data; 2020; 7, 59. [DOI: https://dx.doi.org/10.1038/s41597-019-0346-5]

22. Tropeano, M.; Caldara, M.A.; De Santis, V.; Festa, V.; Parise, M.; Sabato, L.; Spalluto, L.; Francescangeli, R.; Iurilli, V.; Mastronuzzi, G.A. et al. Geological uniqueness and potential geotouristic Appeal of Murge and Premurge, the first territory in Puglia (Southern Italy) aspiring to become a UNESCO Global Geopark. Geosciences; 2023; 13, 131. [DOI: https://dx.doi.org/10.3390/geosciences13050131]

23. Ciaranfi, N.; Pieri, P.; Ricchetti, G. Nota alla Carta Geologica delle Murge e del Salento (Puglia meridionale). Mem. Della Soc. Geol. Ital.; 1988; 41, pp. 449-460.

24. Funiciello, R.; Montone, P.; Parotto, M.; Salvini, F.; Tozzi, M. Geodynamical evolution of an intra-orogenic foreland: The Apulia case history (Italy). Boll. Della Soc. Geol. Ital.; 1991; 110, pp. 419-425.

25. Tropeano, M.; Sabato, L.; Pieri, P. The Quaternary ‘Post-turbidite’ sedimentation in the South-Apennines Foredeep (Bradanic Trough-Southern Italy). Boll. Della Soc. Geol. Ital.; 2002; 121, pp. 449-454.

26. Pepe, M.; Parise, M. Structural control on development of karst landscape in the Salento Peninsula (Apulia, SE Italy). Acta Carsologica; 2014; 43, pp. 101-114. [DOI: https://dx.doi.org/10.3986/ac.v43i1.643]

27. Parise, M. Surface and subsurface karst geomorphology in the Murge (Apulia, southern Italy). Acta Carsologica; 2011; 40, pp. 79-93. [DOI: https://dx.doi.org/10.3986/ac.v40i1.30]

28. Cotecchia, V. Le acque sotterranee e l’intrusione marina in Puglia: Dalla ricerca all’emergenza nella salvaguardia della risorsa. Mem. Descr. Carta Geol. D’italia; 2014; 92, 416.

29. Parise, M.; Inguscio, S.; Marangella, A. Elementi di geomorfologia carsica della Puglia. Geomorfologia Carsica, Proceedings of the Grottaglie, Atti 45° corso CNSS-SSI, Grottaglie, Italy, 2–3 February 2008; Litografia Ettorre: Grottaglie, Italy, 2008.

30. Zumpano, V.; Pisano, L.; Parise, M. An integrated framework to identify and analyze karst sinkholes. Geomorphology; 2019; 332, pp. 213-225. [DOI: https://dx.doi.org/10.1016/j.geomorph.2019.02.013]

31. Pisano, L.; Zumpano, V.; Liso, I.; Parise, M. Geomorphological characterization of the ‘Canale di Pirro’ polje, Apulia (Southern Italy). J. Maps; 2020; 16, pp. 479-487. [DOI: https://dx.doi.org/10.1080/17445647.2020.1778550]

32. Micheletti, F.; Fornelli, A.; Spalluto, L.; Parise, M.; Gallicchio, S.; Tursi, F.; Festa, V. Petrographic and geochemical inferences for genesis of Terra Rossa: A case study from the Apulian Karst (southern Italy). Minerals; 2023; 13, 499. [DOI: https://dx.doi.org/10.3390/min13040499]

33. Parise, M. Flood history in the karst environment of Castellana-Grotte (Apulia, southern Italy). Nat. Hazards Earth Syst. Sci.; 2003; 3, pp. 593-604. [DOI: https://dx.doi.org/10.5194/nhess-3-593-2003]

34. Mossa, M. The floods in Bari: What history should have taught. J. Hydraul. Res.; 2007; 45, pp. 579-594. [DOI: https://dx.doi.org/10.1080/00221686.2007.9521795]

35. Martinotti, M.; Pisano, L.; Marchesini, I.; Rossi, M.; Peruccacci, S.; Brunetti, M.; Melillo, M.; Amoruso, G.; Loiacono, P.; Vennari, C. et al. Landslides, floods and sinkholes in a karst environment: The 1–6 September 2014 Gargano event, southern Italy. Nat. Hazards Earth Syst. Sci.; 2017; 17, pp. 467-480. [DOI: https://dx.doi.org/10.5194/nhess-17-467-2017]

36. Gentile, P.; Iaia, C.; Liso, I.; Parise, M. Eventi alluvionali nell’ambiente carsico pugliese. SIGEA Geol. Dell’ambiente; 2020; 1, pp. 56-63.

37. Parise, M. Geomorphology of the Canale di Pirro karst polje (Apulia, Southern Italy). Z. Für Geomorphol.; 2006; 147, pp. 143-158.

38. Di Geronimo, I. Geomorfologia del versante adriatico delle Murge di SE (zona di Ostuni, Brindisi). Geol. Romana; 1970; 9, pp. 47-57.

39. Liso, I.S. Hydrogeological Characterization and Groundwater Modeling in the Apulian Karst (Southern Italy). Ph.D. Thesis; University of Bari Aldo Moro: Bari, Italy, 2022.

40. Parise, M.; Benedetto, L. Surface landforms and speleological investigation for a better understanding of karst hydrogeological processes: A history of research in southeastern Italy. Advances in Karst Research: Theory, Fieldwork and Applications; Parise, M.; Gabrovsek, F.; Kaufmann, G.; Ravbar, N. Geological Society, London, Special Publications: London, UK, 2018; Volume 466, pp. 137-153. [DOI: https://dx.doi.org/10.1144/SP466.25]

41. Liso, I.; Chieco, M.; Fiore, A.; Pisano, L.; Parise, M. Underground geosites and caving speleotourism: Some considerations, from a case study in Southern Italy. Geoheritage; 2020; 12, 13. [DOI: https://dx.doi.org/10.1007/s12371-020-00424-z]

42. Parise, M.; Benedetto, L.; Chieco, M.; Fiore, A.; Lacarbonara, M.; Liso, I.S.; Masciopinto, C.; Pisano, L.; Riccio, A.; Vurro, M. First outcomes of a project dedicated to monitoring groundwater resources in Apulia, southern Italy, 2020. Eurokarst 2018, Besançon: Advances in the Hydrogeology of Karst and Carbonate Reservoirs; Bertrand, C.; Denimal, S.; Steinmann, M.; Renard, P. Advances in Karst Science Springer: Berlin/Heidelberg, Germany, 2019; pp. 243-249. ISBN 978-3-030-14014-4

43. Gams, I. Tectonic impacts on poljes and minor basins (case studies of Dinaric Karst). Acta Carsologica; 2005; 34, pp. 25-41. [DOI: https://dx.doi.org/10.3986/ac.v34i1.277]

44. Gams, I. The polje: The problem of definition. Z. Fur Geomorphol.; 1978; 22, pp. 170-181.

45. Lopez Chicano, M.; Calvache, M.; Martin-Rosales, W.; Gisbert, J. Conditioning factors in flooding of karstic poljes—The case of Zafarraya Polje (South Spain). Catena; 2002; 49, pp. 331-352. [DOI: https://dx.doi.org/10.1016/S0341-8162(02)00053-X]

46. Gracia, F.; Gutierrez, F.; Gutierrez, M. The Jiloca karst polje—Tectonic graben (Iberian Range, NE Spain). Geomorphology; 2003; 52, pp. 215-231. [DOI: https://dx.doi.org/10.1016/S0169-555X(02)00257-X]

47. Nicod, J. A little contribution to the karst terminology: Special or aberrant cases of poljes?. Acta Carsologica; 2003; 32, pp. 29-39. [DOI: https://dx.doi.org/10.3986/ac.v32i2.334]

48. Gunn, J. Encyclopedia of Caves and Karst Science; Taylor and Francis: New York, NY, USA, 2004; pp. 599-600.

49. Milanovic, P. Dinaride poljes. Encyclopedia of Caves and Karst Science; Gunn, J. Taylor and Francis: New York, NY, USA, 2004; pp. 291-293.

50. Dogan, U. Sariot Polje, central Taurus (Turkey): A border polje developed at the contact of karstic and non-karstic lithologies. Cave Karst Sci.; 2003; 30, pp. 117-124.

51. Durdağ, D.; Drahor, M.G.; Yağlıdere, M.S. Combined application of electrical resistivity tomography and multi-channel analysis of surface waves methods in the tunnel detection: A case study from Kocaeli University Campus Site, Turkey. J. Appl. Geophys.; 2023; 208, 104895. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2022.104895]

52. Roncella, R.; Forlani, G.; Remondino, F. Photogrammetry for Geological Applications: Automatic Retrieval of Discontinuity Orientation in Rock Slopes. Videometrics VIII; 2005; 5665, 17. [DOI: https://dx.doi.org/10.1117/12.587822]

53. Voyat, I.; Roncella, R.; Forlani, G.; Ferrero, A.M. Advanced Techniques for Geo Structural Surveys in Modelling Fractured Rock Masses: Application to Two Alpine Sites. Proceedings of the GoldenRocks 2006: 41st U.S. Rock Mechanics Symposium; Golden, CO, USA, 17–21 June 2006.

54. Westoby, M.J.; Brasington, J.; Glasser, N.F.; Hambrey, M.J.; Reynolds, J.M. “Structure-from-Motion” Photogrammetry: A Low-Cost, Effective Tool for Geoscience Applications. Geomorphology; 2012; 179, pp. 300-314. [DOI: https://dx.doi.org/10.1016/j.geomorph.2012.08.021]

55. Buyer, A.; Schubert, W. Calculation the Spacing of Discontinuities from 3D Point Clouds. Procedia Eng.; 2017; 191, pp. 270-278. [DOI: https://dx.doi.org/10.1016/j.proeng.2017.05.181]

56. Loiotine, L.; Wolff, C.; Wyser, E.; Andriani, G.F.; Derron, M.-H.; Jaboyedoff, M.; Parise, M. QDC-2D: A Semi-Automatic Tool for 2D Analysis of Discontinuities for Rock Mass Characterization. Remote Sens.; 2021; 13, 5086. [DOI: https://dx.doi.org/10.3390/rs13245086]

57. Colomina, I.; Molina, P. Unmanned aerial systems for photogrammetry and remote sensing: A review. ISPRS J. Photogramm. Remote Sens.; 2014; 92, pp. 79-97. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2014.02.013]

58. Derenyi, E. Photogrammetry: The Concepts; Department of Geodesy and Geomatics Engineering Lecture Notes, University of New Brunswick: Fredericton, NB, Canada, 1996; Volume 220, 57.

59. Remondino, F.; El-Hakim, S. Image-based 3D modelling: A review. Photogramm. Rec.; 2006; 21, pp. 269-291. [DOI: https://dx.doi.org/10.1111/j.1477-9730.2006.00383.x]

60. Mikhail, E.; Bethel, J.; McGlone, J. Introduction to Modern Photogrammetry; John Wiley & Sons: Chichester, UK, New York, NY, USA, 2001.

61. Abate, N.; Ronchi, D.; Vitale, V.; Masini, N.; Angelini, A.; Giuri, F.; Minervino Amodio, A.; Gennaro, A.M.; Ferdani, D. Integrated Close Range Remote Sensing Techniques for Detecting, Documenting, and Interpreting Lost Medieval Settlements under Canopy: The Case of Altanum (RC, Italy). Land; 2023; 12, 310. [DOI: https://dx.doi.org/10.3390/land12020310]

62. Nappi, R.; Paoletti, V.; D’Antonio, D.; Soldovieri, F.; Capozzoli, L.; Ludeno, G.; Porfido, S.; Michetti, A.M. Joint Interpretation of Geophysical Results and Geological Observations for Detecting Buried Active Faults: The Case of the “Il Lago” Plain (Pettoranello del Molise, Italy). Remote Sens.; 2021; 13, 1555. [DOI: https://dx.doi.org/10.3390/rs13081555]

63. Oskooi, B.; Parnow, S.; Smirnov, M.; Varfinezhad, R.; Yari, M. Attenuation of random noise in GPR data by image processing. Arab. J. Geosci.; 2018; 11, 677. [DOI: https://dx.doi.org/10.1007/s12517-018-4035-z]

64. Hager, J.; Carnevale, M. The Application of Low Frequency GPR to Stratigraphic Investigations; Hager Geosci Inc.: Woburn, MA, USA, 2006; pp. 2-7.

65. Bristow, C.; Jewell, C. GPR studies in the Piano di Pezza area of the Ovindoli-Pezza Fault, Central Apennines, Italy: Extending paleoseismic trench investigations with high resolution GPR profiling. Near Surf. Geophys.; 2004; 4, pp. 147-153. [DOI: https://dx.doi.org/10.3997/1873-0604.2005040]

66. Finizola, A.; Ricci, T.; Rossi, M.; Praticelli, N.; Giocoli, A.; Romano, G.; Delcher, E.; Suski, B.; Revil, A.; Menny, P. et al. Adventive hydrothermal circulation on Stromboli volcano (Aeolian Islands, Italy) revealed by geophysical and geochemical approaches: Implications for general fluid flow models on volcanoes. J. Volcanol. Geotherm. Res.; 2010; 196, pp. 111-119. [DOI: https://dx.doi.org/10.1016/j.jvolgeores.2010.07.022]

67. Ercoli, M.; Cirillo, D.; Pauselli, C.; Jol, H.; Brozzetti, F. Ground-penetrating radar signature of Quaternary faulting: A study from the Mt. Pollino region, southern Apennines, Italy. Solid Earth; 2021; 12, pp. 2573-2596. [DOI: https://dx.doi.org/10.5194/se-12-2573-2021]

68. Micallef, A.; Saadatkhah, N.; Spiteri, J.; Rizzo, E.; Capozzoli, L.; De Martino, G.; Pace, L.; Burrill, C.; Gatt, M.; Schmidt, M. et al. Groundwater seepage is a key driver of theater-headed valley formation in limestone. Geology; 2022; 50, pp. 686-690. [DOI: https://dx.doi.org/10.1130/G49938.1]

69. Annan, A.P. GPR Methods for Hydrogeological Studies. Hydrogeophysics; Rubin, Y.; Hubbard, S.S. Water Science and Technology Library Springer: Dordrecht, The Netherlands, 2005; 50. [DOI: https://dx.doi.org/10.1007/1-4020-3102-5_7]

70. Alcalá, F.J.; Martínez-Pagán, P.; Paz, M.C.; Navarro, M.; Pérez-Cuevas, J.; Domingo, F. Combining of MASW and GPR Imaging and Hydrogeological Surveys for the Groundwater Resource Evaluation in a Coastal Urban Area in Southern Spain. Appl. Sci.; 2021; 11, 3154. [DOI: https://dx.doi.org/10.3390/app11073154]

71. Paz, C.; Alcalá, F.J.; Carvalho, J.; Ribeiro, L. Current uses of ground penetrating radar in groundwater-dependent ecosystems research. Sci. Total Environ.; 2017; 595, pp. 868-885. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2017.03.210]

72. Anchuela, Ó.; Casas-Sainz, A.; Soriano, M.; Pocoví-Juan, A. Mapping subsurface karst features with GPR: Results and limitations. Environ. Geol.; 2009; 58, pp. 391-399. [DOI: https://dx.doi.org/10.1007/s00254-008-1603-7]

73. Čeru, T.; Šegina, E.; Gosar, A. Geomorphological Dating of Pleistocene Conglomerates in Central Slovenia Based on Spatial Analyses of Dolines Using LiDAR and Ground Penetrating Radar. Remote Sens.; 2017; 9, 1213. [DOI: https://dx.doi.org/10.3390/rs9121213]

74. Margiotta, S.; Negri, S.; Parise, M.; Quarta, T.A.M. Karst geosites at risk of collapse: The sinkholes at Nociglia (Apulia, SE Italy). Environ. Earth Sci.; 2016; 75, 8. [DOI: https://dx.doi.org/10.1007/s12665-015-4848-y]

75. Margiotta, S.; Marini, G.; Fay, S.; D’Onghia, F.M.; Liso, I.S.; Parise, M.; Pinna, M. Hydro-stratigraphic conditions and human activity leading to development of a sinkhole cluster in a Mediterranean water ecosystem. Hydrology; 2021; 8, 111. [DOI: https://dx.doi.org/10.3390/hydrology8030111]

76. Caselle, C.; Bonetto, S.; Comina, C.; Stocco, S. GPR surveys for the prevention of karst risk in underground gypsum quarries. Tunn. Undergr. Space Technol.; 2020; 95, pp. 103-137. [DOI: https://dx.doi.org/10.1016/j.tust.2019.103137]

77. Bermejo, L.; Ortega, A.; Parés, J.; Campaña, I.; Castro, J.; Carbonell, E.; Conyers, L. Karst features interpretation using ground-penetrating radar: A case study from the Sierra de Atapuerca, Spain. Geomorphology; 2020; 367, 107311. [DOI: https://dx.doi.org/10.1016/j.geomorph.2020.107311]

78. Bélanger, C.; Bernard, G.; Erwan, G.; Lefebvre, R. GPR, ERT and CPT data integration for high resolution aquifer modeling. Proceedings of the 13th Internarional Conference on Ground Penetrating Radar; Lecce, Italy, 21–25 June 2010; [DOI: https://dx.doi.org/10.1109/ICGPR.2010.5550182]

79. Al-Fares, W.; Bakalowicz, M.; Guérin, R.; Dukhan, M. Analysis of the karst aquifer structure of the Lamalou area (Hérault, France) with ground penetrating radar. J. Appl. Geophys.; 2002; 51, pp. 97-106. [DOI: https://dx.doi.org/10.1016/S0926-9851(02)00215-X]

80. Baradello, L.; Accanio, F.; Pepi, S.; Vaccaro, C.; Nieto, D. Geophysical Methods for the Assessment of the Ejection of Sand Due to the Liquefaction in the Mirabello Area. Proceedings of the Near Surface Geoscience—20th European Meeting of Environmental and Engineering Geophysics; Athens, Greece, 14–18 September 2014; [DOI: https://dx.doi.org/10.3997/2214-4609.20141961]

81. Mohamed, A.; El-Hussain, I.; Deif, A.; Araffa, S.; Mansour, K.; Al-Rawas, G. Integrated ground penetrating radar, electrical resistivity tomography and multichannel analysis of surface waves for detecting near-surface caverns at Duqm area, Sultanate of Oman. Near Surf. Geophys.; 2019; 17, pp. 379-401. [DOI: https://dx.doi.org/10.1002/nsg.12054]

82. Zhu, J.; Currens, J.; Dinger, J. Challenges of using electrical resistivity method to locate karst conduits: A field case in the Inner Bluegrass Region, Kentucky. J. Appl. Geophys.; 2011; 75, pp. 523-530. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2011.08.009]

83. Iovine, G.; Vennari, C.; Gariano, S.L.; Caloiero, T.; Lanza, G.; Nicolino, N.; Suriano, S.; Ferraro, G.; Parise, M. The “Piano dell’Acqua” sinkholes (San Basile, Northern Calabria, Italy). Bull. Eng. Geol. Environ.; 2016; 75, pp. 37-52. [DOI: https://dx.doi.org/10.1007/s10064-015-0737-6]

84. Robert, T.; Dassargues, A.; Brouyère, S.; Kaufmann, O.; Hallet, V.; Nguyen, F. Assessing the contribution of electrical resistivity tomography (ERT) and self-potential (SP) methods for a water well drilling program in fractured/karstified limestones. J. Appl. Geophys.; 2011; 75, pp. 42-53. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2011.06.008]

85. Bermejo, L.; Ana Isabel, O.; Guérin, R.; Benito-Calvo, A.; Pérez-Gonzàlez, A.; Parés, J.; Aracil, E.; Bermùdez de Castro, J.M.; Carbonell, E. 2D and 3D ERT imaging for identifying karst morphologies in the archaeological sites of Gran Dolina and Galería Complex (Sierra de Atapuerca, Burgos, Spain). Quat. Int.; 2016; 433, pp. 393-401. [DOI: https://dx.doi.org/10.1016/j.quaint.2015.12.031]

86. McCormack, T.; O’Connel, Y.; Daly, E.; GIll, L.; Henry, T.; Perriquet, M. Characterisation of karst hydrogeology in Western Ireland using geophysical and hydraulic modelling techniques. J. Hydrol. Reg. Stud.; 2017; 10, pp. 1-17. [DOI: https://dx.doi.org/10.1016/j.ejrh.2016.12.083]

87. Keshavarzi, M.; Baker, A.; Kelly, B.; Andersen, M. River–groundwater connectivity in a karst system, Wellington, New South Wales, Australia. Hydrogeol. J.; 2017; 25, pp. 557-574. [DOI: https://dx.doi.org/10.1007/s10040-016-1491-y]

88. Tarquini, S.; Nannipieri, L. The 10m-resolution TINITALY DEM as a trans-disciplinary basis for the analysis of the Italian territory: Current trends and new perspectives. Geomorphology; 2016; 281, pp. 108-115. [DOI: https://dx.doi.org/10.1016/j.geomorph.2016.12.022]

89. Tarquini, S.; Isola, I.; Favalli, M.; Battistini, A.; Dotta, G. TINITALY, a Digital Elevation Model of Italy with a 10 Meters Cell Size (Version 1.1) [Data Set]; Istituto Nazionale di Geofisica e Vulcanologia (INGV): Naples, Italy, 2023; [DOI: https://dx.doi.org/10.13127/tinitaly/1.1]

90. Zheng, L.; Liu, Z.; Wang, G.; Zhang, Z. Research on application of Hilbert transform in radar signal simulation. Proceedings of the 7th International Conference on Environmental and Engineering Geophysics & Summit Forum of Chinese Academy of Engineering on Engineering Science and Technology, Beijing, China, 26–29 June 2016; Beijing, China, 26–29 June 2016; Atlantis Press: Paris, France, 2016; pp. 347-349. [DOI: https://dx.doi.org/10.2991/iceeg-16.2016.93]

91. Tso, C.-H.M.; Kuras, O.; Wilkinson, P.B.; Uhlemann, S.; Chambers, J.E.; Meldrum, P.I.; Binley, A. Improved characterisation and modelling of measurement errors in electrical resistivity tomography (ERT) surveys. J. Appl. Geophys.; 2017; 146, pp. 103-119. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2017.09.009]

92. Blanchy, G.; Saneiyan, S.; Boyd, J.; McLachlan, P.; Binley, A. ResIPy, an Intuitive Open Source Software for Complex Geoelectrical Inversion/Modeling. Comput. Geosci.; 2020; 137, 104423. [DOI: https://dx.doi.org/10.1016/j.cageo.2020.104423]

93. Doyoro, Y.G.; Chang, P.-Y.; Puntu, J.M. Uncertainty of the 2D Resistivity Survey on the Subsurface Cavities. Appl. Sci.; 2021; 11, 3143. [DOI: https://dx.doi.org/10.3390/app11073143]

94. Zhou, W.; Beck, B.; Adams, A. Effective electrode array in mapping karst hazards in electrical resistivity tomography. Environ. Geol.; 2002; 42, pp. 922-928. [DOI: https://dx.doi.org/10.1007/s00254-002-0594-z]

95. Diaferia, I.; Barchi, M.; Loddo, M.; Schiavone, D.; Siniscalchi, A. Detailed imaging of tectonic structures by multiscale Earth resistivity tomographies: The Colfiorito normal faults (central Italy). Geophys. Res. Lett.; 2006; 33, pp. 1-4. [DOI: https://dx.doi.org/10.1029/2006GL025828]

96. Chávez, R.E.; Tejero, A.; Cifuentes, G.; Hernández, E.; Aguilar, D. Imaging fractures beneath a residential complex using novel 3-D electrical resistivity arrays. J. Environ. Eng. Geophys.; 2015; 20, pp. 219-233. [DOI: https://dx.doi.org/10.2113/JEEG20.3.219]

97. Kaufmann, O.; Deucester, J. Detection and mapping of ghost rock features in the Tournaisis area through geophysical methods: An overview. Geol. Belg.; 2014; 17, pp. 17-26.

98. Binley, A. R3t Version 1.8 Manual; Lancaster University: Lancaster, UK, 2013; Available online: http://www.es.lancs.ac.uk/people/amb/Freeware/R3t/R3t.htm (accessed on 1 January 2018).

99. Liso, I.S.; Cherubini, C.; Parise, M. Hydrogeological characterization and modeling at two test sites of the Apulian karst (southern Italy), 2023. Eurokarst 2022, Malaga: Advances in the Hydrogeology of Karst and Carbonate Reservoirs; Andreo, B.; Barberà, J.A.; Duràn-Valsero, J.J.; Gil-Marquez, J.M.; Mudarra, M. Advances in Karst Science Springer: Berlin/Heidelberg, Germany, 2023; ISBN 978-3-031-16878-9