Imaging Multi-Age Construction Settlement Behaviour by Advanced SAR Interferometry

: This paper focuses on the application of Advanced Satellite Synthetic Aperture Radar Interferometry (A-DInSAR) to subsidence-related issues, with particular reference to ground settlements due to external loads. Beyond the stratigraphic setting and the geotechnical properties of the subsoil, other relevant boundary conditions strongly influence the reliability of remotely sensed data for quantitative analyses and risk mitigation purposes. Because most of the Persistent Scatterer Interferometry (PSI) measurement points (Persistent Scatterers, PSs) lie on structures and infrastructures, the foundation type and the age of a construction are key factors for a proper interpretation of the time series of ground displacements. To exemplify a methodological approach to evaluate these issues, this paper refers to an analysis carried out in the coastal/deltaic plain west of Rome (Rome and Fiumicino municipalities) affected by subsidence and related damages to structures. This region is characterized by a complex geological setting (alternation of recent deposits with low and high compressibilities) and has been subjected to different urbanisation phases starting in the late 1800s, with a strong acceleration in the last few decades. The results of A-DInSAR analyses conducted from 1992 to 2015 have been interpreted in light of high-resolution geological/geotechnical models, the age of the construction, and the types of foundations of the buildings on which the PSs are located. Collection, interpretation, and processing of geo-thematic data were fundamental to obtain high-resolution models; change detection analyses of the land cover allowed us to classify structures/infrastructures in terms of the construction period. Additional information was collected to define the types of foundations, i.e., shallow versus deep foundations. As a result, we found that only by filtering and partitioning the A-DInSAR datasets on the basis of the above-mentioned boundary conditions can the related time series be considered a proxy of the consolidation process governing the subsidence related to external loads as confirmed by a comparison with results from a physically based back analysis based on Terzaghi’s theory. Therefore, if properly managed, the A-DInSAR data represents a powerful tool for capturing the evolutionary stage of the process for a single building and has potential for forecasting the behaviour of the terrain – foundation – structure combination.


