Natural Sinkhole Monitoring and Characterization: The Case of Latera Sinkhole (Latium, Central Italy)

Abstract

The occurrence of sinkhole phenomena in Italy is a prevalent and very uncertain class of geological hazards that pose a significant threat to human infrastructure and individuals. These events are characterized by their unpredictability and the challenges associated with their accurate forecasting. Both natural and anthropic factors influence the occurrence of these events; therefore, accurate identification of the above factors is critical for effective proactive and predictive efforts. The work presented in this paper refers to a collapse that occurred in a volcanic region in northern Latium (central Italy) on 31 January 2023. The area has been monitored using drones since the early stages of the sinkhole’s formation and has continued to date. Then, the collapse and the neighboring area were examined via geophysical and geochemical investigations to identify potential underlying factors. Geophysical and geochemical data were combined to provide a preliminary hypothesis on the collapse’s genesis. The obtained data indicate that the structural collapse can be attributable to the fluctuation in groundwater levels as well as the development of instabilities along its banks, leading to a growth in its dimensions.

1. Introduction

In the technical literature, a wide range of terminology is used to describe different types of collapse. Over the past several decades, the term “sinkhole” has been widely used to encompass a subcircular depression related to an abrupt collapse of the ground surface into an underlying void or cavity, irrespective of its morphogenetic origin or geometric characteristics [1,2,3,4,5]. Sinkholes are common in Mediterranean countries and their formation and geomorphological evolution can be attributed to both natural and anthropogenic factors [6,7,8,9,10,11,12,13,14,15,16,17,18]. They represent one of the most relevant geohazards due to their high level of danger and unpredictable nature and may cause severe damage to infrastructure as well as casualties, especially in urban areas [19]. Numerous attempts were made to classify sinkholes, considering the processes through which they are formed. These efforts include assessments of the terminology employed to designate various forms of geological cavities. The most recent and thorough study on this topic was conducted by Whaltham [9,10]. However, there are many other documented examples of regional inventories, such as those presented by some researchers, e.g., [11,12,13,14,15,16,17,18]. Natural sinkholes mainly occur in regions characterized by the presence of limestone or evaporitic rock formations. The predisposing factors for their formation include the existence of soluble rocks that undergo karst processes, as well as the geotechnical characteristics of the underlying bedrock and the chemistry and flow rate of groundwater in the overlying loose layers, if present. In addition, the upwellings along faults and fractured systems of acidic deep high-pressure fluids containing significant amounts of CO₂ and H₂S can also have the potential to disrupt the carbonate matrix of the shallow rocks and sediments, thus creating the instability of the shallow sedimentary cover [12,13,14]. Furthermore, extreme weather events linked to climate change can also cause significant variations in groundwater levels and flow rates, which can affect soil dissolution and erosion in both natural and anthropic settings [20,21,22,23]. In contrast, anthropogenic sinkholes are generally caused by human activities such as mining, drilling, and construction. In urban areas, underground cavities are considered a common predisposing factor for anthropogenic sinkhole formation [15]. Regardless of its origin, if detected at an early stage, sinkhole formation and evolution can be successfully observed and monitored using indirect methods [24,25,26,27]. The prediction of sinkhole occurrences is challenging due to the typically minimal surface deformations that serve as precursors [27,28,29]. However, some preventive measures should be adopted to assess and possibly mitigate the sinkhole risk. These include defining the geological subsoil setting at a proper scale, monitoring groundwater levels, surveying the area for any anomalies or fractures in the subsoil, and checking the condition of sewer and water systems (in urban areas).

This paper examines the emergence of a natural sinkhole formation near Provincial Road n.118 in the municipality of Latera (northern Latium, central Italy) in late January 2023. The collapse took place within a volcanic region where no previous occurrences of similar phenomena have been documented in the geological record. Although sparsely inhabited, the study area has experienced significant alterations to its hydrographic network and restoration initiatives, which involved extensive excavation and backfilling. The potential consequence of these modifications is the loss of historical evidence, which complicates the risk assessment approach. In the area, technical and scientific activities were conducted to support the local administration in the safe management of the roadway while assessing the residual risks in the surroundings of the collapsed area.

In particular, the following activities were carried out:

• Monthly based monitoring by using an Uncrewed Aerial Vehicle (UAV) to evaluate soil subsidence and study the rise in the groundwater level in the collapse during the first 30 days after its formation.
• Geophysical surveys to define a reliable subsurface model useful to address direct investigations (boreholes, geotechnical tests) and to image the collapsed area to possibly infer considerations on buried structures potentially responsible for the event.
• Geochemical survey to highlight the presence of anomalous concentrations of soil gases, indicating the presence of advective gas migration from depth along major permeability pathways (e.g., faults and fractures) that may trigger the collapse.

The data given refer to continuous eight-month monitoring activities, leading to the construction of a subsurface model supporting a preliminary hypothesis about sinkhole genesis via suffusion [17].

2. Study Area

2.1. General Geological Setting

The study area is located in the Latera Volcanic Complex (LVC), whose activity dates back to about 0.28–0.14 Ma, consisting of a gently sloping, central volcanic edifice about 30 km wide, truncated by a polygenetic caldera associated with several stages of collapse [30,31]. The LVC is part of the Quaternary Vulsini Volcanic District (VVD), a belt of young volcanoes aligned along the western coast of Italy. The northwestern sector of the Latera caldera hosts the minor collapse structure of the Vepe caldera, interpreted as a nested caldera system developed through multiple collapse phases, each following intervals of moderate-magnitude explosive eruptions (Figure 1). Based on the distribution of post-collapse lava domes, this polycyclic caldera had been interpreted “to involve coherent subsidence of a piston-like block” bounded by a ring fracture system. These occurred in four major pulses that produced moderate-volume ignimbrites, pumice-fallout, pyroclastic-surge deposits, lavas, and largely post-caldera cinder and spatter deposits [32,33,34]. According to many authors (see, e.g., [32]), the products referring to the older Vepe’s activity can be distinguished from those derived from the recent activity, respectively, represented by the “Poggio Pinzo” unit and by the “Vulcanite complessa di Pitigliano” (VCP) unit as in Figure 1B.

