Multidisciplinary Study of Subsidence and Sinkhole Occurrences in the Acque Albule Basin (Roma, Italy)

We present the results of a combined analysis of remote sensing and geophysical‐geotechnical data carried out in the Acque Albule Basin, a sinkhole prone area located close to the city of Roma, where a wide travertine wedge is present. We carried out geophysical measurements and borehole drillings over two test areas to image the subsoil where paroxysmal surficial dynamics occur. One site is marked by subsidence occurring at least since the early 2000s, whereas the other site hosts the “La Regina” and “Colonnelle” sinkhole lakes, which discharge sulfur‐carbonated waters. The stability of these two sites threatens highway, railway, and airport facilities, and this study helps to assess the geological hazard. For example, InSAR and LiDAR data helped define wide scale subsidence over the last 20 years and previously undetected small‐scale morphologies. Geophysical measurements of the latter revealed shallow and deep dissolution affecting the travertine and driving surficial paroxysmal events. Both study sites were found to lie inside a large depression located at the junction between Jurassic carbonate and Plio‐Pleistocene units in association with paleo karst morphologies in the travertine deposits and affected by the present‐past spillage of sulfurous waters. Given these elements, multidisciplinary geophysical observations are crucial for assessing and mitigating the geological risk and guiding land use planning and management.