Introduction
In the frame of geotechnical risks, subsidence is a relevant threat to human activities and urban settlements [1]. Structures and infrastructures can be considered at the same time a triggering/accelerating factor of subsidence and the elements at risk. Among the main human activities causing or accelerating subsidence are the exploitation of underground resources (e.g., water and hydrocarbon withdrawal), shallow water table fluctuations, and external loading caused by urbanization [2][3][4][5][6][7][8][9][10][11].
Subsidence processes often affect flat regions, such as alluvial and coastal plains, particularly suitable for urban expansion. It is worth noting that, from a geological point of view, these environments usually show a very complex geometry of stratigraphic successions in the subsoil, featuring lateral and vertical variations in facies even on a short wavelength.
The main tangible effect of subsidence is the differential settlement and the related impact on structures and infrastructures that can eventually lead to the exceedance of the serviceability limit up to the structural damage. Due to these consequences, proper actions should be undertaken to mitigate subsidence risk, starting with the detection and monitoring of the process.
Currently, several tools and techniques are available for measuring and monitoring the deformational patterns of subsidence with very high accuracy and precision. Among these, in the last 20 years, a revolution has been initiated by the use of satellite Synthetic Aperture Radar (SAR) images. SAR interferometry is particularly suitable for detecting subsidence processes, as they primarily feature vertical deformation that can be well-resolved by combining results from ascending and descending satellite orbits. In addition, the availability of SAR imagery starting from 1992 allows for the recovery of long time series of the deformation pattern, which implies the possibility to perform back-analyses. For this purpose, the most reliable interferometric approach is advanced SAR interferometry (A-DInSAR) based on the Persistent Scatterer Interferometry (PSI) technique [12]. A-DInSAR techniques can be profitably exploited to investigate the spatial and temporal evolution of subsidence processes with millimetric accuracy [2,[13][14][15]. A-DInSAR methods are characterized by analyses of multi-temporal data stacks that allow us to generate several interferograms, thus achieving redundant and time-related Interferometric patterns [12,[16][17][18][19]. PSI is based on the information obtained from pixels of the SAR images characterized by high coherence over long time intervals [12,16,17]. Generally, man-made structures, such as buildings, bridges, dams, railways, or pylons, or natural elements, such as outcropping rocks or homogeneous terrain regions, can represent good persistent scatterers (PSs). As stated by Hanssen [20], A-DInSAR is an "opportunistic" technique, i.e., the presence of potential good scatterers is fundamental to have to achieve a sufficient number of PSs over the investigated region.
InSAR-based techniques can contribute to the analysis of the cause-and-effect relationships between subsidence and predisposing/triggering factors [2,5,8,9,21,22]. Indeed, the spatial pattern of deformation rates is commonly related to the spatial distribution of predisposing (e.g., the thickness of compressible deposits in the subsoil) and/or triggering (e.g., the location of pumping wells or imposed loads) factors. Furthermore, the time series of ground/structural displacements can provide useful hints for interpreting the temporal pattern of the process and, thus, the relationship with boundary conditions: the behaviour over time of the displacements can differ significantly if triggered by, for instance, groundwater withdrawal, which is potentially cyclical and with a remarkable seasonal control, or imposition of external loads, which is continuous with variable rates.
The above-cited literature demonstrates the potential of A-DInSAR techniques for analysing subsidence processes over wide regions, thus giving an overall look at the affected regions and allowing the zonation of the effects in terms of deformation rates. Nevertheless, the possibility of making quantitative use of such data is not at all obvious. Indeed, in the analyses of a single building or a neighbourhood (feasible thanks to high spatial resolution satellite images), in some cases a sort of badly organized pattern of displacement rates over the same structure is observed independently from the geological complexity of the subsoil. This behaviour is caused by potential bias in the interpretation of the deformation, especially in terms of cause-and-effect relationships.
From the perspective of end-users of A-DInSAR data who infer information about the subsidence process in urbanized regions, we think that the main limits to overcome are intrinsic to the PSI technique and in the complexity of the physical process of subsidence, especially if the physical process is triggered by external loads (e.g., urbanisation) and regarded as a consolidation process developing over time. In fact, the PSI technique is strongly based on the presence and "quality" of the signal amplitude and the coherence of persistent scatterers (PSs) primarily featured in structures and infrastructures in an urban context. On the one hand, the presence of human settlements guarantees a reasonable number of PSs that can provide spatially dense information; however, on the other hand, the analysis of PSs lying on structures and infrastructures addresses the deformational response of buildings to a more complex process of subsoil settlement triggered by their load. The behaviour of the PSs on buildings can be regarded as a combination of factors, such as geological and geotechnical properties of the subsoil, the entity of the stress variations (either imposed loads or pore pressures), and the load entity and modalities of emplacement (i.e., types of foundations).
With the present case study, our intention is to exemplify a typical subsidence issue related to external loads (urbanisation) and a local geological and geotechnical setting (i.e., a subsoil composed of alluvial and coastal compressible deposits). The study evaluates the quantitative contribution of satellite SAR interferometry over a wider perspective of land-use planning and risk management, trying to answer questions such as the following examples: Can A-DInSAR help to reconstruct the consolidation process that usually presides over the subsidence? Is it possible from PS data to understand the evolutionary stage of the consolidation process that involves a structure/infrastructure in the subsidence process? Is it possible to predict the future behaviour of the terrain-foundation-structure combination?
The study area, located in the municipalities of Rome and Fiumicino (Italy), is affected by important subsidence phenomena [23], features a high geological and geotechnical variability, and has been built upon for approximately 100 years. Moreover, in the region, a large urban expansion is in progress and additional expansion is already planned. These features allowed us to process, analyse, and compare A-DInSAR data in light of the different combinations of geological structures, building types (loads and foundations), and ages. We adopted an integrated approach reflecting the complexity of the process that can be regarded as a "multivariate" issue. Then, the basic idea was to rely on a sort of stratified sampling of the displacement data inferred by the PS time series that was obtained by means of A-DInSAR analyses. The population of PSs, previously selected on a pure radiometric basis (i.e., signal amplitude and coherence), was split into subsets homogeneous in terms of structural features (foundation type) and class of construction age, which can provide a proxy of the "t0" of a consolidation process. These derived datasets can then be analysed to better understand and constrain the subsidence process, especially in terms of cause-and-effect relationships with the local geological setting (i.e., presence, depth, and thickness of compressible soil layers). Indeed, the geotechnical model of the subsoil, summarizing its mechanical characteristics and the related geometrical distribution, remains the only controlling variable on the magnitude of the process.
Finally, from the perspective of quantitative applications, we compare the time series of selected PSs with the consolidation curves obtained by applying mono-dimensional Terzaghi's theory.

Study Region and Geological Setting
The investigated region (approximately 250 km 2 ) is located southwest of the city of Rome. From a physiographic point of view, the region is dominated by the Tiber River delta and the coastal plain, both of which are characterized by a very flat morphology and an original marshy environment. This region belongs to the municipalities of Rome and Fiumicino between the Ponte Galeria Hills and the mouth of the Tiber River along the Tyrrhenian Sea ( Figure 1).
The region has been settled since the Palaeolithic period [24] during the time of ancient Rome and, more recently, during the reclamation of the ponds of Ostia and Maccarese in the late nineteenth and early twentieth century (Figure 1a). The reclamation works, which started in the late 1800s, consisted of the construction of the banks of the Tiber River, excavation of the two main channels ("Collettore delle Acque Alte" and "Collettore delle Acque Basse"), and the construction of a series of networks of orthogonal channels [25]. In the early 1900s, several farmhouses were built and regularly distributed in the region to service the reclamation phase.
In the middle 1900s, the presence of abundant undeveloped land near to the city of Rome influenced the location of the Leonardo Da Vinci International Airport in the northwestern part of this coastal plain and the subsequent urban expansion towards the sea; this expansion has been particularly intense over the last 30 years as recorded in this paper.
The Tiber River delta has been interpreted to be wave dominated [31,39,40]. Two main sectors have been distinguished: an outer delta, formed by beach ridges, dunes, and sand deposits; and an inner delta, composed of the typical peat and organic rich clay deposits of a lagoonal basin (Figure 1d).
The underlying bedrock consists of Lower and Middle Pleistocene overconsolidated clay: silty/clay shelf deposits and sandy/gravel fluvial deposits. The sandy/gravel fluvial deposits belong to the Ponte Galeria Sequence (PGS). The top of these units is characterized by an erosive surface (unconformity, Figure 1b) that constitutes the base of the TDS. The complex geological architecture of the TDS is the result of the delta's evolution over the past 20,000 years and is extensively summarized in the geological sketches reported in Figure 1c.
According to Milli et al. [39,40], the main bodies of compressible layers of the Tiber delta were deposited during the transgressive and highstand system tracts, when the sedimentation was characterized by a barrier-island system and lagoonal basin. In the last 5000-6000 years, this lagoon became the two marshy coastal ponds mentioned above ("Stagno di Maccarese" and "Stagno di Ostia"), which remained active up to the 1884 reclamation [25].
The hydrogeological setting of the Tiber delta features a deep and artesian main aquifer located in sandy gravels at the base of the TDS; it is sustained at the base by the lower Pleistocene clays acting as an aquiclude and in the upper part partially sealed by the low-permeable silty clay and clay lagoon deposits. Another aquifer is located in the sandy deposits of the TDS; its piezometry has been reconstructed on the left of the Tiber River [41,42]. The piezometric distribution shows that the aquifer top is below sea level in most of the region. The maximum depression is located just south of Ancient Ostia, with elevations reaching −5 m above sea level (a.s.l.). These values are linked to the increased use of pumping wells associated with increasing urbanization [41,42]. Recently, in the Fiumicino region, several underground degassing events occurred during drillings performed for different purposes at sea and inland [43][44][45]. The region of Fiumicino is affected by natural degassing, which was recently monitored by systematic measurements of ground gas emissions. In this context, Bigi et al. [45] illustrated the spatial correlation between deep faults and regions of anomalous concentrations of inorganic, deep-origin CO2 and CH4.
In this context, many recently built structures and infrastructures suffer problems related to differential settlements. This is, for example, the documented case history of airstrip №3 of the Leonardo Da Vinci airport [46], the Rome Fair region, and the Interporto Romano as well as some residential buildings. It has been inferred that the differential settlements are related to the presence of compressible strata in the subsoil belonging to the TDS. The more compressible geological units inside the TDS are the silty clay and clay lagoon deposits, the peat clay lagoon deposits, and the clay alluvial deposits followed by the sand dune deposits and the silty sand beach ridge deposits as classified by Milli et al. [39,40].  [36]). LST, lowstand system tract; TST, transgressive system tract; HST, highstand system tract.

Methods and Materials
In this study, the data referring to the subsoil, i.e., the geological setting and geotechnical characterization, and the data referring to the urbanisation, i.e., the ages and foundation types of the structures and infrastructures, are combined to explain the displacement data obtained via A-DInSAR analyses. The workflow is summarized in Figure 2.

Figure 2.
Workflow of the proposed approach to assessing ground response to urbanisation by combining geological and geotechnical data with satellite interferometric synthetic aperture radar (InSAR) data. PS, Persistent Scatterer; 2D, two-dimensional; 3D, three-dimensional; GIS, geographical information system.
The central pillar of the activity is the reconstruction of the subsidence record for the study region in recent decades. ERS, ENVISAT, and high-resolution COSMO-SkyMed satellite data were used, thus deriving the spatial evolution of the subsidence processes in the time frame 1992-2015. Specifically, ERS and ENVISAT PS data, spanning 1992-2000 and 2002-2010, both in ascending and descending geometries, were provided by the Ministry of Environment, Land, and Sea (MATTM), while a COSMO-SkyMed StripMap HIMAGE data stack (46 acquisitions covering the period from 2011 to 2015) in ascending geometry was provided by the Italian Space Agency (ASI) in the frame of an "open call for science" project and processed using PSI techniques for the specific purpose of this project (Table 1). The second pillar of our methodology was the definition of high-resolution lithostratigraphical and geotechnical models. The main objective was to define the spatial arrangement of the more compressible lithological units and the related geotechnical parameters. More than 480 stratigraphic logs from geognostic boreholes were collected in a depth range between 20 and 80 m ( Figure 3). A geotechnical dataset consisting of site tests (SPT, CPT) and laboratory tests on samples was also acquired as well as piezometric data. All these datasets cover a large part of the investigated region with a non-uniform distribution. The third pillar concerns the urbanisation of the study region in recent decades to identify the load emplacement start period of every building as a presumed proxy of the beginning of the settlement record. It has been identified by change detection analyses using topographic maps and optical images obtained from the Ministry of Environment, Land, and Sea (MATTM), the Italian Military survey office (IMG), and the municipality cartographies listed in Table 2.
The results obtained by the three main actions described above were preliminarily combined in a geographical information system (GIS)-based geo-spatial framework aimed at properly interpreting the spatial and temporal correlations between ground deformations and each supposed triggering/control factor of the subsidence process on a large scale (i.e., the scale of the Tiber River delta plain). Afterwards, further analyses were conducted to study in detail and quantify these relationships. For this purpose, we first selected some pilot regions based on the quantity of available data and, thus, the resolution of the reference geological/geotechnical models. Specifically, the choice criteria for selecting the pilot sites were (1) a high density of geological/geotechnical and PS data and (2) the representativeness of the geological, subsidence, and urbanisation context with respect to the entire study region. As a result, we focused on eight sub-regions identified within the original study region. As a second step, we partitioned the ground deformation datasets (i.e., the results of the central pillar) into "thematic" subsets, i.e., PS data were grouped in subsets differentiated on the basis of the age (11 classes, see Section 3.3) and the foundation type (2 classes, shallow and deep) of the structures on which they are located. This step allowed us to better constrain on a small scale the relationships between entity and rate of ground and structures deformation and the characteristics of the subsoil. To strengthen the inferred relationships, we selected a number of PSs and the related time series of displacements that are (i) representative of different classes in terms of the structure age and foundation type and (ii) located in regions where the subsoil setting allows for the application of Terzaghi's theory of one-dimensional consolidation [47]. This theory was then used, and the resulting theoretical consolidation curves were compared with the subsidence curves inferred by the PS time series. Table 2. List of topographic maps and optical images used to perform manual change detection analysis (sources: I.G.M., Italian Army's Geographic Institute; C.T.R., official topographic map of the Lazio region; P.C.N., "Geoportal of Italian Ministry of Environment, Land, and Sea").