The ancient Vepe’s activity is represented by several kinds of pyroclastic deposits, mostly consisting of pumice deposits and scoriaceous beds, enclosing ballistic blocks and bombs due to strombolian fallout, interbedded with massive plane- and cross-laminated and dune-bedded deposits of ash and poorly vesicular scorias and lapilli from hydromagmatic surges [33,34]. According to many authors, during the initial activity, a first calderic collapse takes place, marking the beginning of the formation of Vepe’s caldera. The last phase of Vepe’s activity results from the emplacement of the ca. 166–155 ka VCP unit, made by different volcanic sub-units with strong vertical and horizontal variability, fed by a chemically zoned magma chamber poor in volatile components [35,36]. In the final stage of Vepe’s caldera evolution, dome-like structures occurred on the pyroclastic deposits at the base of the VCP. The final phase of the Vepe’s activity is documented by the Mezzano maar structure with associated hydromagmatic products that, in turn, were covered by the Monte Spinaio units representing the last dome-like structures put in place in the area, presumably from more than one emission center (Figure 1A). The alluvial and lacustrine deposits are widely outcropping within the caldera with a thickness of several tens of meters that includes peaty levels with poor geotechnical characteristics [36].

2.2. Stratigraphic Setting

The products of Vepe caldera have been intercepted by the Olpetella well, located some hundreds of meters from the sinkhole. The well has been realized for geothermal exploration with a core-destruction drilling technique; therefore, log descriptions mostly referring to the lithology, grain size, composition, and texture of championed geological units. For the first 16 m, loose covers with an abundant sandy–silty matrix, referred to as alluvial and lacustrine environments, have been recorded in the well-log. This mantle rests on porous but compact lava, jointed at the base close to 36 m in depth, where a 4 m thick pyroclastic deposit has been found. From 40 to 45 m in depth, the well encountered compact lava, resting on an undistinguished coarse and altered pyroclastic layer. From 78 m to 102 m, light-grey brittle pozzolanas were drilled, lying over about 40 m of compact effusive unit (grey lava). The well-log ends at a depth of 196 m, and the basal layer recorded consists of undistinguished and heterogeneous pyroclastic materials, coarse, altered and with increasingly fine grain size to depth (Figure 2B). With only synthetic lithostratigraphic descriptions available, no clear hypotheses can be made on the attribution of the drilled volcanic units to geologic formations described in the literature. No indication about the groundwater level is reported in the well log; however, filters were placed at three intervals starting from 55 m in depth, suggesting the possible presence of a multi-level aquifer. In this regard, note that at a regional scale, the main groundwater table is expected to be located in this area between 440 and 420 m above sea level, i.e., from 10 to 30 m in depth [37].

2.3. Morphological Setting

Due to the persistence of a marshy environment in many areas within the caldera, many reclamation activities took place, in particular during the XVIII century, such as reclamation works focused on the Mezzano Lake and its hydraulic network [36]. Figure 1A shows that the ditches northeast of Mezzano Lake appear clearly rectified. These works also concerned the Olpeta River, which runs just north of the study area, to lower the riverbed by 2 m. Further reclamation works [38] were realized for the last time between 1934 and 1936 in the Olpeta area: the final morphological setting of the study area at present is depicted in Figure 2A, where the hydraulic network had been marked with blue color in order to underline the effects of the reclamation activities. After the analysis of available aerial photos in the same picture (dotted blue), the course of a possible old ditch that ran through the area (before it was channelized) has been recognized.

2.4. Climate