Introduction
Among the natural catastrophes afflicting the Italian peninsula, sudden and slow ground sinking (or subsidence and sinkholes) are common phenomena. Sinkholes mostly occur because of the interaction between water and carbonate bedrock, (e.g., Mesozoic limestone and quaternary travertine), salt, and gypsum. Sinkhole risk evaluation is particularly important for forecasting the development of new sinkholes and their rejuvenation processes. In this study, we investigate the Acque Albule Basin (AAB), a wide sinkhole prone area  situated 20 km north-east of Rome, (Figures 1a and 1b). The AAB is an easily accessed highly populated area of roughly 37 km 2 . This area represents a unique laboratory for studying geological risk due to its extensive travertine mining, thermal bath, and ground subsidence involving a major highway and railway.  (Nanni et al., 2014), and the locations of Sites A and B outlined by the black rectangles. The AAB is characterized by hot springs of sulfurous waters known since Roman times and exploited for thermal purposes and its travertine deposits that constitute the main basin fill. It holds world renown travertine deposits from which most of Rome's monuments (e.g., the Coliseum and Saint Peter's Cathedral) were constructed. Currently, the AAB is a highly urbanized plain strongly affected by local ground subsidence and sinkholes. In the AAB, many hydrothermal springs and degassing phenomena are observed, although the underground circulation has been severely affected by lowering the water table to exploit the travertine deposits down to the depth of 30-40 meters (Figure 1c). Faccenna, Funiciello, Montone, et al., (1994 interpreted the AAB as a pull-apart basin located in an N-S trending shear zone controlled by N30°-40°striking transtensional to normal, right-lateral faults that may be located at depth on both sides of the basin, although active faulting has not been observed (Coltorti, pers. comm. Figure 1a). However, a major tectonic escarpment marks the western side of the Apennine fringe that constitutes the eastern border of the basin. The master tectonic element (called the "Sabina Fault") was active during the early Middle Pleistocene. The antithetic faults and fractures associated with the Sabina Fault may be located at a depth beneath the La Regina and Colonnelle Lakes area. These elements also affect the Plio-Pleistocene travertine plateau, where a N40°-striking set of joints has been observed (Faccenna, Funiciello, Montone, et al., 1994;Sagnotti et al., 1994). It is also possible that these secondary tectonic elements may be associated with ancient buried faults and represent a line of weakness that drives karstification, ground collapse, and the subsurface hydrology.
The thickness of the travertine deposits is approximately 80-90 m within the shear zone and decreases steeply westward and eastward to approximately 10 m De Filippis et al., 2013;Figure 1c). The geological, structural, and hydrogeological characteristics of the AAB render this area prone to a variety of karst and hyper-karst effects (e.g., collapsed sinkholes, Gunn, 2004;Sauro, 2007;Nisio et al., 2007;Klimchouk et al., 2017).
In this work, we present a multidisciplinary investigation of the geological and geomorphological predisposing features that help drive paroxysmal phenomena using remote sensing, geophysical, and borehole data. Synthetic-aperture Radar (SAR) and the interferometric synthetic-aperture radar (InSAR) were used to map vertical ground movements, and Light Detection Ranging (LiDAR) data allowed us to construct a Digital Elevation Model (DEM) of the area. The remote sensing data cover the northern and western AAB, including the towns of Guidonia-Montecelio, Villalba, and Bagni di Tivoli, which have long subsidence histories (Salvi et al., 2004(Salvi et al., , 2005. Field measurements were carried out in the northern and northwestern sectors of the AAB (Site A and Site B in Figure 1c). These sites were chosen based on the remote sensing and their strong potential for geomorphological risks (e.g., ground collapse and subsiding phenomena, Nanni et al., 2014).
Site A is farmland affected by strong local subsidence, which apparently began in the early 2000s and suddenly accelerated in 2012. At Site A, preliminary geophysical and geological studies were conducted between 2012 and 2013 (Argentieri et al., 2015) that revealed a buried cavity in the travertine plateau for further investigation. Comparable geomorphological processes affect Site B, which is located on agricultural land approximately 1.8 km southwest of Site A (Figure 1c). In fact, the sinkholes holding the La Regina and Colonnelle Lakes resulted from resurgent geothermal activity in pre-Roman Age times . Therefore, considering the strong similarities between the two sites, the present study investigates and compares their surface and subsurface geological and geomorphological features and processes.

Geological Setting
The AAB is a depression developing in the NS and EW directions. The basin is bordered by the carbonate and silicate-marl Lucretili and Cornicolani ranges to the north and east, the carbonate and detrital carbonate Tiburtini-Prenestini range to the east, and the Colli Albani volcanic district to the south (Cosentino & Parotto, 1986). To the west, the AAB has a hilly landscape developed on Pleistocene coastal clays, sands, and conglomerates, and alluvial Quaternary deposits.
In the south-central part of the basin, the travertine topography is interpolated from borehole stratigraphy (Faccenna et al., 2008). Along the NS direction, the travertine deepens northward and towards the 10.1029/2019EA000870

Earth and Space Science
DE RITIS ET AL. south-central basin. Along the EW direction, a major depression occurs close to the western margin, and a secondary depression is in the central basin. The travertine features have a V-like geometry that suggests the existence of buried paleo-valleys belonging to the pre-travertine landscape (Mazza & La Vigna, 2011). However, the origin of these depressions has been associated with faults that did not seem to affect the base of the travertine (Faccenna, Funiciello, Montone, et al. (1994)  ). Nevertheless, given that the N-S karstic depressions may be following carbonate bedrock weakness zones, the genesis of these depressions continues to be debated.
The geological setting of the basin is related to the tectonic history of the Central Apennine fold-thrust belt, which underwent since Miocene time a sequence of compressive (Tortonian-early Pliocene) and extensional (mid Pliocene-present) tectonic phases (Mattei et al., 1986, Acocella & Funiciello et al., 2006, Faccenna, Funiciello, Montone, et al., 1994. During the Pliocene, marine ingression-regression phases affected the Tyrrhenian margin, and a thick terrigenous sequence of marine and coastal units was deposited over the Meso-Cenozoic carbonate bedrock. This bedrock is approximately 300-400 m deep on average in the AAB and becomes shallower in its northwestern sector (Maxia, 1950).
During the middle-upper Pleistocene, the extensional tectonic regime determined the onset of vulcanism that produced the Mt. Sabatini and the Albano Hills (ca. 700-800 kyr ago) and the deposition of thick pyroclastic cover. The last phase of the Albano Hills activity (170-180 kyr, Giordano et al., 2006) was characterized by hydromagmatic eruptions. Faccenna et al. (2008) suggested that the conspicuous inflow of thermal mineral water may have added to the circulation in the carbonate system generating the travertine formation, which occurred between 120 and 30 kyr ago. In fact, Annunziatellis et al. (2009), Mazza and La Vigna (2011), and Di Salvo et al. (2013 suggested that the travertines may have resulted from precipitation after the degassing of the sulfur carbonic waters with the mixed salt content of bicarbonate and alkaline sulphate upwelling from the Alban volcanism's deep geothermal circulation. Radiometric data, the surface geometry, and the reconstructed isopach indicate that travertine deposition started in the northern sector of the AAB, moved then southwards, and interrupted by erosive episodes (Faccenna et al., 2008;Funiciello et al., 2005). The basin's depocenter and resurgences shifted similarly. However, considerable time passed between the early Middle Pleistocene onset of volcanic activity and the travertine formation that occurred during the last Interglacial-Glacial cycle.
The AAB's stratigraphy from top to bottom includes Holocene flood deposits of the Aniene River and its tributaries; travertine formations; Albano pyroclastic deposits; Plio-Pleistocenic coastal clays, sands and conglomerates; and Mesozoic carbonate bedrock ( Figure 1a). The top of the travertine deposits ranges from 80-85 m asl (northern sector) to 40-45 m asl along the Aniene river (southern sector).
Hydrogeologically, the basin hosts two main aquifers with water flowing southwards and eastwards. The deeper circulation lies inside the Mesozoic carbonate bedrock in hydraulic continuity with the Cornicolani-Lucretili and Tiburtini-Prenestini ranges. However, the shallower aquifer lies inside the travertine body. The two aquifers, which are separated by the Pliocene clay aquiclude, appear to be connected through resurgences and normal-transcurrent tectonic elements, which have disarticulated the Mesozoic series (Di Salvo et al., 2013;Faccenna et al., 2008).
Limited geological information regarding Site A comes from the 26-m deep borehole located approximately 400 m southeastwards, where and the lithoid travertine occupies roughly the borehole's bottom 12 m (NT1 in Figures 2a and 2b), Funiciello et al., 2005). Moreover, the isopach trend shows an eastward thickening of the travertine (Funiciello et al., 2005;De Filippis et al., 2013), and a few hundred meters to the east, the travertine sequence some 30 m thick is exposed in quarries showing a slight clino-stratification towards the southwest (Figure 2c). The collected information suggests that the travertine thickens and deepens towards the quarries.
At Site B, the complex subsurface features between the La Regina-Colonnelle and S. Giovanni Lakes were modeled by Annunziatellis et al. (2009) on the base of nine boreholes up to 60 m deep and of the geological survey of the area. The results for the western sector suggested that travertine is approximately 5-20 m thick and rests on pyroclastic units approximately 40-50 m thick (SG-5 in Figures 2d and 2e) related to the Albani Mts. In the eastern sector, the SG-10 borehole intercepted some 40 of travertine but no pyroclastic units.

Data and Methods
Different methodologies were applied to discern the possible geophysical signatures of the surface and subsurface geological processes common to both study sites. The sections below further describe the details of the remote sensing and geophysical analyses constrained by the borehole data.

InSAR
Interferometric Synthetic Aperture Radar (InSAR) is based on the phase correlations between two images acquired at different times over the same area. This analysis provides topographical and ground movement estimates along the sensor's Line Of Sight (LOS) by exploiting the interferogram of the phase difference in each homologous pixel between two the images. The potential of SAR data is fully exploited using a large stack of images processed with multitemporal SAR interferometry techniques (commonly named A-InSAR). In this work, we adopt the Permanent Scatterers (PS) methodology as implemented by Ferretti et al. (2001).
To evaluate the AAB's ground deformations using the PS technique, three SAR images were considered as acquired by the European Space Agency's ERS1-2 and Envisat SAR sensors over the 1992-2000 and 2003-2010 periods, respectively, and the Italian Space Agency's COSMO-SkyMed constellation (CSK) over the March 2012-January 2016 period. Regarding the ERS and Envisat sensors, the LOS is the same and almost vertical (~23°), whereas for the CSK data, the LOS is~27°. Specifically, we analyzed 65 ERS and 46 Envisat images acquired along descending orbits and 50 Envisat images acquired along ascending orbits. We identified~71,000,~35,000, and~66,000 PS of good quality. Additionally, we considered 58 CSK images along the ascending orbit, which yielded more than 440,000 PS. As previously mentioned, most PS are located over buildings or urban feature with the PS-density that is function of the urbanization level. The PS were selected based on quality parameters related to the radar signal's stability. The InSAR outcomes include a ground velocity map and cumulated deformation time series (Figures 3 and 4). The displacements are related to reference points located outside the AAB, which are not affected by the subsidence phenomena affecting the AAB.
The displacements are mainly vertical and may include a small horizontal contribution (Fernández et al., 2018). However, the main focus of this study is to highlight the evidence of subsidence in the investigated areas. The Envisat data sets in the areas where the ascending and descending orbits overlap were decomposed into the E-W and vertical ground displacement components. Furthermore, the retrieved ground velocities along the LOS's from the adopted PS approach were mainly related to vertical movement. The InSAR analyses showed that several sites are strongly affected by subsidence phenomena.

LiDAR
LiDAR (Light Detection and Ranging) data provide high accuracy measurements of the distances to the ground from multiple returns of a laser beam. This method allows the rapid measurement of topography over large and even vegetated regions. Consequently, LiDAR data facilitate the detection and quantification of submetric relief to help characterize hydrogeomorphological features and tectonic processes (Meigs, 2013;Passalacqua et al., 2015).
The AAB area was surveyed by airborne LiDAR in 2008 by the Ministry of Environment and of Land and Sea Protection (MATTM). This data set has vertical and horizontal accuracies of approximately 15 and 30 cm, respectively. We used the LAStools software suite (n.d., https://rapidlasso.com/lastools/) to classify, tile, and filter the raw LiDAR-derived point cloud. Then, we created 1 m-pixel raster elevation models of the first returns of the ground's surface. A first returns' surface includes the tree canopy and buildings and is often referred to as the digital surface model (DSM). However, the ground devoid of the vegetation canopy and buildings is frequently called the digital terrain model (DTM). For the AAB analysis, we used the slopeshade image, which is a digital derivative map showing the slope steepness variations to simplify the perceiving and interpreting of the topographic features ( Figure 5).
The LiDAR (Light Detection and Ranging) data provide high accuracy measurements of the distances to the ground from multiple returns of a laser beam. This method allows the rapid measurement of topography over large and even vegetated regions. Consequently, LiDAR data facilitate the detection and quantification of submetric relief to help characterize hydrogeomorphological features and tectonic processes (Meigs, 2013;Passalacqua et al., 2015).

Borehole and Geotechnical Analyses
In September 2014, two boreholes were drilled at Site A to investigate the northern and eastern sides of the subsiding area (Figures 2a and 2b). Geotechnical analyses of the boreholes were conducted to help constrain the geophysical investigations. The vertical S1 borehole was located near the road on the northern side  beyond the subsiding area and equipped with an open pipe piezometer. The S2 borehole was located beyond the eastern side of the subsiding area and dipped approximately 30°to the west. Geotechnical laboratory analyses (Table S1 in the supplemental materials) on nine undisturbed samples were also carried out. The groundwater level in the S1 borehole was measured periodically during 2014-2016 ( Figure 2f). In addition, the reported results for borehole NT1 (Funiciello et al., 2005 in Figures 2a and 2b) and boreholes SG5 and SG10, which are near the lakes on the southwestern side of Site B (Annunziatellis et al., 2009 in Figures 2d and 2e), are also investigated below.

Gravimetric Prospecting
Gravimetric prospecting is a useful geophysical technique for investigating lateral density variations of the subsurface. Sinkholes originate from the dissolution of the travertine deposits by groundwater circulation. The cavities in the travertine deposits can be empty or filled with low-density materials to produce relatively negative local gravity anomalies. Thus, a gravimetric survey was carried out in the southern sector of Site A (Figures 6a and 6b) to complement the region investigated by Argentieri et al. (2015) and outline the spatial extent of the sinkholes. The microgravity network consisted of 53 stations spaced 25 meters apart on average covering an area of 250 × 150 m 2 to study the near-surface geology (e.g., Hinze et al., 2013). The gravity data were acquired using a LaCoste & Romberg model G gravimeter with a nominal accuracy of 0.01 mGal. The geographic positioning of the microgravimetric stations and their precise ellipsoidal and orthometric heights were acquired using a Leika Viva CS 15 GPS device (Leica Geosystems AG-Part of Hexagon). The real-time kinematic (RTK) data collection exhibited a vertical accuracy of approximately 15 mm.
The instrumental drift curve was evaluated by a pair of reference stations located in the survey area (purple triangles in Figure 6a) and linked to the absolute reference gravity station (001294-International Gravity Bureau database-BGI) in Rome from the International Gravity Standardization Net IGSN71 (Morelli et al., 1972). The simple Bouguer gravity anomaly map was computed after removing (i) the instrumental drift curve sampled at intervals of 1-2 hours (Mares & Tvrdý, 1984); (ii) the free-air correction of 308.6 μGal/m; (iii) the Bouguer slab effect computed for the regionally appropriate density of 1.9 g/cm 3 (Argentieri et al., 2015); and (iv) a first-order regional surface. Figures 6a and 6b show the residual simple Bouguer anomalies interpolated by the Kriging method onto a grid at an interval of 20 m.

Laboratory Measured Magnetizations
Sixty-eight discrete samples were collected at an approximately 1-m interval from the S2 borehole for analyses at the INGV paleomagnetic laboratory (Figures 2b and 7 and Figure S1 in the supplemental materials). The analyses were performed to investigate variations with depth in the concentration and coercivity of the magnetic minerals and to define the sediment's magnetic, particularly in the intervals with higher intensity

10.1029/2019EA000870
Earth and Space Science natural remanent magnetization (NRM) and magnetic susceptibility (k), which promote the production of magnetic anomalies. The NRM was measured using a 2G Enterprises (model 755) cryogenic magnetometer (Weeks et al., 1993) and stepwise demagnetized and measured at 10, 20, 30, 40, 50, 60, and 100 mT alternating field (AF) peaks. Low field magnetic susceptibility measurements were acquired using a Bartington magnetic susceptibility meter with a loop sensor (MS2C) in line with the magnetometer. An anhysteretic remanent magnetizations (ARM) were imparted by applying a 100 mT AF peak with a superimposed 0.05 mT bias. The ARM was stepwise demagnetized up to 40 mT and measured after each step.

Magnetic Prospecting
Magnetic anomalies develop mainly beneath the subsided or depressed sector D1 in Figures 5b and 5c within the poorly magnetized lithologies (e.g., travertine, clays, and clay silt). Magnetic prospecting was undertaken to constrain the possible lithologies, depths, and geometries of the magnetic anomaly sources to help provide insights on the D-depression's origins and development. The magnetic measurements were carried out using a Geometrics G858 scalar magnetometer at a sampling sensitivity of 0.01 nT. The magnetometer was also equipped with a GPS sensor for precisely estimating the intensity of the Earth's main magnetic field at each survey station. Magnetic survey Mag1 (Figures 5b, 6b, and 6c) consists of a 220 × 90 m 2 area using 70 survey lines striking NNE-SSW and spaced approximately 2 m apart ( Figure 6c). Survey Mag2 involves 20 survey lines striking WNW-ESE and spaced approximately 2 m apart that covers an area of about 300 × 200 m 2 .
The processing of the magnetic field measurements aims to isolate the anomalous magnetic effects of the near-surface (e.g., Hinze et al., 2013) by removing (1) the main geomagnetic field component estimated from the 2016 model of the international Geomagnetic Reference Field (IGRF2016 from http://www.noaa.gov); (2) diurnal magnetic field variations; and (3) the regional longer wavelength component. The resulting residual or Total Magnetic Intensity (TMI) anomaly was gridded to the anomaly map shown in Figure 6c. The TMI anomalies were reduced-to-the-pole (RTP) to mitigate the lateral displacement of the anomaly along geomagnetic declination from the source's subsurface position due to the non-vertical inclination of the main inducing geomagnetic field (Baranov & Naudy, 1964;Blakely, 1996). This technique aims to compute the anomaly field produced by the near-surface magnetic sources at the geomagnetic poles, where the field inclination is vertical, and the geomagnetic field induced anomalies are centered on their subsurface sources ( Figure 6d).

Magnetic Modeling
The main aim of magnetic modeling is to quantify the spatial properties of the causative sources and interpret them within the applicable geological framework. In this work, we implemented the Li and Oldenburg (2003) inversion method, to estimate the magnetic volume susceptibilities k for a 3D Voxel network of prismatic cells filling out the topography using the MGinv3D software available online (http://www. scicomap.com/MGinv3D.html, MGinv3D 1.24, 2013). The modeled subsurface volume was discretized with a core mesh of 2.5 × 11 × 4 m 3 prismatic cells in the respective x, y, and z directions with the bottom fixed at a depth of −60 m. The input data for the inversion were the RTP-TMI anomalies (Figure 6d), where each datum was voxel-centered (Walls & Hall, 1998, Shirman et al., 2014. Finally, the smoothing weights in the x, y, and z directions were set to unity, with the solution's target misfit obtained after 20 iterations.

2D and 3D Electrical Resistivity Tomography (ERT) 3.8.1. ERT 2D
The ERT measurements were carried out using the IRIS Syscal Pro Switch 72 resistivity meter. At Site A (Figures 8), seven profiles were acquired, each of resistivity data mapped by Wenner-Shlumberger and Pole-Dipole arrays with the parameters given in Figure 8 caption. Sections of the effective resistivity distribution at depth were generated from the inversion of the apparent resistivity measurements using the TomoLab code (http://www.geostudiastier.it/). The smoothed least squares inversion algorithm uses quasi-Newton optimization to minimize the misfit between the model's response and the observed data values (Loke & Barker, 1996;Loke & Dahlin, 2002). At Site B (Figure 9), three 2D ERT profiles were mapped using the Wenner-Schlumberger and pole-dipole arrays. The NE-SW oriented profile 18 runs a few hundred meters eastwards of the La Regina and Colonnelle Lakes. This profile images variation of the travertine thickness and the resistivities of the underlying lithological units to about 130 m BGL (Figures 5d and 9a). The parallel E-W oriented profiles 19 and 20 extend a few meters to the north and south of the two lakes (Figure 5d and Figures 9b and 9c).

ERT 3D
At Site A, ERT 3D surveys were carried out over blocks B1 and B2 (Figure 5b and Figure 10a) to the northeast and southeast of the subsiding area. The two blocks image areas where abrupt variations in resistivity are evident in the ERT 2D results. In both blocks, the measurements were carried out using the pole-dipole (pd) array. The cable layout was optimized to acquire true 3D data both along and across the cables. Therefore, a pseudo-volume or cloud of apparent resistivity values was generated and related to 3D resistivities estimates of a finite element model of the subsurface by least square inversion with smoothing weights using the ErthLab code (http://www.geostudiastier.it/). At Site B (Figures 5d and 11), the southwestern ERT 3D block measured the circular depression mapped in the LiDAR-derived DTM (F1 in Figure 10a) that aligns between the two lakes. The 3D ERT resistivity model was obtained by interpolating four 2D parallel ERT profiles. The profiles used a WS array with 48 electrodes spaced 5 m apart, whereas the profile spacing was approximately 10 m. The 3D block strikes WNW-ESE and is approximately 240 m long by 40 m wide and 42 m deep.

InSAR
To enhance the details in the InSAR results, the PS values were gridded with a cell dimension of approximately 1,000 × 1,000 m. The LOS velocity and accumulated near-vertical deformation maps are shown in Figures 3 and 4. The investigated area (12 × 8 km 2 ) encompasses a large portion of the AAB's northern sector. Both sites A and B are mostly located within agricultural lands with a limited number of PS, whereas the neighboring villages, by contrast, are characterized by more densely clustered PS values.  Figures 3a, 3b, and 3c. Additionally, the 2012-2016 period's CSK data set in Figure 3c shows the

Earth and Space Science
deformation velocity of about −12 mm/year for the subsiding SSA and CSSA regions. During this observation period, the peak subsidence velocities were along the highway in Site A and the railway in Site B. Figure 4a indicates the coverages of the three InSAR data sets within the SSA and CSSA, whereas the ground velocity profiles 1 and 2 in Figures 3d and 3e show the deformation pattern across the SSA.
The cumulated CSK ground deformation map (Figure 4b) shows higher subsidences in the southern SSA about the lakes, as well as between the Bagni di Tivoli and Villalba villages, and along the railway. For example, the subsidences reach −45 mm some 700 m southeastwards of Regina Lake, −107 mm approximately 480 m southeastwards of Regina Lake, and −200 mm in the Villalba village. In addition, the maximum subsidence values near to and within the A and B study sites generally trend NS and NE-SW parallel to the strikes of the regional structures and the SSA.

LiDAR
According to the LiDAR-derived raster elevation model (DEM), the north side of the AAB lies in three main sectors (Figure 5a). The first sector includes the highest average elevations (75 < z < 100 m) that are approximately aligned NE-SW south of Guidonia and then arc NS along the northern edge of the basin. The second sector involves transition values (70 < z < 75 m) between the first sector and the flatter and more subdued Analysis of the DEM at Site A reveals the elliptical D1 depression elongated in the NNE-SSW direction with the circular D2-depression some 200 m northwestward on strike with the D1's long axis and the circular D3 depression that merges with D1's western margin (Figures 5a, 5b, and 5c). However, further study of the D2 depression was not possible because it is in the military area with restricted access. The DEM for the Site B reveals the circular 250 m diameter F1 topographic depression that is less than 200 m east of the La Regina lake (Figure 5d).

Site A 4.3.1. Boreholes
The stratigraphic sequence in the S1 borehole involves alternating lithotypes over the depth range of 0-23 m BGL (Figure 2b). The sequence consists of compact travertine benches with thicknesses greater than 3 m and loose clay-to-sand and gravel layers. A massive travertine was found from 23-to-30 m BGL at the bottom of the borehole. During the geophysical campaign, the S1 piezometer showed a water level oscillation of approximately 14 m (Figure 2f).
The S2 borehole sampled the lithologies lying below the subsiding area reaching the −60 m BGL depth. The stratigraphic sequence consisted of soil and gravel on a silty matrix in the upper three meters (Figure 2b), fractured travertine vacuolar-to-massive benches down to 8 m BGL, and a homogeneous sequence of clay, dark varved silty clays, and peat to the borehole's bottom. The S2 terrigenous lithotypes have horizons with high-water content ranging from 37% to 54%, whereas the degree of saturation varies from 78% to 100% as shown in the supplemental Table S1.

Gravimetry
The residual Bouguer amplitude values (Figures 6a and 6b) range from −0.29 to 0.16 mGal with the positive values slightly lower from those measured by Argentieri et al. (2015). This difference most likely reflects the use of different gravimeters (intrinsic error) and different absolute gravity stations of the International Gravity Standardization Net (IGSN71). The map shows the well-defined L 1 anomaly minimum or low covering a substantial portion of the surveyed area map. It represents the southward extension of the large minimum Argentieri et al. (2015) detected over the subsiding area's D1 depression. The positive H 1 and H 2 anomalies rim the L 1 minimum, which is not completely sampled by the survey. Additionally, the L 2 minimum is also incompletely sampled in the southeastern border of the survey. Figure 6b highlights NE-SW elongation of the L 1 gravimetric minimum, which contrasts with the more circular minimum mapped by Argentieri et al. (2015).

S2 Core Magnetic Properties
The detected magnetic anomaly intensity (Figure 6c) ranges from −150 nT to 320 nT, implying the existence of magnetized units in the study area, particularly inside the buried depression. Analysis of the magnetic properties of the S2 borehole samples was carried out to investigate the origin and characteristics of the magnetic sources derived mainly by the distribution of allogenic magnetic minerals in the terrigenous sequence (e.g., the pyroclastic outcrops surrounding the AAB) and authigenic mineralization in the clays.
The magnetic volume susceptibility (κ), the natural remanent magnetization (NRM), and the anhysteretic remanent magnetization (ARM) intensities are commonly used in environmental magnetism as magnetic concentration parameters because these provide information on the concentration of magnetic minerals in the sediment (Peters & Dekkers, 2003;Liu et al., 2012). Other parameters are also commonly used like (1) ARM/κ ratio that provide information on grain size because ARM is more sensitive to changes in the stable single domain (SD) components (≈0.1-0.03 μm) (King et al., 1982;Maher, 1988), while κ is more sensitive to changes in the coarser, multi domain (MD) components, besides SP ferrimagnetic particles (Dearing et al., 1996); (2) the median destructive field of the NRM (MDFNRM), defined as the AF value necessary to reduce the NRM to half of its initial intensity, indicates a variable resistance to the AF demagnetization treatment (computed using the DAIE software Sagnotti, 2013); and (3) the ARM30/ARM0 ratio that provide information on magnetic coercivity (Dankers, 1981;Venuti et al., 2005).
In this study, the magnetic volume susceptibility (κ) and the natural remanent magnetization (NRM) and anhysteretic remanent magnetization (ARM) intensities were determined from samples collected from the S2 core to help in constraining the possible near-surface sources of the surveyed magnetic anomalies. Figure 7 displays enhanced magnetization at the top of the S2 core and in the interval between 8 and 37.5 m of silty clay and clay. These parameters have similar trends, indicating an overall broadly ferromagnetic contribution, although small differences between κ and the remanent magnetizations exist and suggest the presence of diamagnetic and/or paramagnetic minerals in some thin intervals. The low MDFNRM (<30 mT) suggests that a low coercivity mineral like magnetite or maghemite may be the main magnetic carrier. Differences between ARM and NRM intensity reflect magnetic domain size differences, where ARM is enhanced in smaller sized (single domain [SD] and pseudo-single domain [PSD]) magnetite. In the up-core higher magnetization interval, the ARM/κ and ARM30/ARM0 logs suggest a slight increase in the magnetic grain-size and decreases in coercivity, which likely correspond to an increase in the MD magnetite grains with respect to the SD magnetite in the sediment. The few samples at 0. 43, 5.62, 8.05, 43.73, 50.66, and 59.32 m acquire magnetic remanent components during NRM demagnetization. This finding is typical of Gyro Remanent Magnetization (GRM) and indicates the presence of sulfides such as greigite (Snowball, 1997;Stephenson & Snowball, 2001). These results constrained the inverse modeled source depth, values, and magnetic mineralization estimates.

Magnetic Anomaly Fields
The complex TMI anomalies highlight the distribution of magnetic sources in the subsiding area beneath the D1 depression (Figures 6c and 6d). The anomalies exhibit predominantly induced N-S dipolar attributes. Within the D1 depression anomaly wavelengths range from a few meters to 70-80 m, and they and their strikes also follow the D1 depression's edges and fracturing. The prominent anomalies at the eastern border are due to the effects of an iron fence. They truncate the E-W striking A1 anomaly and the N10E trending A2 anomalies on the eastern side of the subsiding sector (Figure 6c). The NW-SE trending A3 and A4 anomalies lie in the central part of the D1 depression and merge with the A1 and A2 anomalies. The N10E trending A5 anomaly parallels the A2 anomaly and intersects the concentric A6 and A7 anomalies. The A8 anomaly forms an arcuate band along the southeastern border of the D1 depression. The RTP-TMI anomaly transformation combines the A1 and A2 TMI anomalies into the RTP-A1 anomaly along the eastern side of the subsiding area (Figure 6d). The RTP process also has revealed a well-defined circular anomaly on the northern side of the RTP-A1 anomaly. The RTP-A2 and RTP-A3 features are a few meters north of their TMI anomalies. In addition, the linear attributes of the RTP-A5 anomalies are better defined, and the circular geometries of the RTP-A6 and RTP-A7 anomalies are enhanced.

Magnetic Modeling
The obtained magnetization model of the subsurface (Figure 12) was evaluated for magnetic source depths and intensities in the context of the S2 borehole samples (Figures 7). Using an applied field intensity of F = 46,374 nT that operated during the surveying (https://www.ngdc.noaa.gov/geomag/calculators), as well as the measured (on the samples) maximum susceptibility k max = 0.0004 SIu and the NRM of 0.06 A/m, obtains a total magnetization estimation of approximately 0.0748 A/m. The inverse modeling of the magnetic anomaly, on the other hand, obtained susceptibility estimates ranging over 0.0011 to 0.0035 SIu corresponding to an estimated magnetization range of 0.040 to 0.1 A/m, assuming parallel induced and remanent magnetizations (Paine et al., 2001). Consequently, the measured and modeled magnetic properties are generally consistent, where the misfit between the observed and predicted anomalies ranges between −5 and 12 nT (Figure 6e) or roughly 8% of the amplitude range of the RTP-TMI anomalies. The inversion results in Figure 12 include (a) the area of the magnetic data inversion overlaid on the Site A's DEM, (b) estimated susceptibilities over the x-y-planes at depths of 10 and 30 m BGL, (c) an isometric plot of the susceptibility envelop for k = 0.06 SIu, and (d) the shallowest depths of the 0.06 SIu susceptibility. Modeling suggests that the shallower susceptibility map in Figure 12b highlights the possible distribution of magnetic sources over depths ranging from 0 to 10 m BGL with geometrical attributes that correlate with the topography in and around the subsiding area of the D1 depression. The magnetization model also suggests that the sources for the region from 10 to 40 m BGL follow the deeper susceptibility distribution in Figure 12b and coalesce into a concentrated positively magnetized feature in the eastern portion of the D1 depression. The magnetization model did not require sources deeper than 40 m BGL to accommodate the observed magnetic anomalies.
The isometric plot of the magnetic susceptibility envelop for k = 0.06 SIu in Figure 12c has a thickness of approximately 20 m and extends along a SE-dipping plane. Figure 12d shows the depth to top surface of the isometric plot in panel (c) and infers a N20°E trending channel that flows from the shallower t2 trough into the magnetic t1 depression, where two small sinkholes also appeared in the 2013 (Argentieri et al., 2015). In addition, the magnetic depression t3 underlines the D3 topographic depression highlighted in the DTM of Figure 5c.

2D ERT
Interpretation of the ERT sections at Site A ( Figure 8) is constrained by the borehole stratigraphy. The effective electrical resistivity may be broadly classified into the blue-colored unit of the relatively conductive (3-100 ohm-m) continental-marine terrigenous sequence (e.g., blue-clays, silt, silty clays, clayey silt, and peaty clays) and the red-colored unit of the relatively resistive (100-10,000 ohm-m) travertine.
The A, C, and AC features are highly resistive (up to 3,500 ohm-m) with an average thickness of approximately 45 m. These travertine features behave variably along the edges of the subsiding area. For example, in profile P4, the resistive (up to 200 ohm-m) t3 feature extends from approximately 2 to 25 m BGL (Figure 8b). Beneath t3, the A unit forms a resistive wedge with an eastwards dipping upper surface. In profile P6, the B1 conductive channel that separates the resistive AC and A1 features is approximately 30 m deep. Here, the A1 feature forms a near vertical dipping block with the vertical f discontinuity (Figure 8c) between the A1 and B features. In this tomography, the thickness of the AC feature gradually decreases NNW from roughly 60 to 40 m towards the subsiding area. In general, the A feature has an undulating upper surface that is inclined to the NW by about 10°-15°in profiles P6-14 (Figures 8c and 8d) or forms a stair-step structure as in the western edge of the subsiding sector (profile P4 in Figure 8b). The topography reflects the travertine development of stairway-like slopes, fractures, and sinkholes along the edges of the D1 depression (i.e., the d1, d2, and d3 step-like depression in profile of the P.2 and the d2, d3, and d4 depressions in profiles P.6 of Figures 8a and 8c).
The conductive B feature of resistivities <80 ohm-m below the subsiding area (profiles P2-4-6-14 in Figures 8a-8d) has an average thickness of approximately 60 m. This terrigenous feature is constrained by the S2 borehole and exhibits quite variable geometric attributes in the sections.
In profiles P15 and P16 at the southern edge of the subsiding area, the A, B, and C features are differentially distributed. The A and C features maintain average thicknesses of approximately 60 and 100 m, respectively. However, the resistive AC feature appears in the central upper part of the sections, interrupting the electrical frameworks in the other sections (Figures 8e and 8f). The almost flat AC feature located at approximately 10 to 35 m BGL suggests thin travertine, which thickens southwards from approximately 20 to 65 m. The thickening of the travertine southwards (hat is also visible in the profile P6 in Figure 8c possibly suggests different sedimentary condition by this side of the D1-depression. The conductive D feature is approximately 5-10 m thick and lies in the upper portions of all the tomographies. This feature is poorly resolved by the ERT because its thickness is typically smaller than the electrode spacing. However, the upper portion of feature D includes the agricultural soils (3-4 m thick), whereas the underlying electrical layers contain the sand, gravel, and clays sampled by the S1 and S2 boreholes (Figure 2b). Figures 10a-10c show the three-dimensional (3D) electrical resistivity patterns for the D1 depression at its northeastern and southeastern edges. The travertine A feature with resistivities of 1 to 5,000 ohm-m has an arcuate boundary with a corrugated 100 ohm-m surface envelop (Figure 10b) approximately 2 m BGL that exhibits grooves, fractures, and disjoined blocks separated by low resistivity channels. The area below the subsidence in the central part of Figure 10c has very low resistivity (<10 ohm-m). Figures 10d and 10f show the three-dimensional resistivity tomography for the nonsubsiding region southeast of the D1 depression that includes the D4 depression (Figure 5b). The surface of travertine A feature with 220 to 500 ohm-m resistivities deepens from 5 m BGL near its arcuate edge to roughly 7 m BGL at the center that lies approximately 5 m BGL and deepens at the feature's center to about 6.5 m BGL (Figures 10d  and 10e). This geometry is also reflected in the topographic profile of Figure 10e, which shows a ridge (R) and a depression (D) sequence in the NW-SE direction related probably to karsting in the near-surface travertine. The resistive B feature of 10 to 190 ohm-m resistivities lies in the 3D block's northwestern corner facing the cavity.

Site B
At Site B, the subsidence analysis used both the 2D and 3D electrical tomographies as constrained by the SG-5 and SG-10 borehole stratigraphy data (Figures 2a and 2d) on the presence of terrigenous, travertine, and volcanic lithological units (Annunziatellis et al., 2009). The SG-10 borehole intercepted 40 m of travertine to the bottom of the borehole, whereas the SG-5 borehole sampled 10 m of travertine and approximately 50 m of volcanic rocks. Based on these results, the effective resistivity variations at Site B may be broadly classified as conductive, low-resistivity 3 to 100 ohm-m features related to continental-marine lithologies involving blue-clays, silt, silty clays, clayey silt, and peaty clays; low conductivity, high resistivity 100 to10.000 ohm-m features related to the travertine; and moderately conductive and resistive 50 to 500 ohm-m features related to loose pyroclastic and ash deposits and other volcanic lithologies. The travertine and pyroclastic lithologies have overlapping resistivity values that may be partially differentiated based on the borehole depths of these lithologies.

ERT 2D
For the region between the Regina-Colonnelle and S. Giovanni Lakes, Figure 9a shows the electrical resistivity in profile P18, which involves the thick resistive H 1 -18 feature with considerable bottom depth variation ranging from 26 to 82 m and the thin resistive H 2 -18 feature with an approximate thickness of 25 m and an almost flat lower surface. The two resistive features are separated along the red-dashed vertical boundary labeled f in Figure 9a. The upper and lower portions of these two features respectively may consist largely of travertine and pyroclastic lithologies. Below these shallower resistive units, the conductive L 1 -18 feature is evident throughout the profile's length and encloses the resistive H 3 -18 feature with resistivities up to 213 ohm-m at the profile's NE end. The conductive L 1 -18 feature may be related to the terrigenous continental-marine sequence. The altimetry data in Figure 9a show possibly collapsed cavities that include the t4 feature at the mean altitude of 69 m asl, and the t1, t2, and t3 inferred troughs at the average altitude of 67.5 asl. The inferred t4 depression coincides with the F1 feature observed in Figure 5d, and beneath F1 in the resistive upper surface is the inferred d3 cavity (Figure 9a).
The tomography profile P19 in Figure 9b shows that the conductive low resistivity L 1-2 -19 features with resistivities of 7 to 35 ohm-m dominate the section's eastern three quarters. These features may reflect fractured travertine fed with hydrothermally circulated groundwater. At the eastern side of the section between the lakes, the resistive H 1 -19, H 2 -19, and H 3 -19 features with resistivities of 100 to 500 ohm-m occur down to depths of roughly 52 m BGL that may represent massive anhydrous travertine. The altimetry profile of Figure 9b bifurcates into the western side d1-level near Colonnelle Lake and the slightly lower d2-level on the eastern side.
The tomography profile P20 in Figure 9c shows the electrical characteristics of the area boarding the lakes on the south, where the wide resistive H-20 feature occurs with resistivities up to 1,900 ohm-m. This pervasive feature extends some 50 m from the surface to the bottom of the section. It is disrupted by the conductive L 2 -20, L 1 -20, and L 3 -20 features with resistivities up to roughly 25 ohm-m developed at the eastern and western side of the sections, respectively. The vertically dipping L 1 -20 unit (10 ohm-m) lies a few meters south of Regina Lake. Relative to the resistive H-20 feature that may reflect compact travertine the conductive L 1 -20, L 2 -20, and L 3 -20 features may involve fractured travertine in contact with thermally circulated groundwater.

ERT 3D
The three-dimensional (3D) ERT block covers the eastern margin of the F1 depression (Figures 5d and 11). Here, the volume is mostly occupied by the resistive H 1 -3D and H 2 -3D blocks that may represent massive travertine. The H 1 -3D block with resistivities up to 5,700 ohm-m has a bowl shape that corresponds with the eastern flank of the F1 feature (Figure 11c). It is about 40 m thick and extends from a few meters BGL to the bottom of the tomography. The low resistivity L 2 -24 feature (Figure 11b) is approximately 7 m thick and may represent the thin sedimentary fill of the F1 depression. The H 2 -3D block unit on the section's eastern end with resistivities up to 2,000 ohm-m and minor depth extent of about 20 m may represent the thinning of the travertine plateau perhaps via the d1 channel mapped as the H 1 -18 feature in the ERT 2D profile P18 in Figure 9a. By this interpretation, the low resistivity L 1 -3D feature (L 1 -24 feature in Figure 11b) may represent partially saturated karst debris fill in the D1 depression.

Discussion
The InSAR and LiDAR data for the northern AAB that includes the A and B sites show extensive subsidence over the northern AAB and additional buried paleo-morphologies such as Site A's D1, D2, D3, and D4 depressions and Site B's F1 cavity ( Figure 5).
The subsidence of the SSA over the 1990-2016 period evolved in the NS and NE-SW directions (Figures 3a-3c and 4a and 4b) along regional tectonic trends (Faccenna et al., 2008) to a maximum of about 43 cm near the SSA's center. The subsided area corresponds to the gravimetric minimum that Di Nezza and Di Filippo (2015) interpreted for a local deepening of Mesozoic carbonate bedrock. Furthermore, SSA subsidence corresponds to mapped subsurface patterns of the electrical conductivity, temperature, and water table in the travertine (Nanni et al., 2014). Accordingly, subsidence in the northern AAB may coincide with the movement of trans-tensional basin segment created by the Sabina Fault (Faccenna et al., 2008), even though the AAB appears to lack outcrops with clear evidence of tectonic movement. Alternatively, the observed subsidence in the SSA may be caused by karst dissolution in the deeper pre-travertine Mesozoic bedrock that rearranged the overlaying sedimentary sequence (Coltorti, personal communication).

Site A
At Site A, subsidence is restricted to the 200 × 180 m 2 area of the D1 depression (Figures 5b and 5c). The borehole stratigraphy and the respective 2 and 3D ERT resistivity models in Figures 2b, 8, and 10 relate the D1 depression to the subsiding terrigenous sequence filling a cavity in the travertine. Indeed, the resistivity contrast across the margins of the conductive subsiding area into the resistive region is quite apparent. However, along the margins of the cavity, the travertine thins to a thickness of about 40 m in disarticulated blocks with marked top surface variations that appear to cause irregularities on the ground's surface topography. For example, in profile P2, the top surface of the electro-layers (t features) correlates directly with the topographic surfaces (d-features); see the t1-d3 and t2-d4 features in Figure 8a; in profile P4 see the t1-d3, B1-d1, and t2-d2 features in Figure 8b; and in profile P6 the t1-d1, t3-d3, and t2-d2-d4 features in Figure 8c.
On the D1 depression's eastern and northern sides, the ERT 2D (Figures 8a-8d) suggest a vertical contact between the travertine (A feature) and the terrigenous sequence (B feature). The ERT 3D also shows disarticulated travertine blocks located roughly 5 m BGL in the B1 block (Figures 10b and 10c). The corresponding ground topography features 1-2 m undulations, small sinkholes, and fractures (Argentieri et al., 2015).
The electrostratigraphy of the western side suggests the travertine C feature with a stair step-like surface towards the center of the D1 depression (e.g., t1 in Figures 8a and 8c). On this side, fractures and sinkholes seem to be lacking, with the notably lower subsidence yielding the 3 m diameter doline-like D3 depression (Figure 5c).
On the southern side of the D1 depression, the electrostratigraphy is completely different. Here, the travertine surrounding the D1 depression gradually increases in thickness southwards and westwards up to 60 and 98 m, respectively (Figures 8c, 8e, and 8f). Near the D1 depression, travertine may occur the thinner, more conductive AC feature at shallower depths of about 30 m. These results also suggest that the travertine gradually pinches out between the ERT profiles P14 and P16. The B feature lying beneath the AC feature is significantly less conductive than the B unit within the D1 depression (Figures 8e and 8f). The travertine thickness variation of the D1 depression's southern edge suggests different former sedimentary conditions 10.1029/2019EA000870 Earth and Space Science or alternatively, an erosive process which might have affected the travertine. The pinch out geometry of the travertine deposition suggests that the first hypothesis is more realistic.
In addition at the D1 depression's southeastern side, the topography forms a bulge (Figure 5c) as detailed by the stair-step like surfaces of the travertine in Block 2 of the ERT 3D data (Figure 10e). The bulge probably represents part of an approximately 7 m deep doline located a few meters outside the area of the subsiding D1 depression. This result together with the D3 depression suggests that karst cavities may be coalescing into a much larger depression at Site A. More generally, these results suggest that the topography correlates strongly with karsting of the shallower travertine.
The new gravimetric survey shows that the residual L 1 minimum overlies the D1 depression's area beyond the limits identified by Argentieri et al. (2015). In fact, the L 1 anomaly is elongated in the NE-SW direction (Figure 6b) in contrast to the circular geometry observed by the earlier study. The L 1 low represents the superimposed gravity effects of both the travertine and terrigenous lithologies that underline the D 1 depression (e.g., gravimetric profile above the 2D ERT profile P4 in Figure 8b). The L 1 minimum is related to the negative density contrast between the compact higher density travertine surrounding the D 1 cavity and the lower density terrigenous sequence and the fractured and karstified travertine beneath and along the D 1 margins. The residual L 2 gravity low that was not completely sampled may represent the signature of the D4 doline depression observed in the 3D ERT data of Figures 10e and 10f.
The magnetization model in Figure 12 constrains the subsurface distribution of the magnetic sources underlying the D1 depression. In addition to the sources depth and susceptibility estimates, the borehole sampled magnetic minerals in the 8 to 37 m depth range (Figures 7 and 12b-12d) help to clarify the magnetic lithology of the subsurface. Consequently, the magnetic sources appear to be caused mostly by the higher concentration of magnetic minerals in the terrigenous sequence with mean susceptibility κ mean = 0.00022 SIu) along mainly the eastern border of the D1 depression (Figures 6c and 6d). This concentration was likely deposited by circulating groundwater in the deep fractures.
Finally, the coincidence of the t3 depression in Figure 12d with the D3 doline in Figure 5c suggests that the magnetization processes along the eastern side of the D1 depression may also have been active within the smaller D3 doline. Therefore, the magnetic mineralization appears to be related to the water circulating through the karst formations.
The magnetic property measurements from the S2 borehole in Figure 7 show that the main magnetic carrier is low coercivity magnetite/maghemite with some "ferromagnetic" Fe-sulfides mixed. During diagenesis, small amounts of authigenic magnetic minerals may be altered both chemically and physically. During early diagenesis, iron-reducing bacteria often contribute to the dissolution of Fe-bearing mineral (Liu et al., 2012;Roberts, 2015). During the travertine deposition, erosion of the surrounding volcanic hills and the local hydrological network mobilized significant quantities of magnetic minerals. When the circulating groundwaters reached the travertine deposits, they flowed through the karstic cavities depositing the coarser mineral fractions. The S2 core, for example, contains ferrous mineral grains with millimeter to centimeter dimensions as shown in the supplemental Figure S2. Therefore, the differences in magnetic particles sizes may be related to episodic cavity deposition of the sediment eroded from the study area's volcanic terrains. Figures 13a and 13b outline the inferred attributes of the travertine D1 depression, which include the vertically dipping northern and eastern boundaries, the stair-step western boundary, and the pinched out southern boundary. Intense karstification is common in the area, where the D1 depression's fill consists of some 40 m of highly compressed marshy-lacustrine deposits (Figure 2b), and magnetic minerals concentrated along an inclined channel at the eastern boundary (Figures 13b and 12d). These results reveal an ancient buried morphology for the D1 depression that lies a few meters BGL. In addition, the D1 depression lies in the northeastern AAB, which is the main geothermal resurgence area for the basin (De Filippis et al., 2013).
In general, the D1 depression's paleo-morphology may represent the remnant of an ancient hydrothermal source developed along a zone of weakness like a NE-SW fault of the carbonate basement (Klimchouk et al., 2017) that generated a lake and surrounding marsh environment. Travertine encrustations in borehole S2 between the depths of 59.20 and 63.90 m supports this hypothesis. With the completion of resurgence, karst and hyperkarstic processes developed filling the D1 depression with terrigenous sediments. In

10.1029/2019EA000870
Earth and Space Science general, sin-sedimentary sinkholes (Gutiérrez et al., 2008) formed on paleo-topographic surfaces beneath the Holocene cover are common in the AAB Funiciello et al., 2005 andKaufmann, 2014). Thermal springs aligned in the NE-SW direction near the study area are also observed (De Filippis et al., 2013).
The subsidence within the D1 cavity may reflect the high piezometric drop of the travertine aquifer invoked by the local paroxysmal subsidence. The stratigraphy from the S2 borehole shows a sequence of units with high-water content (see supplemental Table S1) that corresponds to the conductive B feature in Figure 8. Ground displacement histories in the AAB recorded by mining operations highlight widespread soil compression and dewatering events (Bozzano et al., 2015). In fact, oscillations of groundwater table in the past few years were about 10 m and up to approximately 14 m in the S1 piezometer measurements of 2016 ( Figure 2f). Thus, when the water is pumped, a strong pressure consolidates the compressible soils within the cavity. Moreover, the confinement conditions within the 200 × 180 × 40 m 3 D1 cavity further amplify compressing and consolidating the soils.

Site B
The electrical tomography of Site B reconstructs the paleo geographical landscape as constrained by the data from boreholes S5 and S10 and the geological model proposed by Annunziatellis et al. (2009) and summarized in Figure 2e. The ERT profile P18 (Figure 9a) shows the resistive H 1 -18 feature with six undulations (A1-6) along its interface with the more conductive L 1 -18 feature. The upper part of the H 1 -18 feature involves travertine with a westward decreasing thickness as interpreted in Figure 13c, whereas the underlaying L 1 -18 feature is composed of volcanic lithologies as suggested by the SG5 borehole data in Figure 2d.

Earth and Space Science
Thesix undulations at the bottom of the H 1 -18 feature also may provide electrical resistivity evidence for buried narrow valleys filled with Alban pyroclastics.
The resistive H 2 -18 layer, by contrast, represents travertine of relatively uniform thickness and horizontal interface with the underlying more conductive L 1 -18 feature shown in Figures 9a and 13c. Here, the thermal resurgences of the two lakes affected the travertine and spread a large number of electrolytes throughout the aquifer (Figure 13d).
In Figure 9a, the resistive H 1 -18 and H 2 -18 features are separated by the vertical f-discontinuity that may provide evidence of the presumed faulting across the AAB at Site B as interpreted in Figures 9a and 13c. Further evidence for the faulting has been inferred from the discovery of radon, thorium, and helium gas anomalies aligned along the NS direction (Chiodini et al., 2004;Faccenna et al., 2008;Pola et al., 2014;Soligo et al., 2002). The interpretation of the section beneath 60 m BGL lacks borehole constraints. However, the widespread conductive L 1 -18 feature probably represents the distribution of the continental and marine terrigenous sequences, although the boundaries between the sequences are not electrically discernible. The resistive circular H 3 -18 feature with resistivities up to 215 ohm-m likely indicates the lenticular disruption of in the terrigenous sequence probably by conglomerate deposits.
In addition, the 3D ERT profile P24 in Figure 11a samples roughly 35 m of the subsurface across the F1 depression imaged by the high-resolution DTM of Figure 5d. These results show that the depression's topographic rim coincides with an upward protrusion of the resistive travertine H 1 -24 and H 1 -3D features in Figures 11b and 11c, respectively). Thus, the F1 depression may represent karst effects approximately in line with the EW trend of the Colonnelle and Regina Lakes. Moreover, the low resistivity d1 vertical channel (Figure 9a), which corresponds to the conductive L 3 -20 and L 1 -24 features of the respective ERT 2D profiles P18 and P24 in Figures 9c and 11b, and to the conductive L 1 -3D feature in the ERT 3D block of Figures 11c  and 11d, suggests that the travertine thickness was strongly affected by hyperkarst and karst phenomena, most likely from hydrothermal ground water circulation. Figure 13c summarizes the modeling of Site B in terms of a near-surface sequence of narrow valleys in the continental terrigenous units filled with pyroclastics and covered by an almost flat travertine layer. The upper sector represents the pre-eruptive badland landscape high-gradient clay terrains (e.g., the "Calanchi" landforms of the Lazio region) surrounding the Colli Albani volcano (Howard, 1994, Del Monte, 2017, whereas the lower sector depicts the almost flat landscape of the AAB.
At Sites A and B, the topography correlates strongly with both shallow and deep travertine structures. However, the shallower travertine processes clearly have the greatest influence on the topographic and related environmental hazards.

Conclusions
The synthesis of the remote sensing, gravity, magnetic, and electrical resistivity tomography data mapped vertical ground movements for a wide area of the AAB about Sites A and B and allowed the identification of buried paleo-forms located a few meters below the ground surface. The results show a direct relationship between these paleo-forms and the surface topography and its modifications of the above topography. Furthermore, the study sites were found to be located in the area of the AAB with the highest subsidence rate over the past 20 years, where pre-travertine bedrock helped to drive the travertine plateau's evolution at both regional and local scales. Karst and hyperkarst processes drove and continue to drive the production and infilling of sinkholes and other buried paleo-forms that likely extend beneath the AAB's urban areas.
The above results describe an active geomorphological system, where several factors affect geological risk. For example, the travertine a few meters BGL drives changes in the surficial topography even small scales. From the beginning of travertine deposition, karst and hyperkarst processes have been producing, buried springs, sinkholes, and vertical dissolution chimneys at depth that may not be perceptible at the topographic surface. These elements critically affect land use planning for AAB, especially where urbanization may alter the geomorphological and hydrological dynamics of the travertine plateau and thus trigger accelerated development of subsurface cavities and regional and local subsidence activities.
This investigation demonstrated a novel approach for studying subsidence and sinkhole phenomena by integrating data from remote sensing and geophysical surveying. The LiDAR and InSAR observations helped to identify ground movements and characterize small-scale morphologies to help constrain in combination with the borehole data, the subsurface geophysical interpretations of the gravity, magnetic, and electrical resistivity survey data. More generally, the proposed workflow is very effective for mapping and monitoring the geohazards attributes of karst-prone regions like the AAB.