Type
Acquisition

SAR Data
For the present case study, ERS and ENVISAT results were collected from the Ministry of Environment, Land, and Sea. They had already been analysed through A-DInSAR techniques and provide information on the entire deltaic plain. ERS and ENVISAT were processed using PSInSAR and PSP-DIFSAR techniques in the frame of the Italian National InSAR project financed by MATTM. Furthermore, PSI analyses were performed to acquire information about recent deformations that affected the ground surface using a stack of 46 COSMO-SkyMed images in ascending orbital geometry (spanning 2011 to 2015). Ground deformation measurements were obtained using the PS-InSAR technique [12,16,17] and proprietary procedures were implemented in SARPROZ [48], a software developed for multi-image InSAR analyses with the PS-InSAR technique.
Large-scale analyses were performed on a portion of the SAR images frame of approximately 15 × 15 km 2 . For this study, a conventional approach was adopted. Specifically, to reduce the decorrelation effects, all images were related to a single image, which was selected by considering the normal and temporal baselines. After all images were co-registered to the master image (1 November 2013), the radar parameter maps were generated: the reflectivity map (a multi-temporal amplitude value for each pixel) and the amplitude stability index map (the coefficient of variation of the amplitude). These maps were used as quality estimators for the selection of PS candidates (PSCs) in the PSI workflow. Then, the standard "star" graph was used ( Figure 4) to connect all slave images to a single master image to generate interferograms. A network of PSCs was created to estimate preliminary height and displacement parameters. This step is needed to detect and remove the atmospheric phase screen (APS) starting from the residual phase components. After APS estimation and removal, a second parameter estimation was performed on a much larger set of points, which was selected on a spatial coherence and amplitude stability index combination criterion as the final step of the PSI procedure. A digital elevation model with 30-m resolution was used to compute the differential interferograms (e.g., to subtract the topographic phase component from the interferometric phase) and to geocode the PS results. At the end of the PS analyses, all PSs with a coherence above 0.65 were selected. For each PS, the line of sight (LOS) velocity, displacement time series, and height were computed (relative to a reference point identified in a stable region outside the Tiber River delta plain).