The territory of Latera is characterized by a warm–rainy temperate climate (Csa, according to [39]) with average temperatures between 18° and −3 °C. The dry season lasts from June to August, and the temperatures in the hottest month reach above 22° (these characteristics can also be deduced from the climatograph obtained from the data of the reference rain gauge station for the municipal area, namely the “Latera-Centro Florovivaistico” station of the Regional Agency for Agricultural Development and Innovation of Lazio—ARSIAL (open data available at http://dati.lazio.it/catalog/it/dataset?organization=arsial (accessed on 10 November 2023), located just 2.5 km southeast of the study area). The climatograph shown in Figure 3 was reconstructed taking into consideration the meteorological series reported on the ARSIAL website and referred to the period 2004–2023. The 2005 data were not considered, due to the lack of measurements during the last four months of the year. The average temperature of this area is 13.7 °C, while the annual rainfall slightly exceeds 1000 mm. The rainiest months are November and December, while constant rainfall levels are recorded in the months of January, February and March. As for 2023, the month of January was characterized by rainfall of 154.9 mm, with many events in the first twenty days.

3. Materials and Methods

The characterization of the subsoil in the study area was accomplished by implementing the Electrical Resistivity Tomography (ERT) method first, combined with surface-wave seismic surveys and single station microtremors recordings for Horizontal-to-Vertical Spectral Ratio technique ([40] hereafter referred to as HVSR). After the geophysical campaigns, a preliminary subsoil model was defined, and additional surveys were conducted, consisting of geochemical parameters measurements (natural radioactivity), in order to constrain the subsoil model and to investigate the predisposing factors, e.g., fractures/faults, that could have contributed to the formation of the sinkhole. In the following sections, the geological setting of the study area is briefly described, along with the expected local stratigraphic setting as derived from the well-log for geothermal exploration near the site. The results of photogrammetric monitoring of the sinkhole area are then illustrated, including the estimation of the groundwater balance within the cavity over eight months. The indirect methods applied are briefly described, and the outcomes obtained from geophysical and geochemical surveys are integrated and eventually discussed.

3.1. UAV Photogrammetric Surveys

UAVs surveys were carried out monthly, from February to September 2023, aiming to monitor any variations in the shapes and dimensions of the sinkhole. Two types of Uncrewed Aerial Vehicle (UAV) were used: the Italdron EVO 4hse with a 9 kg take-off configuration and armed with a Sony Alfa 7r mark II with 35 mm optics; and the DJI Mavic Pro, with a 900 g take-off configuration and on-board camera. To obtain topographic data from UAV surveys (see below), all flights were supported by placing Ground Control Points on the ground (three of which fixed and non-removable). The results obtained were also used as a topographical basis for geophysical investigations. The flights were carried out, imposing a photographic overlap of 75–80% between frames at an altitude varying between 30 m and 65 m (depending on the UAV used) to obtain a minimum nominal definition on the ground between 0.9 and 1.3 cm/pixels. The nadiral flight was integrated with passages of camera image acquisitions at 45° and lower altitudes to reconstruct the structure of the sinkhole and facilitate the cleaning and classification of the dense cloud. The images and coordinates thus obtained were processed using the commercial software Agisoft Metashape Professional v.1.6.3 (https://www.agisoft.com, accessed on 10 November 2023), following the steps described below:

• Image control and quality assessment to identify high-quality images to be used;
• Image alignment and preliminary sparse cloud elaboration;
• Attribution of a reference system (EPSG 3003) and positioning of marks in photographic images, with subdivision into ground control point (GCP) and control point (C.P.);
• Process of alignment optimization and gradual selection for error reduction;
• Heights map from high-quality, dense cloud generation;
• Elimination of low confidence points and classification of the dense cloud into terrain points;
• Generation of the final model using terrain points (DTM) and quality reduction from 1.5 cm/pixel to 20 cm/pixel;
• Generation of contour lines (50 cm interval);
• Creation of the textured 3D topographic model;
• Creation of the orthomosaic (1.15 cm/pixel).

3.2. Geophysical Methods

The geophysical surveys in the study area were implemented as depicted in Figure 4, including single-station ambient vibrations recordings for HVSR, active seismic (surface-wave surveys) and geoelectrical profiles for tomographic 2D and 3D imaging (ERT). The location and the length of the seismic and resistivity profiles have been fixed after the early field surveys. The eastern end of the investigated area includes the location of the Olpetella geothermal well, from which information on the expected stratigraphic setting has been derived. Seismic surveys were extended to the Olpetella well area to possibly extend stratigraphic information in the study area through geophysical models.

3.2.1. Electrical Resistivity Surveys

This method is recognized as the most suitable electrical method to image the subsoil by exploiting the correlation of electrical resistivity with hydraulic and geotechnical parameters [41]. Investigations are routinely conducted from the soil surface by using a four-pole system according to which the soil is energized with a known amount of electric current by the application of a potential difference between two metal pins conventionally defined as current electrodes. As a result of the induced current, a potential field is generated within the ground, the intensity of which is measured by means of a second pair of electrodes (potential electrodes). The four-electrode configurations used for resistivity measurements follow field arrangements called “array”, e.g., Wenner (W), Wenner–Schlumberger (WS), Dipole–Dipole (DD), etc. (see, among others, [42]), for which the main characteristics of resolution, sensitivity, and depth of investigation are known [42,43]. Subsoil resistivity distribution along a 2D section is estimated through the use of geophysical software, by means of the so-called “inversion” procedure of the experimental data which proceeds through successive iterations to search for the distribution of resistivity values that minimize the differences between measured and calculated resistivity values, within acceptable error limits (usually expressed as root mean square error, RMS). The area affected by the sinkhole was investigated through four resistivity profiles, crossing the study area to possibly model the subsoil setting along the NS and EW directions. Resistivity data were also acquired in a real 3D configuration, using L-shaped segments of cables deployed around the sinkhole, to achieve resistivity measurements referring to the subsurface of the area now occupied by the hollow. The MAE X612EM+ multichannel geo-resistivimeter (www.mae-srl.it, accessed on 10 November 2023), a multi-electrode system with an integrated P.C. capable of handling up to 96 electrodes, was used for geoelectrical prospecting. Dipole–Dipole (DD) and Wenner–Schlumberger (WS) arrays were used for the characterization of the area along all profiles to combine the needs of horizontal coverage, vertical and lateral resolution, and depth of investigation. Based on the available space, a maximum of 72 electrodes with 3 m spacing were deployed along each profile or L-shaped ground configuration around the sinkhole. Data referring to the 2D profiles were processed with the 2D inversion software Res2DInv [43,44]. As a subsequent step, the 3D resistivity modeling referred to the subsoil of a squared area (about 45 × 45 m²), approximately centered on the sinkhole lake, up to 10 m in depth. Data were processed by means of the ViewLab software (Geostudi Astier srl (Livorno, Italy), Multi-Phase Technologies LLC (Sparks, NV, USA)) that implements the Occam’s regularization as described in [45,46]. The inversion algorithm is efficient in managing noisy data, as the robust inversion procedure allows a data re-weighting for the adaptive changes in the variance matrix after each iteration (see [47] for further details on this topic). A value of Gaussian noise equal to 5% (i.e., the error associated with measurements, in percentage) was chosen before proceeding to the inversion, as well as the initial apparent resistivity value for the starting homogeneous half-space, chosen to be equal to 20 Ωm (corresponding to the mean apparent resistivity value).

3.2.2. Surface-Wave Surveys

Seismic methods have been widely used for many decades for subsurface characterization in support of exploration and geotechnical design due to the link between seismic parameters, such as the compressional (P) wave velocity, shear wave (S) velocity, signal attenuation, energy of reflections/reflections, etc., and the resistance of materials to normal and tangential stresses [48]. Among the seismic methods, the surface waves can be efficiently applied for characterizing soils in terms of shear wave velocity [49]. Surface waves are called “dispersive” because their propagation in inhomogeneous media occurs according to a phase velocity that varies with frequency [50,51,52]. Data acquisition can be performed with seismic survey equipment by employing an array of receivers (for phase velocity estimation, as in [51,52]) or a single sensor (as for the case of amplitude and group velocity analysis, e.g., [53]). According to the acquisition procedure and analysis of a multichannel seismic record proposed by [52], i.e., the so-called Multi-Channel Analysis of Surface Waves (MASW), a series of regularly spaced vertical sensors, usually 4.5 Hz geophones, are deployed along a straight line and a vertical impact source is used to generate Rayleigh waves. The recorded seismic traces are then transformed into the velocity–frequency domain to obtain the (Rayleigh or Love) phase velocity spectrum. Active seismic sources (e.g., sledgehammer or drop weight) are mostly used to generate artificial surface waves, but microtremors (e.g., [54,55]) are also exploited as passive seismic sources for site investigation due to their energy content at very low frequency, thus allowing us to extend the velocity profile to depth. In the classical approach, the velocity spectrum is interpreted in terms of modal dispersion curves that are inverted to eventually obtain the VS vertical profile [51]. For the sites under investigation, the data were processed according to the Full Velocity Spectrum approach (FVS, see [56]) through the WinMASW Academy 2019 software (www.eliosoft.com, accessed on 10 November 2023), whose modeling/inversion genetic algorithm consists of a process of progressive optimization of the search for models that best fit the experimental data via the Pareto front analysis [57,58], represented by one or more objective functions (e.g., single dispersion (R or L); Rayleigh dispersion + Love dispersion; Love dispersion + HVSR, etc.). The data considered in the modeling process for the investigated sites are summarized in

Table 1. Measurement parameters adopted, and surface-wave components recorded. Specification of the sensors used (H = horizontal geophone; V = vertical geophone) and the seismic phases sampled named according to [57] (ZVF: Rayleigh vertical component; RVF: Rayleigh radial component; THF: Love phase) used to estimate velocity profiles.

Site	Source	Impact	Geoph.	SW Comp.	Recording Time (s)	Samp. Rate (ms)
MASW1	Active	V	H, V	ZVF, RVF	2	0.125
MASW2	Active	V, H	H, V	ZVF, THF	2	0.125
MASW3	Active	V	V	ZVF	2	0.125
ESAC/MASW	Active	V	V	ZVF	2	0.125
ESAC/MASW	Passive	V	V	Z	420	8

.

The active seismic surveys at the ESAC/MASW site were carried out using 16 vertical geophones (natural frequency 4.5 Hz) equally spaced 3 m, setting the off-end shots at a distance equal to 5 and 10 m, in a mutual configuration. Data at the ESAC/MASW site were acquired using a 16-channel Sara Instruments seismograph connected to 16 vertical geophones (4.5 Hz natural frequency). The active seismic data were coupled via passive seismic recordings, which were processed according to the ESAC technique [59]. At the MASW 1 and MASW2 sites, acquisitions were made by using the “Tromino” equipment (www.moho.it, accessed on 10 November 2023), a 3-Component velocimeter capable of recording seismic signals along the horizontal plane in two directions (one parallel and the other perpendicular to the spreading for Love waves) and along the vertical plane. The minimum offset distance of 10 m was adopted, simulating a multichannel array using one receiver [60] by conducting multiple shots at increasing distances from the 3-C receiver. The dataset for these two sites was subsequently recomposed to perform the multichannel analysis.

3.2.3. Ambient Vibration Recordings

For many decades, the Horizontal to Vertical Spectral Ratio (HVSR) method popularized by [40], has been widely used for estimating the resonance frequency of a site. The HVSR method relies on the analysis of the spectral ratio between the horizontal (H) and the vertical (V) components of ground motion originated by ambient vibrations and recorded using a three-component (accelerometers or velocimeters) single station. According to many authors, the frequency corresponding to the maximum value of the HVSR function was shown to have a strict correspondence with the local resonance frequency f0 of a site (see, among many others, [61,62,63,64,65]). The HVSR curve may exhibit several peaks, caused by the presence at depth of sedimentary layers contrasting in lithology and, thus, also in acoustic impedance (i.e., the product of geologic material density and its seismic velocity). Therefore, in principle, by applying the HVSR methods, it is possible to constrain the local seismo-stratigraphical setting. Assuming the hypothesis of a 1D planar-layered stratigraphy, e.g., a soft sedimentary cover overlying a rigid (rocky) bedrock, the resonance phenomenon results from the entrapment of seismic waves in the loose covers. The resonance of the sedimentary layer occurs at the frequency corresponding to the ratio between the average shear wave velocity (VS) and its thickness (H) according to the following equation (see [61]):

f₀ = VS/4H

(1)

Moreover, it can be seen that if the subsoil average seismic properties (VS) are known, it is possible to reconstruct a simplified subsurface model by estimating the thickness of the resonant layer from (1). In the study area, single-station ambient recordings were collected by employing the Tromino instrument described before. The duration of the noise recordings ranged from 20 to 220 min, and the sampling frequency was fixed at 128 Hz.

3.3. Geochemical Surveys

Natural radioactivity, originating from nucleosynthesis and biogeochemical cycling, is widespread in the lithosphere and various environmental compartments. It constitutes a significant part of the ambient radiation dose for humans, mainly from primordial radionuclides like 40K and the decay series ²³⁸U, ²³²Th, and ²³⁵U. Each of these decay series in-corporates an intermediate radon isotope, distinguished by its noble gas properties, which render it prone to facile mobilization from the mineral matrices where it originates through radioactive decay. The naturally occurring radon isotopes encompass ²²²Rn (with a half-life of t1/2 = 3.8 days), ²²⁰Rn (t1/2 = 55.6 s), and ²¹⁹Rn (t1/2 = 3.96 s), arranged in decreasing order of half-life. Of these, only ²²²Rn possesses a half-life of sufficient duration to hold geological significance. Indeed, in contrast to ²²⁰Rn, the source of ²²²Rn can be attributed not solely to local contributions from shallow soil and rocks but may also result from advection processes, i.e., being transported to the surface from specific depths [65,66,67,68]. The measurement of soil gases and Terrestrial Gamma Dose Rate (TGDR) is of fundamental importance in environmental monitoring and understanding the natural and anthropogenic processes that affect air composition and ionizing radiation levels in the surrounding environment [65,69,70,71]. Over time, various measurement techniques have been developed and successfully applied to collect reliable and accurate data. Previous geochemical studies of the Latera caldera, aimed at carbon dioxide storage and geo-thermal exploration, have revealed the existence of near-surface deep thermogenic fluids, with a specific focus on pronounced CO₂ anomalies [72,73]. The sampled parameters included soil gas (CO₂, O₂, ²²²Rn, ²²⁰Rn) and terrestrial gamma dose rate (TGDR). Regarding soil gas sampling, this well-defined technique involves collecting soil air samples by inserting a 6.4 mm thick stainless-steel probe into the ground to a depth ranging from 0.7 to 0.9 m to minimize the influence of meteorological variables. A portable gas analyzer (Draeger Xam 7000) was directly connected to the probe to measure CO₂, CH₄, O₂, H₂S, and H₂ concentrations in situ. Furthermore, no CH₄, H₂S, or H₂ were detected. Radon (²²²Rn) and thoron (²²⁰Rn) activities in the soil were measured using a portable alpha detector RAD7 (Durridge Company Inc., Sheffield, UK) connected to the sampling probe through a drying tube used to maintain the instrument’s relative humidity below 10%. The alpha detector consists of an ion-implanted silicon semiconductor calibrated to measure in the range of 4 to 400,000 Bq/m³. A single measurement had an average duration of 20 min, with partial readings taken every 5 min. As for TGDR measurements, measurements were conducted using a portable sodium iodide (NaI) spectrometer, known as the “nuke ALERT II Model 1703M,” capable of directly measuring the Ambient Dose Equivalent Rate (ADER). The results were corrected using a calibration coefficient, with an error of approximately 5% and an overall uncertainty of 8%. Field activities were conducted from 8 to 11 August 2023, and included in situ measurements of ²²²Rn, ²²⁰Rn, CO₂, CH₄, O₂, H₂S, and H₂ concentrations in the soil, as well as gamma radiation measurements at each sampling site. The geochemical survey was extended to approximately 30,000 square meters around the sinkhole by means of 89 measurements points. Measurements were made during dry weather conditions with stable atmospheric pressure, minimal or absent wind, and the early morning hours were avoided to prevent the effect of radon progeny accumulation in dew, which could increase the normal soil emission rate by up to 15%. To characterize the population of soil gas and TGDR samples, descriptive statistics and graphical representations were employed. In order to examine the spatial variability of the measured parameters, data were processed into contour maps using a geostatistical approach, which includes an analysis of the experimental variogram and the construction of contour maps through ordinary kriging. Anisotropy was determined for kriging estimation, as spatial correlation can vary depending on the direction in which experimental variograms are calculated. The experimental variograms thus obtained were modeled using common spatial correlation models, and the overall error of each model (root mean square error, RMSE, and/or mean standard error, MSE) was calculated using the leave-one-out cross-validation technique. The parameters of the selected model were utilized to estimate soil gas concentrations at unsampled locations using ordinary kriging as geostatistical interpolator. All resulting maps have a grid size of 5 × 5 m², representing one-quarter of the average point measurement distance.

4. Results

4.1. Photogrammetry

The Digital Terrain Model (DTM) obtained from the survey conducted at the end of March can be observed in Figure 5. It is used here to describe the overall morphological setting of the study area. Two distinct small, depressed areas were identified within a larger one, delimited westward and northward by the 447.5 contour line. The sinkhole occupies the first depression that, based on the 3D topographic model of 1 February, is elliptical in shape with a surface of about 695 m², with a major axis length of 34.5 m. The second depression, with a minimum altitude slightly lower than 446.5 m a.s.l., is located just northwest (Figure 4). Based on the monitoring activity conducted so far, since the end of February 2023, no significant morphological variations were detected by comparing subsequent 3D digital terrain models.

For safety reasons, it was impossible to approach the hollow edge to perform direct water level measurements; thus, we exploited the UAV surveys to retrieve information about the water level variation with respect to the ground surface. In

Table 2. Water depth in the sinkhole (meters from ground surface).

Survey Date (Month/Day/Year)	Water Level [m]
02/01/2023	3.0
02/25/2023	0.4
03/30/2023	0.4
04/27/2023	0.4
06/09/2023	1.0
08/09/2023	2.7
09/18/2023	3.3

, the estimates of the water level inside the sinkhole are reported.

The initial size and depth of the sinkhole were not measured; the hollow was almost completely filled by groundwater in about eight days, between 31 January and 7 February, as reported by a few eyewitnesses. Based on the differential measurements from 1 February, i.e., when the water level was measured 3 m below the ground level, to 25 February, an estimated additional volume of water equal to about 1900 m³ filled the cavity as the water level reached 0.4 m from the ground surface. During the same time period, we detected a slight enlargement of the sinkhole in the direction of Provincial Road n.118 due to some minor slips along the banks (compare aerial Figure 6A,B); thus, the total surface reached 720 m². The subsequent measurements on 30 March and 27 April did not indicate any significant changes in water level and in sinkhole shape. Conversely, after the survey on 9 June, we determined a decrease in the water level of approximately 0.65 m (Figure 6C), corresponding to a water volume loss of about 470 m³, despite the days between May and the beginning of June being characterized by intense rainfall events (see Table 2). The survey of 9 August showed a decrease in the water level that was measured at about 2.7 m from the ground level (similar to that recorded on 1 February), corresponding to approximately 1800 m³ of water being leaked from the cavity. From the last survey, conducted on 18 September (Figure 6D), the water level was measured at 3.3 m from the ground level, for a total loss of about 2000 m³ since the beginning of May. The beginning of the emptying phase was estimated around 12 May (about 104 days from the sinkhole’s formation) through a linear regression analysis of the measurements inferred from the UAV surveys (Figure 7).

To possibly identify the cause of water level fluctuations, monthly precipitations between February and September were also considered to estimate the water volume loss for evapotranspiration. The data in

Table 3. Data from the Latera rain gauge station and the estimated volume loss due to the evapotranspiration.

Month	Rainfall	PET	Volume Loss (m3)
February	28.8	38.3	26.81
March	67.4	61.6	43.12
April	59.4	89.3	62.51
May	87.2	107.3	75.11
June	75.6	138.9	97.23
July	10.2	37.2	26.04
August	43.3	157.9	110.53
September 1	12.6	88.5	61.95

¹ Only 20 days.

were retrieved from the ARSIAL meteorological station “Latera-Centro Florovivaistico”. Available data include the potential evapotranspiration (PET) values calculated according to the Penman formula [74] after [75]. Water volume loss due to evapotranspiration was estimated by considering a constant surface of about 700 m² for the sinkhole.

The ARSIAL database reports a monthly rain equal to 37.6 mm between the end of January and the end of February; this amount of rain was considered not to be enough to explain the infill of the sinkhole by surface water. On the other hand, the loss of water volume due to evapotranspiration estimated between May and September is equal to c.a. 371 m³, which is not able to explain the volumetric decrease of approximately 2000 m³ starting from May. Therefore, even without considering the rainfall, the decrease in the water level (still ongoing) was most likely due to changes in the groundwater flow.

4.2. Resistivity Surveys

The 2D models along the ERT1 and ERT3 profiles are presented in Figure 8 as representative of the resistivity values of the soils along the two main directions of investigation. In general, the association of DD and WS measurements resulted in better horizontal sensitivity (laterally at the model edges) than that could have been obtained with WS or DD data alone: the models were reliable and useful for characterizing the subsurface up to 30 m in depth (elevations of 425–430 m above the sea level). The 2D model for ERT1 is proposed in Figure 8A: within the first 3–4 m depth, a medium resistive surface layer (40–70 Ωm) corresponds to the superficial layer with fine-to-medium grain size. From about 5 m in depth, the resistivity drops to values below 15 Ωm almost evenly along two-thirds of the section.

These low values are compatible with the presence of fine grain size deposits (and/or reworked pyroclastic deposits), most likely saturated by groundwater. Below and laterally, the resistivity increases and therefore it is possible to assume a different lithological composition/a remarkably lower water content. To be noted, starting from the progressive x = 144 m in the SSE direction, there can be observed the presence of a resistivity gradient that marks the contact of the aforementioned deposits with a geological body characterized by resistivity even higher than 150 Ωm. These values are compatible with coherent volcanic deposits and/or fractured lavas with water. The model adopted for the ERT 3 profile is presented in Figure 9. The resistivity distribution clearly recalls that of the ERT1 profile in terms of the number and thickness of resistivity layers, buried geometries and values. Exceptions are the two sectors at the edges of the model, where two areas of high resistivity could indicate, similarly to what was observed for ERT1, the possible presence in depth of coherent volcanic deposits and/or fractured lavas with water. The final 3D resistivity distribution revealed an anomalous sector characterized by resistivity values significantly higher than the average of soils in the study area. Resistivity values suggested the volume could be empty or only partially filled with water. The longitudinal axis of the anomaly runs oriented NW-SE, intersecting the roadway, as shown in Figure 9.

4.3. Surface-Wave Surveys

In general terms, for all the sites, the inversion converged with an associated misfit of less than 10% and the velocity profiles were reliable for depths of several tens of meters. As shown in Figure 10, the seismic S-wave velocities close to the sinkhole do not exceed 200 m/s for nearly 20 m in depth, suggesting the presence of deposits with little stiffness for many meters in depth. Along the southern margin, near the Provincial roadway, the soil condition seems to be slightly better as seismic velocities appear to be higher, albeit by a small amount. Note incidentally that, based on [76], the time-averaged VS value to a depth of 30 m (VS30), i.e., thus including almost the whole thickness of the superficial covers, was expected to be 340 m/s on average, so that a strong impedance contrast was also possible, considering the presence of the lava unit just below.

The VS30 values are considerably smaller than expected, in particular moving from the MASW3 site close to the Olpetella well to MASW2 and ESAC/MASW sites near the sinkhole. In addition, based on the velocity values in the upper 5 m, the VS values barely reached 200 m/s.

4.4. Ambient Vibration Recordings

The procedures proposed in [62] have been implemented for the processing and the statistical validation of the Horizontal-to-Vertical Spectral Ratio (HVSR) curves. The data were analyzed and processed using the WinMASW Academy 2019 software (www.elisoft.it, accessed on 10 November 2023); to retrieve HVSR curves, each recording was first sub-divided into time windows of 40 s, selecting the most stationary part of the registered signals and excluding transient parts since these are potentially associated with close sources. Then, for each time window, the signal was cosine tapered (5%), and the Fourier spectrum was calculated within the 0.20–20 Hz frequency interval. Smoothing was applied using a triangular window at a width equal to 10% of the central frequency. The single-component spectra (NS, EW and Z) are the averages resulting from all analyzed windows; the final curve for each recording site was then calculated between the two horizontal components and the vertical one according to Nakamura’s procedure [40]. As an example, Figure 11 shows the data processed in terms of HVSR curves from noise recordings placed east and west of the sinkhole so as to highlight the marked shift toward low peak frequencies in the same direction. The spectral ratio curves have been further processed in order to obtain 2D profiles of the HVSR values [77,78,79] by considering the average value of VS equal to 150 m/s. Along these two sections, oriented WSW-ENE (Figure 12A) and NNW-SSE (Figure 12B), the upper color scale indicates positive acoustic impedance contrast, while cool colors are associated with low-to-intermediate contrasts that could correspond to an absence of seismic velocity changes with depth.

The 2D distribution of HVSR values highlights a progressive deepening of the impedance contrasts from SSE to NNW, while an even sharper deepening can be observed from SSE to NNW very close to the sinkhole (Figure 12C and Figure 12D, respectively); both sections suggest a sudden and unexpected increase in the thickness of the superficial deposits right in the center of the study area up to 25–30 m, confirming what was observed along the ERTs.

4.5. Geochemical Survey

The investigated area has been centered on the sinkhole, including the whole sector previously investigated via the ERTs (Figure 13); thus, the results have been integrated to obtain a more comprehensive and reliable subsoil model. The geometry of the sampling point distribution and the distance between samples were chosen considering the extent of the sinkhole and the accessibility of the area. The distribution of sampling points is not regular due to dense vegetation in certain areas. The sampling points were spaced approximately 10 m apart near the sinkhole and about 20 m apart moving away from it. Additional sampling points were measured far from the sinkhole to assess the general background of the caldera area.

In

Table 4. Main indexes of the studied variables (N: number of samples; Std. Dev: standard deviation).

Parameter	N	Mean	Median	Min	Max	Std. Dev.
CO2 vol%	89	2.8	2.8	1.0	5.4	1.0
O2 vol%	89	19.6	19.7	18.3	20.9	0.6
TGDR μSv/h	89	0.218	0.217	0.187	0.261	0.014
222Rn kBq/m3	89	173.2	161.7	62.9	303.1	54.2
220Rn kBq/m3	89	399.4	349.5	102.0	862.8	142.4

, the main statistical indexes for the measured variables are reported.

Specifically, in relation to the range and maximum values (minimum for O₂) of CO₂ and O₂, no anomalous values have been observed beyond what is considered normal due to biological activity and the morphology of the area. The TGDR (despite having a limited range) exhibits elevated yet normal values for a volcanic area, which can be attributed to the geological characteristics of the lithotypes in the study area. As for the soil activity of ²²²Rn and ²²⁰Rn, anomalous values have been detected, suggesting a deep-seated source for ²²²Rn and a shallower origin, likely related to high-permeability zones for ²²⁰Rn [80,81,82]. Regarding ²²²Rn, a comparative analysis with other studies realized within the Vulsini volcanic complex immediately highlights the significantly elevated values uncovered in this study. For instance, within the municipality of Celleno (Vt), [70] recorded an average ²²²Rn concentration of 60 kBq/m³, with a maximum value of 253 kBq/m³ (vs. 303.1 kBq/m³ found in this study). Analyzing the spatial distribution of CO₂ in detail, as shown in Figure 14A, concentration anomalies are evident to the west and north of the sinkhole. These anomalous values, with a maximum concentration of 5.4 vol%, fall within the range attributable to soil biological activity (1–10 vol%), according to [83]. The distribution of anomalies is likely correlated with the presence of morphological depressions with dense vegetation and moist soil. In Figure 14B, a map of the spatial distribution of TGDR is depicted. Several small areas with positive anomalies are observed to the north and south of the sinkhole, as well as to the east of it. TGDR is related to the content of radionuclides in rocks, and negative anomalies in this radiation should indicate the presence of less dense areas, such as cavities, as well as faults and fractures. The negative anomaly just east of the sinkhole, parallel to the ravine, could suggest the presence of a structure with a north–south direction. In addition, the discontinuity among the positive anomalous values found in the sinkhole’s southeast sector is compatible with the presence of the sub-horizontal conduit (see Figure 9) inferred from the 3D resistivity model. Figure 14C depicts the spatial distribution of ²²²Rn. There are primarily two areas of positive anomalies located to the south and northeast of the sinkhole. On the other hand, the area to the west of the sinkhole exhibits lower values of ²²²Rn, confirming the surface biological production mentioned for CO₂. Regarding ²²⁰Rn, in Figure 14D, there is a large positive anomalous zone at the easternmost part of the study area, while both around the sinkhole and westward, the values are much lower. The presence of this anomalous zone could be related to surface factors such as soil permeability in the upper meters and water content. In this area, dry and superficially fractured soils (mud cracks effect) were found that could promote the upward migration of ²²⁰Rn in the last meters.

Finally, in Figure 15, the ²²²Rn/²²⁰Rn ratio is represented. This ratio highlights the areas with a deeper source of ²²²Rn (values ranging from 0.48 to 0.72). Indeed, only ²²²Rn has a sufficiently long half-life (3.8 days) to originate from greater depths together with a surface origin. This can only occur through advection transport, using a carrier gas such as CO₂ and the presence of sufficiently permeable fractures/faults to allow for a rapid ascent [65,66,67]. It is evident, examining the spatial distribution map of ²²²Rn/²²⁰Rn, that the main positive anomalous zone is located east of the sinkhole in a north–south direction. This leads us to hypothesize the presence of a discontinuity, also supported by CO₂, which acts as a carrier gas for ²²²Rn, also found in high concentrations along the ditch.

5. Discussion

Since its formation on 31 January 2023, the sinkhole has been monitored on an airborne surveys basis, with attention focused on its morphological evolution (shape and size) as it opened very close to a provincial road. In this regard, no significant changes in the morphological setting have been registered during the monitoring period (still active), nor along the Provincial Road located very close to the hollow, or across the whole study area. Minor changes only occurred in the shape of the sinkhole due to the occurrence of soil slips along the banks; as a consequence, its elliptical shape was slightly modified with an increment in the total surface of less than 5% (from 695 m² to 720 m²). UAV data were also processed to retrieve indirect measurements of the water level within the sinkhole; it varied sensibly across the months and this was most likely due to the groundwater flow changes, indicating that, until present, the fluid circulation beneath the sinkhole is far from stable. The analysis of orthophotos from 1984 to 2021 referring to the study area (not shown here) has allowed for the identification of a sub-circular sinkhole paleoform in the same site. The paleoform was already visible in 1984 and has become progressively more evident over the years. It was particularly pronounced in 2015. Therefore, it is possible to consider the formed sinkhole as a reactivation phenomenon. Additionally, information collected regarding the local traditions of the study area has highlighted the narrative of a sinkhole transforming into a small lake in the nineteenth century. Non-destructive methods were also applied to investigate the subsoil. Based on the indirect methods, the sinkhole opened inside a small basin filled with fine grain-size deposits, with a thickness of several tens of meters at least (to be constrained by using boreholes). The edges of the basin were inferred along two main directions investigated via ERTs from the presence of resistivity gradients that mark the contact of the conductive deposits with geological units that, based on their higher resistivity (up to 150 Ωm and more), were interpreted as coherent volcanic deposits and/or fractured, saturated lavas. In particular, due to effusive deposits at different depths, such depth intervals should correspond to high seismic velocity layers and medium-to-high resistivity ones, at least within the first tens meters in depth (i.e., above the groundwater level). The HVSR results supported the preliminary hypothesis of a strong impedance contrast between the covers and underlying rigid layer (lava units based on the well-log), and the frequency shift from 4.5 Hz–8 Hz to 1 Hz–2 Hz (Figure 11) approaching the sinkhole from NE was interpreted as the effect of the deepening of the main impedance contrast. As highlighted by the 2D HVSR sections, the increase in the soft covers occurred very abruptly, mostly in the north–south direction; thus, the presence of a discontinuity (fracture/fault zone) could be inferred very near the sinkhole. In this context, geochemical investigations complemented the geophysical data, as shown by considering the respective spatial distributions. Positive anomalies of TGDR and the ²²²Rn/²²⁰Rn map support the hypothesis of a buried discontinuity running along the NS direction, responsible for the abrupt increase in loose deposits just in the center of the investigated area. This seems to indicate the presence of a north–south oriented discontinuity located just east of the sinkhole (without, however, the definition of a clear geometry). The ²²²Rn and ²²⁰Rn distributions indicate a possible lithological change between the western sector and the eastern sector within the study area. Differences in average shear-wave velocity values approaching the sinkhole area can be explained considering (i) the possible presence of terrains softer than expected; and (ii) a non-negligible modification in the strength parameters of terrains close to the sinkhole area, possibly induced by the collapse. The main contribution of the ERT to the subsoil model was the identification of the 3D anomaly located beneath the sinkhole that runs below the Provincial roadway. Incidentally, it can be noted that this method has proven to be the most cost effective among those implemented here since it is able to provide reliable 2D and 3D information about the subsoil. On the other hand, the constraints introduced into the subsurface model through the seismic and geochemical methods, which require more processing, proved essential for effective characterization of the area. Considering the hazardous conditions, it is also worth mentioning the smart use of the UAV data for monitoring the water level inside the hollow as a proxy parameter to the possible evolution of the phenomenon.

6. Conclusions

Airborne surveys helped in periodically monitoring the morphological setting of the study area, also allowing for the indirect measurements of the water level inside the cavity. In this regard, after the initial filling of the cavity by groundwater, starting from May 2023, a clear groundwater discharge was detected and monitored. Based on the data already collected, at present, there is concern that the sinkhole may expand laterally and potentially affect the Provincial Road, creating further issues. The geophysical and geochemical investigations conducted so far allowed for a preliminary geotechnical characterization of soils and the image of the main targets, consisting of discontinuities and/or buried structures, potentially responsible for the collapse at the surface. The presence of a paleoform in the sinkhole area since 1984, along with the legends recounted by the local inhabitants about the sudden formation of a cavity and the creation of a lake in the nineteenth century, demonstrate that the phenomenon under investigation is a natural sinkhole. Furthermore, given the current low level of urbanization in the area and its rural character in the nineteenth century, it is possible to exclude a triggering cause related to human activities (such as mining or water regulation). An erosive mechanism can be proposed for the phenomenon under evaluation. At present, the subsoil model accounts for the presence of a buried conduit that could have played a key role in the processes of erosion and syphoning of the covering materials, leading to the final collapse of the superficial covers. The predisposing causes of the phenomenon must be sought in the geological–structural and hydrogeological context of the area: an intra-calderic volcanic plain characterized by the circulation of gaseous fluids and pressurized mineralized waters. Furthermore, the presence of a karstified carbonate basement is excluded. Additional investigations will be needed in order to better define the process that led to the formation of the sinkhole as the deep piping process cannot be excluded, due to the presence of the discontinuity revealed via indirect methods that could drive the leakage of fluids from the collapsed area. As an overall consideration, the geophysical models, the geochemical evidence, and UAV-derived surface data highlighted the significant complexity of the area where the sinkhole has developed, and can motivate local authorities to carry out direct investigations to (i) improve the geology at a local scale; (ii) calibrate the geophysical models; (iii) ascertain the subsurface conditions under the roadway, in correspondence to the high-resistivity anomaly; and (iv) collect the descriptive parameters of the mechanical behavior of the soils to properly assess the residual risk for the Provincial Road n.118, as the sinkhole structure appears to be unstable due to the water discharge (still on course) at present and to soil slips that have occurred along the banks.