Geological and Geotechnical Modelling
Stratigraphic information, such as borehole logs from many different sources (such as bibliographies, local companies, the Fiumicino Airport technical office, municipalities, and the National Roads Department) were used for the geological characterisation of the eight sub-regions selected for investigation. As a first step, to detect the more significant lithological features and related variations in the subsoil, it was necessary to harmonize information coming from many different sources: the stratigraphic descriptions were standardized by attributing a unique code according to a lithological criterion, which implicitly refers also to the technical properties.
The re-interpreted stratigraphic logs were included in a geodatabase that provides, for each borehole, information on the description, thickness, and depth (in meters above sea level and with respect to ground level) of the lithological units. For each lithological unit, the geotechnical parameters were defined using the results of many field and laboratory geotechnical tests derived from a bibliography of foregoing geognostic surveys. Some lithological units were grouped if similar geotechnical behaviour resulted and, in contrast, some lithological units were subdivided if different geotechnical behaviours were inferred.
This database was the cornerstone for obtaining three-dimensional (3D) geological and geotechnical models, which were generated starting from punctual information, and then interpolated via the inverse distance weighting (IDW) algorithm in the focus sub-regions using the RockWorks16™ software package (RockWare, Inc., Golden, CO, USA). The most useful products of the 3D models for our purposes were geological cross sections, maps of the thickness of compressible lithologies, and maps of selected surfaces by, for example, iso-depth lines.
Furthermore, bibliographic research was performed to obtain a detailed or indicative description of the foundation type of some structures. As the final stage of this phase of the study, in the most interesting region in terms of heterogeneous age of urbanisation and ground deformation behaviour (specifically the Rome Fair and Commercity region), the one-dimensional (1D) consolidation theory was adopted [47] to obtain theoretical settlement rates in the region due to load emplacement.

Urbanisation
Our attention was focused first on the age of urbanisation as, for a region surrounding one or more coeval buildings, it represents the well-known time "t = 0" for the consolidation process induced by an external load causing their settlement.
The urbanisation phases were rather discontinuous and inhomogeneous from a temporal and spatial point of view, especially in the period after the reclamation. The building construction dates in the study region range from circa 1900 (e.g., the numerous farmhouses spread across the entire delta plain built during the reclamation phase) to the present day (e.g., the industrial and commercial structures such as the Da Vinci shopping mall and the Rome Fair and Commercity region). To derive the ages of the buildings inside the eight sub-regions, we conducted a manual change detection analysis to minimize errors coming from automatic processes, identifying for each homogeneous structure or group of structures among the images listed in Table 2 the one in which it appears for the first time and creating a corresponding polygon in a GIS-based platform. Eleven "time intervals" were distinguished and each polygon was classified according to this time scale. A total of approximately 350 polygons were created in the sub-regions of interest. Figure 5 contains a sequence of images to show an example of the urban expansion that occurred during the monitoring time.

Data Post-Processing and Integration
ERS, ENVISAT, and COSMO-SkyMed results were combined with geological and geotechnical data and information regarding urbanisation (building ages and foundation types) in a GIS-based platform. A geodatabase of the stratigraphic and geotechnical data was generated containing the location of every borehole, the top/bottom height (in m a.s.l.) of every lithotype, and their geotechnical parameters. The two-dimensional (2D) thicknesses contour maps of the compressible layers derived by Rockworks software were geocoded and imported into the GIS platform. Moreover, the results of the change detection analysis of the urbanisation were imported in the form of polygons with the building age attribute. The same task was performed for the information about the foundation types of the structures. A spatial overlay between these data was performed: the entire achieved A-DInSAR dataset was sampled and partitioned into homogeneous subsets based on relevant features of the structures on which they are identified, such as (i) soft-soil thickness; (ii) type of foundation; and (iii) age of structure, allowing us to obtain statistical correlations between ground deformation and the presence and thickness of compressible soil layers.

Results
The results of the large-scale A-DInSAR analyses for ERS, ENVISAT, and COSMO-SkyMed are reported in Figure 6, in which they are plotted on the map by Amenduni [25], which is a reference map depicting the region in the year 1884, before urbanisation.
Different numbers of PSs were obtained among the satellite missions. Approximately 25,000 PSs were obtained from the ERS ascending and descending data stacks, about 110,000 from the ENVISAT ascending and descending data stacks, and finally approximately 210,000 from the ascending stack of COSMO-SkyMed (covering a reduced area as showed in Figure 3). This difference is due to the different spatial and temporal resolutions of the satellite missions and to the period of acquisition of each satellite mission: for the more recent COSMO-SkyMed, it was possible to catch a larger number of PSs also due to the intense urbanisation that had already accumulated during the acquisition period.
A 22-year deformational history of the region has been derived by combining the ERS, ENVISAT, and COSMO-SkyMed datasets. The PSs show that a large portion of the Tiber Delta Plain is affected by relevant deformation, with displacement rates along the LOS of up to 30-40 mm/year away from the sensor. Relevant displacements approaching the sensors are absent and the measured deformation is primarily due to vertical displacements and subsidence/settlement processes with a good correlation with the geological setting. The main evidence is at the large scale ( Figure 6) where the deformational pattern highlights that the outer delta portion (Figure 1d), characterized by coarse and low compressible soils, is affected by slight or negligible movements, whereas the inner delta portion (Figure 1d), characterized by thick lagoonal/marshy compressible soils, is affected by a large deformation, clearly detected via the A-DInSAR technique. Furthermore, in the northern sector of Figure 6, a region affected by very low deformation is observed: it is located at the foot of the hills bordering the coastal alluvial plain. Here, the subsoil is composed primarily of over-consolidated clays and sandy soils.
The geological and geotechnical differences between the inner and outer deltas are detailed in Figure 7: the two representative small-scale 3D geotechnical models refer to the inner delta (3 km 2 wide in Figure 7a) and outer delta (1.3 km 2 wide in Figure 7b). Their locations are reported in Figure  6f (Sector 1 and Sector 2, respectively). In general terms, the geological reconstruction highlights complex vertical and lateral lithological variations. In some cases, very sharp lateral facies heteropies are detected, implying a complex interpretation of ground response to urbanisation.
In the inner delta region (Figure 7a), a low-compressibility stratum (GU1) constituted by hard-clay overlays a 30-m thick body characterized by medium to very highly compressible deposits: organic clay, peat, and silty clay (GU2, GU3 in the shallow part, and GU4 in the deep part). A more detailed reconstruction is reported in Figure 8. Table 3 shows a synthesis of the corresponding geotechnical parameters. A number of buildings in this region suffered considerable damage due to differential settlement. The geotechnical cross sections reported in Figure 8 show the geometric relationships among the geotechnical units, highlighting the strong lateral and vertical variability also in a small-sized region.  . Geotechnical 3D models of (a) sector 1 (Rome Fair region and Commercity, inner delta) and (b) sector 2 (offices of the Leonardo Da Vinci airport, outer delta). The figure also shows the processing chain for the 3D geotechnical model, which started from the collection and interpretation of borehole data (first images on the top), followed by their correlation as fence diagrams, and finally interpolated as continuous surfaces (last images at the bottom).  The region shown in Figure 7b refers to the outer delta. Its geological and geotechnical setting is characterized primarily by sand deposits, from fine to coarse-grained and from loose to dense, and by the presence of centimeter-thick silt levels or lenses of fine-grained, more compressible soils.
In Figure 9, the different deformational patterns of the two sectors belonging to the inner and outer deltas, due to the different subsoil characteristics, are clearly visible.  Figure 6f) and outer delta (right column, sector 2 in Figure 6f) regions. Figure 10 shows the foundation type (shallow or deep) and the building ages derived from acquisition of technical information and manual change detection. In both sectors, it is possible to observe the presence of buildings with shallow and deep foundations; moreover, the pavement has been distinguished. A few buildings were constructed in approximately 1900 and other buildings were constructed in the last few years.
To focus on the role of the geological and geotechnical setting of the subsoil on the local scale, the airstrip №3 of the Leonardo Da Vinci Airport is analysed (Figure 11). The airstrip lies on an abrupt lithological variation: the subsoil of the southern part of the airstrip is primarily characterized by an approximately 30-m thick body of highly compressible organic silty clay, whereas the northern sector is characterized by a low-compressibility silty sand body [46]. The airstrip suffers large problems caused by the differential settlements of these two sectors and some important remediation works have been conducted over time.
The time series of displacement obtained by merging the data processed from the three satellite missions confirm that the A-DInSAR data can "read" clearly these sharp differences in the subsoil.
The accuracy is such that the high-resolution spatial data from the COSMO-SkyMed mission also reveals more detail, such as the presence of a shallow organic clay body included in the sandy unit of the northern sector (marked by an orange triangle in the time series of the displacements), locally causing an increase in the displacement trend of the northern sector of the airstrip. Returning to the large scale of Figure 6, over the monitoring time interval (1992-2015), an apparent constant spatial distribution of the deformation rate is observed, with the exception of some regions where urbanisation has occurred during this time, namely, the reclaimed ancient Maccarese and Ostia ponds identified by Amenduni [25] and the Piana del Sole portion located east of the Maccarese pond, where the urbanisation was conducted in the period between the ERS and ENVISAT missions. In the region urbanised over the ancient Ostia ponds, a decrease in the deformation rate is observed in the more recent monitoring period. Figure 11. Ground deformation of airstrip 3 of the Leonardo Da Vinci airport analysed by ERS ascending (a); ERS descending (b); ENVISAT ascending (c); ENVISAT descending (d); and COSMO-SkyMed (e). In (f) are reported the cumulative time series of displacement of the three main sectors detected at the airstrip. In (g) a geological cross section of the airstrip is shown (from Manassero and Dominijanni [46], redrawn). The location of this region is reported in Figure 6f.

Data Integration and Discussion
A qualitative general sounding relationship between the ground deformation spatial pattern and the geological and geotechnical setting has been shown in the previous paragraph. Unfortunately, if we try to strengthen this relationship with a more rigorous approach, some surprises are encountered. In the region of Figure 12a, the velocity of displacement of all PSs is plotted against the thickness of the more compressible soils (GU2 and GU3) derived by the geological and geotechnical models: no relationships can be demonstrated. Moreover, along the section A-A', a large total displacement is expected in the northern part with respect to the southern one due to the greater thicknesses of the GU2 and GU3 layers: the displacement data seems to contradict this working hypothesis. Some improvements are needed to infer clues to interpreting the acting processes. This is the reason why, in our opinion, the use of InSAR data and analyses for detecting regions featuring a high thickness of compressible soils, as performed by Del Ventisette et al. [49], can be reliable only for a rough large-scale characterization, but is potentially misleading for an accurate interpretation of the subsidence process if other relevant factors, such as foundation type and age of construction, are not considered.
Focusing on the Rome Fair/Commercity region, the first improvement of the above-discussed sounding clue is obtained if the type of foundation is considered, as the reader can appreciate by observing the white rectangle in Figure 12 (in the northern sector of the A-A' section), where buildings primarily with pile foundations some tens of meters deep are located (see Figure 10c). Because these deep foundations do not interact with the recent compressible deposits, their corresponding deformation is significantly low also with respect to the surrounding pavement ( Figure 12). As a consequence, their displacement rate cannot be compared with the one estimated for the buildings arranged with shallow foundations interacting with the compressible units GU2 and GU3 and affected by higher displacements. In other words, in Figure 12b, PSs corresponding to buildings with deep foundations should be eliminated.
In Figure 13, the attention is focused on buildings with shallow foundations and on the highway embankment and pavement, all interacting with the shallowest and more compressible part of the TDS. By reconstructing the time history of displacement (Figure 13c) combining the three A-DInSAR datasets, a displacement time series lasting more than 20 years is obtained. As one can see, from the geometric point of view, it is very similar to the one expected by the theory of consolidation, similarly to what was observed also by Stramondo et al. [15] even if with a smaller time series. This result is opening interesting perspectives about the potential of A-DInSAR-derived time series as a tool for predicting and calibrating the settlement trend caused by the construction of new buildings. Prediction capabilities based on A-DInSAR-derived time series have been recently hypothesized for landslides by Moretto et al. [50].
The introduction in the selection criteria of the age of the shallow-founded buildings improves more significantly the statistical relationships we are looking for. This is demonstrated by the plots of Figure 14b,c, where the velocity of displacements of these buildings estimated by the analysis of the highest spatial resolution dataset from the COSMO-SkyMed mission is plotted versus the thickness of the compressible units GU2 and GU3 separately for the buildings constructed during 1994-1998 (Figure 14b) with respect to those constructed during 1998-2002 (Figure 14c). From a theoretical point of view, this separation has been introduced to consider the nonlinear trend of the consolidation process over time. Therefore, because a decrease in the displacement rate over time is expected during the consolidation process (see Figure 13), in a rigorous analysis the displacement velocities of differently aged buildings estimated in the same period cannot be compared or mixed to infer other information, such as the thickness of the compressible units.
Finally, if the behaviour of buildings constructed during 1994-1998 and 1998-2002 is examined separately, the relationship for the deformation rate (measured in the timeframe 2011-2015) versus the thickness of the compressible units GU2 and GU3 is statistically reliable for buildings with shallow foundations and pavements (Figure 14b,c).  The relevance of the "thickness of compressible units" parameter is also inferred by analysis of single buildings with shallow foundations using the COSMO-SkyMed PS dataset ( Figure 15). The geotechnical cross section B-B' shows an increase in the compressible soil thicknesses (GU2 and GU3 units) that generate differential deformation to the structures. Bivariate analyses show a good relationship between the velocity of the persistent scatterers and the compressible soil thickness.
However, the same correlation is weaker over another geotechnical cross section C-C'. We suppose that at the scale of single buildings, as in this extreme case, some more details should be considered, probably related to the design and construction details. Consequently, at the scale of a single building, we can also conclude that a detailed definition of the geological setting of the subsoil, knowledge of the type of foundation, and the period of constructions are not enough to explain in all the details the distribution of the deformation rate. These results are demonstrating that, in most of the cases, the use of A-DInSAR as a stand-alone tool for risk assessment of buildings affected by subsidence is still far from being considered a feasible option. However, it has the potential to support analyses performed by conventional building damage analyses [51] as recently demonstrated by Milillo et al. [14]. Regarding the importance of the age of the building for the subsidence analysis, a great opportunity is offered by the diffuse presence of the approximately 100 buildings (ancient farmhouses) built during the reclamation period of the Tiber delta coastal plain, i.e., more than 80 years ago ( Figure 16). These rural structures are (1) built at similar dates and with similar procedures, i.e., similar loads transmitted to the subsoil; (2) uniformly distributed as a regular pattern over the delta plain in different geological conditions; (3) all arranged with a shallow foundation; and (4) ancient, so that the expected settlement induced by these structures is null or very low. The deformation rates of these farmhouses, with the exception of those that have been renovated or influenced by the emplacement of new loads in the near surroundings, indicates an overall absence of relevant deformation independent of whether they are located on the outer or inner delta (approximately −1/−1.5 mm/year away from the sensor). In the inner delta region, where we estimate a very high subsidence rate for the recently built structures and infrastructures as described in the results section, these farmhouses show very low displacement. They can be related to the last phase of the consolidation process or to the secondary consolidation process induced by their load or, in a larger vision, to a "background subsidence", for example induced by secondary consolidation processes or by tectonic causes, or to a combination of effects (natural subsidence as defined by Tosi et al. [11] for the Venice area). However, the measured displacement values are close to the A-DInSAR accuracy. Figures 6, 9 and 11-15 demonstrate that A-DInSAR analysis can read the subsidence process induced by urbanisation at large and small scales, especially if high-resolution SAR data are used. The spatial and temporal distribution of the deformation inferred by this powerful tool allowed us to interpret realistically the process, considering from one side its intrinsic complexity and from the other side the characteristics of the object that we used to measure its development. Because this process should follow the well-known 1D consolidation theory as inferred by Figure 13, in the inner delta region we attempted to test if the estimation obtained by applying this tool agrees with the theory. We chose the inner delta sector of the study region where the Rome Fair and Commercity are located. We selected two representative borehole logs (S1 and S2 in Figure 7a) and the corresponding geotechnical parameters of the recognized units (Table 3). We considered the load induced by standard industrial structures (stresses ranging from 80 to 150 kPa, spot and grade beam foundation) and by industrial and road pavement (stress of 26 kPa; plate foundation) and we estimated (i) the total settlement induced by oedometric theory and (ii) the time to deplete 50% and 100% of the induced total settlement, denoted t50 and t100, respectively. The results are shown in Table 4: Table 4. Theoretical settlement entity and rates, estimated from borehole S1, using standard loads, collected foundation types, and geotechnical parameters derived from 3D modelling.
We roughly observed a generally good agreement with the values measured via A-DInSAR analyses: displacement of shallow-founded buildings and pavement on the order of decimetres in the inner delta sector and lengthy related consolidation processes lasting tens of years.
Finally, we compare our estimations of settlement of buildings on shallow foundations and pavement (built on different ages) with the corresponding displacement measured by COSMO-SkyMed in the period 2011-2015. Specifically, three classes of buildings, i.e., 5 to 15, 15 to 30, and >80 years old and one class of pavement 0 to 5 years old with respect to the COSMO-SkyMed time-frame (2001-2015) were selected. Figures 17 and 18 show the very good and surprising agreement between the theoretical consolidation values and the PSI measured settlement (once projected along the vertical direction). Such a good correspondence looks promising in the attempt to perform a simplified evaluation of settlement induced by construction at the large scale if good geological and geotechnical information are available.

Conclusions
At a first glance, the pattern of ground deformation inferred by A-DInSAR analyses in the study region roughly agrees with the spatial variations in presence, depth, and thickness of soft soils. However, if examined at high resolution, this correspondence is revealed to be quite weak because several "outliers" can be observed. PSs moving rapidly downward are sometimes located in favourable geotechnical conditions (i.e., where the subsoil is primarily composed of low compressibility deposits). In contrast, steady or slow-moving PS can be found in unfavourable geotechnical conditions, where the subsoil features thick layers of soft sediments. This can be primarily attributed to the deformational response of structures where the PS lay, which is not only due to the subsoil's geotechnical properties but also to the type of foundation and the age of construction. In the case of subsidence due to imposed loads, the subsidence can affect the actual layer interacting with the extra load whereas the imposed load regulates the rate of subsidence according to the nonlinear, time-dependent consolidation process.
Hence, for effective risk mitigation purposes, an effort should be made to properly manage and process A-DInSAR data. In our experiences, the most important and decisive step consists of partitioning the entire PS dataset into homogeneous subsets based on relevant features of the structures on which they are located, such as (i) the type of foundation and (ii) the age of construction. This procedure allowed us to analyse the cause-and-effect relationship between subsidence rates and the constitution of the subsoil starting from consistent and comparable information. Furthermore, once the portioning has been performed, it is possible to better constrain the subsidence as measured by A-DInSAR, which is a powerful tool for quantifying the process. The PS time series can be considered a proxy of the consolidation process: the soft-soil thickness and the foundation type control the subsidence magnitude, whereas the age of construction governs the subsidence rate. This evidence is strengthened by the good fitting between the experimental consolidation curves (i.e., those achieved by PSI back-analyses) and the theoretical ones (i.e., those obtained by applying mono-dimensional Terzaghi's theory where appropriate). As a consequence, it is possible to state that, if properly managed, A-DInSAR analyses-especially if performed with high spatial resolution imagery-can help to reconstruct the consolidation process that usually presides over the subsidence, thus allowing the evolutionary stage of the consolidation process involving a given structure/infrastructure to be captured. This evidence suggests that remotely sensed data can therefore be useful for forecasting the future behaviour of a complex "system" made of terrain-foundation-structure, without prejudice to other relevant controlling factors, such as shallow water table fluctuations. As regards controlling factors, effects of possible localized groundwater withdrawals could be related to "punctual anomalies" in the subsidence trend, but in general terms it is possible to state that the overall fluctuations of the main shallow aquifer, whose piezometric level is constantly a few meters below the ground level as testified by the still-active drainage channels related to reclamation works, are confined to seasonal variations in a very narrow range.
In conclusion, the results we obtained provided useful clues to elucidate the potential of A-DInSAR for subsidence risk mitigation. However, its reliability is strictly related to the scale of analysis (the small scale is essential), and to the knowledge of boundary conditions that in turn requires an integrated approach consisting of detailed geological and geotechnical modelling, data acquisition and storage, and multi-temporal analysis of land cover.