Regression Analysis of Subsidence in the Como Basin (Northern Italy): New Insights on Natural and Anthropic Drivers from InSAR Data

Natural and anthropogenic subsidence such as that in the Como urban area (northern Italy) can cause significant damage to structures and infrastructure, and expose the city’s lakefront to an increasing risk of inundation from Lake Como. This phenomenon affecting the Como basin has been studied by several researchers, and the major drivers of subsidence are known. However, the availability of historical InSAR data allowed us to reconsider the relationship between subsidence predisposing factors (i.e., the thicknesses of reworked and compressible layers, overburden stress, and the piezometric level) and ground surface displacements with higher precision over the entire basin. Benefiting from the deep knowledge of the hydromechanical setting of the Como basin and the availability of InSAR measurements from 1992 to 2010, in this paper we model subsidence-related movements using linear and nonlinear regression methods in order to determine the combination of natural and anthropic factors that have caused subsidence in the Como basin over the past decades. The results of this study highlight peculiar patterns of subsidence that suggest the influence of two further causes, namely tectonic control of the sedimentary architecture and diversion of local streams, which have never been considered before. This analysis aims to assess the spatial distribution of subsidence through InSAR analysis in order to enhance the knowledge and understanding of the phenomenon in the Como urban area. The interferometric data could be used to better plan urban risk management strategies.


Introduction
Land subsidence is the progressive lowering or sinking of the ground surface [1], which occurs primary in lowland, coastal, and delta areas, and is caused by various natural or anthropogenic processes [2][3][4][5]. These triggering causes are often investigated independently from each other, although land subsidence may result from their combination or interaction, thus amplifying the subsidence-related hazards in urban areas [6]. However, only a few studies integrate diverse causes when interpreting the pattern of subsidence at large scales [6][7][8].
Traditionally, land subsidence is monitored via in situ techniques, for example using GPS or geodetic leveling [9,10], which allow pointwise, precise, and reliable measurements of the ground subsidence rates [5]. These instruments are usually positioned on a few selected benchmark areas belonging to local or regional networks because of the cost of their installation and maintenance, thus resulting in a lack of measurement density over large areas [5]. The ensuing low spatial resolution of traditional measurements may be insufficient to capture the small-scale variations of ground movements in subsiding areas, which are usually heterogeneous and can extend for hundreds of square kilometers. As a consequence, underestimations of settlement rates are likely to occur when relying exclusively on traditional measurements over large areas.
During the last decades, Interferometric Synthetic Aperture Aadar (InSAR) techniques have provided a valuable alternative for monitoring land subsidence, benefiting from having a relatively low cost, millimetric accuracy, wide area coverage, short revisiting time, and from the availability of historical datasets [11][12][13]. The standard InSAR approach based on single interferograms [14] has been rapidly overcome by more advanced techniques exploiting stacks of Synthetic Aperture Radar (SAR) images (multi-interferograms). Persistent Scatterer Interferometry (PSI) algorithms have been developed to estimate displacement time series and deformation rates over selected targets, defined as Persistent Scatterers (PS) [15]. The first PSInSAR algorithm was proposed by [16] to identify coherent radar targets with high reflectivity and phase stability over the entire observation period. To improve the spatial density of radar targets in non-cultivated or desert areas with moderate coherence due to small scattering objects, other techniques using Distributed Scatterers (DS) were then proposed, such as the Small BAseline Subset (SBAS) [17] and SqueeSAR [18] algorithms (for an overview, see [15]). Since their introduction, InSAR techniques have been largely applied for monitoring land subsidence, both in Italy [5,8,11,[19][20][21][22][23][24][25][26][27] and worldwide [2,6,12,13,[28][29][30][31][32][33][34][35][36]. The InSAR techniques allow the reconstruction of the deformation history of radar targets (PS) and for information about ground settlements to be obtained with high precision [23]. However, the interpretation of surface deformation signals may be challenging [29] in areas experiencing movements in the range of several millimeters to several centimeters. Thus, a deep knowledge of the local hydrogeological and stratigraphic setting is extremely important when interpreting InSAR-derived displacement maps [27] and investigating ground deformations induced by geological processes [11].
In this paper, we use historical InSAR data from 1992 to 2010 to investigate the spatial distribution of subsidence in the Como basin (northern Italy), induced by the combined action of several natural and anthropic causes in the investigated time range. The city of Como is naturally prone to subside because of a thick sequence (more than 180 m) of unconsolidated silty clay sediments with poor mechanical properties [37], which cause a long-term Late Pleistocene to Holocene subsidence rate of 1-2.5 mm/year [3]. Several previous studies have investigated the subsidence in the Como basin and contributed to the current good knowledge of the hydrogeological, geotechnical, and geomorphological properties of the basin [3,37,38]. This knowledge provides the basis for interpreting the recently available historical InSAR measurements, thus exploring the subsidence in the Como basin from a new perspective. By means of multivariate regression analysis, we aim to (i) establish which combination of natural and anthropic factors most likely determines subsidence in the study area; (ii) investigate the relationships between subsidence predisposing factors (SPF) and InSAR-derived ground movements; (iii) interpret the spatial pattern of subsidence in the Como basin over 18 years. This innovative analysis of the Como area allowed us to give new insights on possible natural or anthropic influences on subsidence, which have never been investigated before in this area.

Study Area
The city of Como is located on the southwestern shores of the lambda-shaped Lake Como (also called Lario; Figure 1a), which currently occupies a Pleistocene bifurcated glacial valley at an elevation varying from 197 to 210 m above sea level (a.s.l.). The Como branch of Lake Como has been hydrologically closed since ca. 18 thousand years (kyr) ago, making it a perfect sedimentary trap for fine-grained, often organic lacustrine and palustrine deposits with high sedimentation rates [37]. The city center lies on an alluvial plain drained by the Cosia and Valduce streams, and bordered by steep mountain slopes composed of Mesozoic pelagic carbonates to the NE (i.e., the Medolo Group-Early Jurassic) and by deep sea fan conglomerates to the SW (the Gonfolite Group-Oligo-Miocene) thrust on top of the Mesozoic units ( Figure 1c). The Gonfolite Group includes a thick sequence of conglomerates (Como and Villa Olmo Conglomerates) on top of sandstones and mudstones (the Chiasso Formation).
During the Quaternary, the morphology of the basin was repeatedly shaped by the erosional and depositional activity of glaciers that covered the area with a thickness of hundreds of meters of ice [39][40][41][42].
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 21 3 the erosional and depositional activity of glaciers that covered the area with a thickness of hundreds of meters of ice [39][40][41][42]. Sketchy geological cross-section of AA' (note vertical exaggeration); RM stands for reworked material, SG for sands and gravels, OS for organic silts, SD for sediments with dropstones, and CD for coarser deposits (modified from [37]). (c) Simplified geological map of the Como basin (modified from [42]). (d) Distribution of borehole logs (in orange) and piezometers (in blue), and outline of the geological cross-section.
Since 1998, our research group has studied the geological, environmental, and urban evolution of the Como area by systematically collecting stratigraphic, hydrogeological, geotechnical, and archaeological data [3,37]. Thanks to this extensive and long-lasting effort, we have reconstructed the three-dimensional sedimentary architecture of the Como basin, although only the upper ca. 180 m of the late Pleistocene to Holocene sedimentary sequence is known [37]. Specifically, Figure 2 shows the typical stratigraphy of the basin composed, from top to bottom, of: heterogeneous reworked material with archaeological remains up to 10 m depth (unit 1-RM); alluvial sands and gravels up to 15-24 m (unit 2-SG); palustrine organic and highly compressible silts up to 30 m (unit 3-OS) with sandy (unit 3a) or clayey facies (unit 3b); glaciolacustrine sediments with dropstones (unit 4-SD) and coarser proximal deposits up to 40-60 m (unit 5-CD). In the lakeshore area, silty and highly compressible subunits within anthropic sediments were also recognized [37]. The index properties of each stratigraphic unit are also summarized in Figure 2. (b) Sketchy geological cross-section of AA' (note vertical exaggeration); RM stands for reworked material, SG for sands and gravels, OS for organic silts, SD for sediments with dropstones, and CD for coarser deposits (modified from [37]). (c) Simplified geological map of the Como basin (modified from [42]). (d) Distribution of borehole logs (in orange) and piezometers (in blue), and outline of the geological cross-section.
Since 1998, our research group has studied the geological, environmental, and urban evolution of the Como area by systematically collecting stratigraphic, hydrogeological, geotechnical, and archaeological data [3,37]. Thanks to this extensive and long-lasting effort, we have reconstructed the three-dimensional sedimentary architecture of the Como basin, although only the upper ca. 180 m of the late Pleistocene to Holocene sedimentary sequence is known [37]. Specifically, Figure 2 shows the typical stratigraphy of the basin composed, from top to bottom, of: heterogeneous reworked material with archaeological remains up to 10 m depth (unit 1-RM); alluvial sands and gravels up to 15-24 m (unit 2-SG); palustrine organic and highly compressible silts up to 30 m (unit 3-OS) with sandy (unit 3a) or clayey facies (unit 3b); glaciolacustrine sediments with dropstones (unit 4-SD) and coarser proximal deposits up to 40-60 m (unit 5-CD). In the lakeshore area, silty and highly compressible subunits within anthropic sediments were also recognized [37]. The index properties of each stratigraphic unit are also summarized in Figure 2.  [37]). Extensive radiocarbon dating, which was also performed on other nearby boreholes [37,59], supports the notion that this stratigraphic column spans the last ca. 20 kyr, which is the post-Last Glacial Maximum time windows.  [37]). Extensive radiocarbon dating, which was also performed on other nearby boreholes [37], supports the notion that this stratigraphic column spans the last ca. 20 kyr, which is the post-Last Glacial Maximum time windows.
Besides the stratigraphic and geological setting determining a natural subsidence in Como basin, anthropic activities also played an important role over the last ca. 2 millennia. Since the founding of Como city in 59 BC, the coastline and water courses have been deeply modified through outstanding hydraulic planning, including the diversion of Cosia and Valduce streams performed by Romans in the second half of the I century BC [37,43]. More recently, between 1950 and 1970, the city experienced a ten-fold increase in the ground lowering rates due to extensive exploitation of water from a deep confined aquifer [3], land reclamation, and changes in the urban configuration of the city [37]. In 1976 the Como Municipality established a dedicated commission (called "Commissione Subsidenza") to study the subsidence phenomenon and assess its predisposing and triggering factors in a risk management perspective [44]. With the ceasing of the aquifer exploitation, the unconsolidated silty sediments experienced a slight recovery between 1980 and 1997, as highlighted by the uplift of several benchmarks located in the downtown area [3]. In 1996, the Como Municipality commissioned the construction of fixed and movable bulkheads, detention tanks, and jet grouting barriers along the lakefront [38], aiming to mitigate the occasional flooding that the city still experiences. This project was approved in 2004 and the construction began shortly after; nevertheless, in 2011 the subsiding rate of the lakefront suddenly increased due to the perturbation of the ground water regime, determining the suspension of construction of the anti-flooding facilities.

Materials
Ancillary thematic maps were collected from the Geoportal of Lombardia region [45] and consist of a topographic map (Carta Tecnica Regionale-CTR Lombardia) at a 1:10,000 scale and a digital terrain model (DTM) with 5 m spatial resolution.
Information about the hydrogeologic, stratigraphic, and geotechnical setting of the Como basin was directly acquired from University of Insubria or collected from Como Municipality and other private and public companies. We compiled 261 borehole logs (Figure 1d), reaching a maximum depth of 180 m. Logs have been recorded since the 1950s during drilling campaigns for public and private projects (e.g., building construction, water extraction, restoration of the Como lakefront). The depth of the phreatic water table is currently measured at 28 active piezometers belonging to the municipal network ( Figure 1d). Groundwater flows toward the lake, which represents the local base level; seasonal fluctuations of the piezometric level are of ca. 1.5 m near the lakeshore, which decrease moving from the lake to the SE part of the basin. At a distance of 300 m from the lakeshore, the seasonal oscillations of the piezometric level are less than 0.5 m. The hydraulic gradients in the urban area have been studied in detail since the 1970s and are published in the report by the "Commissione Subsidenza" [44].
The InSAR data (Table 1) were supplied by the Italian Ministry of Environment, Land, and Sea Protection via the "Piano Straordinario di Telerilevamento" [46]. Specifically, 102 images from ERS 1 and 2 (April 1992-December 2000) and 103 from Envisat (July 2003-July 2010) in both ascending and descending orbits were processed via E-geos PSP-DIFSAR technique [47,48], provided as a set of sparse points (called Persistent Scatterers (PS); Figure 3) with coherence greater than 0.6. Each point target has the following specific attributes: PS identifier, geographic coordinates, coherence, mean annual velocity along the line of sight (LOS), standard deviation, and displacement time series. These measurements are temporally dependent from the first acquired image and spatially relative to a reference point located on a bedrock site in Camerlata [3].

Methods
The collected data mentioned before were then used as inputs for a three-phase procedure ( Figure 4) consisting of: (1) data preparation, (2) data analysis, and (3) regression analysis.
As a preliminary step, we derived the slope map from the Digital Terrain Model (DTM; 5 m spatial resolution) to select the flat basin of Como (steepness < 5°) as the study area and discard the Persistent Scatterers (PS) recorded on areas with slopes greater than 5°. Assuming horizontal movements to be negligible in subsiding areas [22], the velocity measured along the satellite line of sight (LOS) of each PS was projected along the vertical direction to obtain the vertical component of movement and computed as follows [49]: where Vvert and VLOS (mm/year) are the PS velocity along the vertical direction and LOS, respectively; θ is the LOS incidence angle. From the borehole logs, we extracted the thickness of each stratigraphic unit and calculated the overburden stress imposed on the deeper compressible unit (i.e., unit 3-OS) as: where σv is the overburden stress (N/m 2 ) at depth z (m); is the soil unit weight (N/m 3 ); and Tn is the thickness (m) of the n-th stratigraphic unit. As in [37,38,44], the isopiezometric curves representing the mean level of the surficial aquifer in

Methods
The collected data mentioned before were then used as inputs for a three-phase procedure ( Figure 4) consisting of: (1) data preparation, (2) data analysis, and (3) regression analysis.
As a preliminary step, we derived the slope map from the Digital Terrain Model (DTM; 5 m spatial resolution) to select the flat basin of Como (steepness < 5 • ) as the study area and discard the Persistent Scatterers (PS) recorded on areas with slopes greater than 5 • . Assuming horizontal movements to be negligible in subsiding areas [22], the velocity measured along the satellite line of sight (LOS) of each PS was projected along the vertical direction to obtain the vertical component of movement and computed as follows [49]: where V vert and V LOS (mm/year) are the PS velocity along the vertical direction and LOS, respectively; θ is the LOS incidence angle.
7 We then rasterized the flat study area over 6207 cells of 20 × 20 m and uniformly interpolated both InSAR and stratigraphic point data with cells of the same dimensions in order to create a continuous surface from sample points at sparse locations. Specifically, we used ordinary kriging [50] to interpolate the InSAR vertical velocity (Vvert) and empirical Bayesian kriging [51] to interpolate the thicknesses of the reworked material (unit 1-RM) and organic silts (unit 3-OS), the σv value, and the piezometric level. Details on the performed interpolations are given in Appendix A. Different interpolators were necessary because of the nonuniform spatial distribution and density of the sample dataset. For each pixel (20 × 20m) of the study area, the thicknesses of unit 1 (RM) and unit 3 (OS), the σv value, and the piezometric level were summarized and considered as the possible subsidence predisposing factors (SPF).
In the second phase, we performed a bivariate correlation analysis to determine if input variables were linearly dependent using Pearson's correlation coefficient (PCC). Based on the PCC values, the bivariate linear correlation values were classified as null (PCC < |0.25|), weak (|0.25|< PCC < |0.50|), moderately strong (|0.50|< PCC < |0.75|), and strong (|0.75|< PCC < |1|). Specifically, we performed two correlation analyses: the first correlation was between the InSAR vertical velocity (Vvert) and SPF, which suggested the type of model to be used in the third phase; the second correlation was among SPF values, which although being a sufficient but not necessary condition, recommended on the need to check for multicollinearity. Any regression model can indeed be biased by multicollinearity, which can reduce the statistical significance of independent variables (i.e., subsidence predisposing factors) that are instead relevant to explain or predict the dependent variable (i.e., InSAR vertical velocity). To this end, we performed the Farrar-Glauber test (F-G test) [52], which consists of a set of three tests: (i) Chi-squared test to detect the existence and severity of multicollinearity; (ii) F-test to locate the multicollinearity; (iii) t-test to determine which variable causes multicollinearity.
As the third phase, we performed linear and nonlinear regression analyses to assess the relationships between dependent (i.e., InSAR vertical velocity) and independent variables (i.e., subsidence predisposing factors), starting with a bivariate model and progressively implementing the independent variables. First, we assumed the InSAR Vvert to be linearly correlated to SPF and performed the regression analysis via the equation: From the borehole logs, we extracted the thickness of each stratigraphic unit and calculated the overburden stress imposed on the deeper compressible unit (i.e., unit 3-OS) as: where σ v is the overburden stress (N/m 2 ) at depth z (m); γ is the soil unit weight (N/m 3 ); and T n is the thickness (m) of the n-th stratigraphic unit. As in [37,38,44], the isopiezometric curves representing the mean level of the surficial aquifer in the absence of anthropogenic perturbations were reconstructed from the active piezometers.
We then rasterized the flat study area over 6207 cells of 20 × 20 m and uniformly interpolated both InSAR and stratigraphic point data with cells of the same dimensions in order to create a continuous surface from sample points at sparse locations. Specifically, we used ordinary kriging [50] to interpolate the InSAR vertical velocity (V vert ) and empirical Bayesian kriging [51] to interpolate the thicknesses of the reworked material (unit 1-RM) and organic silts (unit 3-OS), the σ v value, and the piezometric level. Details on the performed interpolations are given in Appendix A. Different interpolators were necessary because of the nonuniform spatial distribution and density of the sample dataset. For each pixel (20 × 20m) of the study area, the thicknesses of unit 1 (RM) and unit 3 (OS), the σ v value, and the piezometric level were summarized and considered as the possible subsidence predisposing factors (SPF).
In the second phase, we performed a bivariate correlation analysis to determine if input variables were linearly dependent using Pearson's correlation coefficient (PCC). Based on the PCC values, the bivariate linear correlation values were classified as null (PCC < |0.25|), weak (|0.25|< PCC < |0.50|), moderately strong (|0.50|< PCC < |0.75|), and strong (|0.75|< PCC < |1|). Specifically, we performed two correlation analyses: the first correlation was between the InSAR vertical velocity (V vert ) and SPF, which suggested the type of model to be used in the third phase; the second correlation was among SPF values, which although being a sufficient but not necessary condition, recommended on the need to check for multicollinearity. Any regression model can indeed be biased by multicollinearity, which can reduce the statistical significance of independent variables (i.e., subsidence predisposing factors) that are instead relevant to explain or predict the dependent variable (i.e., InSAR vertical velocity). To this end, we performed the Farrar-Glauber test (F-G test) [52], which consists of a set of three tests: (i) Chi-squared test to detect the existence and severity of multicollinearity; (ii) F-test to locate the multicollinearity; (iii) t-test to determine which variable causes multicollinearity.
As the third phase, we performed linear and nonlinear regression analyses to assess the relationships between dependent (i.e., InSAR vertical velocity) and independent variables (i.e., subsidence predisposing factors), starting with a bivariate model and progressively implementing the independent variables. First, we assumed the InSAR V vert to be linearly correlated to SPF and performed the regression analysis via the equation: where y(x) is the predicted variable (i.e., InSAR V vert ); β 0 is the parametric component of the linear predictors; β i (i = 1, . . . , n) are the unknown regression coefficients that explain the relationship between dependent and independent variables; x i (i = 1, . . . , n) are the predictors (i.e., subsidence predisposing factors-SPF). Here, we used a generalized linear model (GLM) [53] with Gaussian probability distribution of the predicted variable and an identity link function between the predicted variable and predictors. The GLM is fitted via maximum likelihood estimation. Then, to give more flexibility to the relationship between InSAR V vert and SPF, which could be either linear or nonlinear, we also tested a generalized additive model (GAM) [54] via the equation: where s(·) is a nonparametric smooth function that can capture the nonlinearity between the predicted variable and predictors. To prevent overfitting, the level of smoothness of the s(·) function is regulated by the smoothing parameter, which is estimated with generalized cross-validation criterion; and the GAM, which is fitted via penalized iteratively reweighted least squares. The adjusted R-squared (adj-R 2 ) and Akaike Information Criterion (AIC) values calculated for GLM and GAM allowed us to then assess the best performing model in terms of explaining the relationship between InSAR V vert and SPF. Figure 5 shows the interpolated InSAR vertical velocity (V vert ) and subsidence predisposing factors (i.e., thicknesses of unit 1 and unit 3, overburden stress (σ v ), and the piezometric level). Both ERS 1 and2 (1992-2000) and Envisat (2003-2010) recorded their highest vertical velocities (V vert ) along the lakeshore and in the Cosia stream delta area in the NW sector of the basin (Figure 5a,b).

Results
The lakeshore subsides at a rate of~5 mm/year almost constantly in both analyzed periods; unit 1 is thicker (up to 10 m) in this area than in the entire basin (Figure 5c), and is accompanied by~20 m of unit 3 (Figure 5d), with an overburden stress (σ v ) of 200-300 N/m 2 (Figure 5e). Here, the piezometric level varies from 197 to 199 m a.s.l. (Figure 5f).
The NW sector of the basin shows widespread vertical movements of~10 mm/year in 1992-2000; in 2003-2010 the subsiding area is less spread, although it is affected by the same rate of movement. The Cosia stream delta area covers~3 m of unit 1 (Figure 5c) and~25 m of unit 3 (Figure 5d), with a σ v value of 300 N/m 2 on average ( Figure 5e) and a piezometric level ranging from 197 to 202 m a.s.l. (Figure 5f).

PCC and Type of Correlation
Although the previous analysis suggested a nonlinear relationship between the InSAR V vert and SPF, for completeness we performed both linear and nonlinear analyses via GLM and GAM, respectively. All tested combinations of predictors, from bivariate to multivariate models, are given in Appendix B. Table 4 summarizes the results of the best fitting models, where InSAR V vert is a function of all SPFs (i.e., thicknesses of unit 1 (RM) and unit 3 (OS), σ v value, and piezometric level). The statistical significance of predictors is determined by their p-values (Table 4). As expected, the AIC decreases when passing from GLM to GAM when modelling both ERS 1 and 2 (1992-2000) and Envisat (2003-2010), thus suggesting that more flexible fitting of predictors is necessary for better results.
Finally, to check the effectiveness of these regression models in optimally explaining the relationship between InSAR V vert and predictors (i.e., thicknesses of unit 1 (RM) and unit 3 (OS), σ v value, and piezometric level), we checked the spatial distribution of the regression residuals. These are defined as the differences between observed and predicted values and can be: (i) >0 (i.e., observed > predicted), indicating underestimation; (ii) =0 (i.e., observed = predicted), indicating a perfect fit; (iii) <0 (i.e., observed < predicted), indicating overestimation. Ideally, the residuals should be small and unstructured in order to properly explain the regression model, as is the case for most of the Como basin ( Figure 6). (3) and (4)). The adjusted R-squared (adj-R 2 ) and Akaike Information Criterion (AIC) values are also listed. (iii) <0 (i.e., observed < predicted), indicating overestimation. Ideally, the residuals should be small and unstructured in order to properly explain the regression model, as is the case for most of the Como basin ( Figure 6).

Discussion
The previous section presented the multivariate linear and nonlinear regression analyses performed in the Como basin over two time ranges: 1992-2000 and 2003-2010. The results demonstrate that the relationship between InSAR vertical displacements and subsidence predisposing factors cannot be described using linear models, which may lead to underestimations of the spatial distribution of the hazards. Indeed, land subsidence needs to be modeled with more

Discussion
The previous section presented the multivariate linear and nonlinear regression analyses performed in the Como basin over two time ranges: 1992-2000 and 2003-2010. The results demonstrate that the relationship between InSAR vertical displacements and subsidence predisposing factors cannot be described using linear models, which may lead to underestimations of the spatial distribution of the hazards. Indeed, land subsidence needs to be modeled with more complex nonlinear functions (e.g., hyperbolic, quadratic, seasonal), because it depends on the stage of soil consolidation [11,32]. Rapid deformation events or changes in predisposing or triggering factors (e.g., groundwater fluctuations, urban excavation, or constructions) can generate nonlinear behaviors of the subsoil, thus inducing acceleration or anomalous motion during the monitoring period of a single sensor [11]. Such abrupt events may remain undetected when the InSAR data are unwrapped with a linear deformation model. However, the analysis of PS displacement time series could help in interpreting the subsidence-induced deformations when multitemporal hydrogeological and geotechnical data are available for ground-based comparison. Additionally, a widespread hydromechanical database covering the entire Como basin is needed to improve the regression analysis presented here and encompass the dynamic evolution of groundwater and soil compaction.
In this paper, the subsidence predisposing factors (i.e., thicknesses of unit 1 (RM) and unit 3 (OS), overburden stress, and piezometric level) are assumed invariant between 1992 and 2010 because (i) the study area is completely urbanized and the thicknesses of the sedimentary units cannot be modified; (ii) any anthropic activity (e.g., the construction of new buildings) would be too localized to affect the overburden stress at the basin scale; and (iii) besides the seasonal fluctuations of the piezometric level, its overall trend can be considered constant.
Although we considered the mechanical setting of the Como subsurface as constant over the investigated time ranges, the AIC of nonlinear regression models decreases when moving from 1992-2000 to 2003-2010, thus suggesting that the considered subsidence predisposing factors better explain the ground movements as detected by Envisat (2003-2010) rather than ERS 1 and 2 (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000). From the analysis of regression results, we also noticed a change in the subsidence pattern from one period to the other ( Figure 6). This may be due to (i) improved SAR sensors mounted on Envisat (2003-2010) compared to ERS 1 and 2 (1992-2000); (ii) movements of the InSAR reference point, which may affect the signal response in the whole Como basin; or (iii) unrecorded anomalies caused by unknown natural or anthropic events. Further investigations are needed to properly determine the possible cause of the variation in average subsidence rates from 1992-2000 to 2003-2010; however, a comparison between the two InSAR series is beyond the scope of this paper.
The spatial distribution of regression residuals in Figure 6b,d suggests that the InSAR vertical movements (V vert ) in the NW sector of the basin might be due not only to the thicknesses of unit 1 (RM) and unit 3 (OS), the overburden stress (σ v ), and the piezometric level, but also to additionally factors. Therefore, we explored the possibility that such movements are due to unconsidered natural or anthropic factors. As the NW sector of Como basin is close to the Gonfolite backthrust and crossed by the artificially confined Cosia stream, therefore we introduced two further variables, namely the distance from the backthrust and the distance from the Cosia stream, which might improve our regression models.
Recent investigations [41,42,55,56] recognized a neotectonic activity of the Gonfolite backthrust, leading to a possible influence on the stratigraphic architecture of syn-growth sequences nearby the active range front [57,58]. The Pliocene to Holocene displacements generated an asymmetric geometry in the Como basin filling (Figure 1), with possibly the thickest sequences of compressible Pliocene clays occurring along the SW border of the plain, at the base of the range front carved in the Gonfolite Group conglomerates. A similar structure occurs in the outcropping, for instance a few kilometers west, near Novazzano (Ticino) [41,42,55,56].
On the other side, the current course of the Cosia stream is the result of anthropic modifications that occurred during the Roman age [43,59]. Indeed, before the founding of Como city, both Cosia and Valduce streams wandered throughout the plain, and flooding events were frequent (Figure 7a). In the absence of additional forces, delta aggradation occurred when sediment deposition rates overpassed subsidence rates. The plain was filled with sandy and gravelly deposits originating from the surrounding slopes or transported by the streams. The analysis of stratigraphic data identified natural inlets along the coastline (Figure 7a). When Romans founded the city of Como, the Cosia stream was confined to the west of the town, while the Valduce stream was confined to the opposite side (Figure 7b). This setting is clearly anthropic, because the creek courses are not in the depocenter of the basin but rather run parallel to the SW and NE mountain slopes. From this moment onwards, the overall landscape evolution of Como basin has been dominated by anthropic processes; the sedimentary supply from the creeks decreased due to watershed management. River and lake flooding continuously affected the town and are documented in the stratigraphic records and in literary sources [60,61], but coastal progradation since Roman times has been due to land reclamation and filling of coastal areas. Today, the Cosia and Valduce streams reach the lake after flowing underground beneath the town ( Figure  7b). These artificial courses are located below the topographic surface and locally alter the direction of groundwater flow. The alluvial deposits (unit 2-SG) indeed are coarser than unit 3 (OS) and have a level of higher permeability ( Figure 2).
As the distances from the backthrust and from Cosia stream are spatially confined in the same area, these variables are strongly correlated. A model integrating both variables was not considered because of their collinearity, which was detected through the Farrar-Glauber test. The correlation matrix is given in Appendix B. However, GAM models that introduced the new predictors were improved with respect to the previous analyses, as shown in the histograms in Figure 8. Particularly, the InSAR vertical movements in the NW sector of the basin seem to be related more to the course of the Cosia stream, rather than the distance from the backthrust. Indeed, the performances of both models with ERS 1 and 2 (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) and Envisat (2003-2010) improved slightly (by 20%) when introducing the distance from the backthrust, but the spatial distribution of the regression residuals remained almost unchanged. When introducing the distance from the Cosia stream, the results were enhanced by about 50% with respect to the previous analyses, and the pattern for the regression residuals was locally improved. AIC values and maps of the residuals for these latter regression analyses are given in Appendix B. However, more in-depth analyses are necessary to further investigate the natural and anthropic factors that may influence the InSAR vertical movements in the NW sector of the basin. When Romans founded the city of Como, the Cosia stream was confined to the west of the town, while the Valduce stream was confined to the opposite side (Figure 7b). This setting is clearly anthropic, because the creek courses are not in the depocenter of the basin but rather run parallel to the SW and NE mountain slopes. From this moment onwards, the overall landscape evolution of Como basin has been dominated by anthropic processes; the sedimentary supply from the creeks decreased due to watershed management. River and lake flooding continuously affected the town and are documented in the stratigraphic records and in literary sources [60,61], but coastal progradation since Roman times has been due to land reclamation and filling of coastal areas. Today, the Cosia and Valduce streams reach the lake after flowing underground beneath the town (Figure 7b). These artificial courses are located below the topographic surface and locally alter the direction of groundwater flow. The alluvial deposits (unit 2-SG) indeed are coarser than unit 3 (OS) and have a level of higher permeability ( Figure 2).
As the distances from the backthrust and from Cosia stream are spatially confined in the same area, these variables are strongly correlated. A model integrating both variables was not considered because of their collinearity, which was detected through the Farrar-Glauber test. The correlation matrix is given in Appendix B. However, GAM models that introduced the new predictors were improved with respect to the previous analyses, as shown in the histograms in Figure 8. Particularly, the InSAR vertical movements in the NW sector of the basin seem to be related more to the course of the Cosia stream, rather than the distance from the backthrust. Indeed, the performances of both models with ERS 1 and2 (1992-2000) and Envisat (2003-2010) improved slightly (by 20%) when introducing the distance from the backthrust, but the spatial distribution of the regression residuals remained almost unchanged. When introducing the distance from the Cosia stream, the results were enhanced by about 50% with respect to the previous analyses, and the pattern for the regression residuals was locally improved. AIC values and maps of the residuals for these latter regression analyses are given in Appendix B. However, more in-depth analyses are necessary to further investigate the natural and anthropic factors that may influence the InSAR vertical movements in the NW sector of the basin.

Conclusions
In this paper, we analyzed the spatial distribution of subsidence in the Como basin from 1992 to 2010, benefiting from the deep knowledge of local hydrogeological, geotechnical, and geomorphological properties, and the availability of InSAR measurements. Bivariate and multivariate linear and nonlinear regression analyses (summarized in Appendix B) were performed to investigate the relationship between InSAR-derived vertical movements and subsidence predisposing factors, i.e., thicknesses of compressible soils (named unit 1 (RM) and unit 3 (OS)), overburden stress, and piezometric level, using both a generalized linear model (GLM) and

Conclusions
In this paper, we analyzed the spatial distribution of subsidence in the Como basin from 1992 to 2010, benefiting from the deep knowledge of local hydrogeological, geotechnical, and geomorphological properties, and the availability of InSAR measurements. Bivariate and multivariate linear and nonlinear regression analyses (summarized in Appendix B) were performed to investigate the relationship between InSAR-derived vertical movements and subsidence predisposing factors, i.e., thicknesses of compressible soils (named unit 1 (RM) and unit 3 (OS)), overburden stress, and piezometric level, using both a generalized linear model (GLM) and generalized additive model (GAM). This approach allowed us to determine the nonlinear combination of factors that most likely induce subsidence in the Como basin, thus proving that land subsidence, here measured via InSAR vertical movements, is a nonlinear event subjected to variations in soil compaction and hydromechanical properties.
The best performing GAM model explained 70% of InSAR-derived movements using the above-mentioned subsidence predisposing factors and highlighted some uncertainties in the NW sector of the basin. As possible reasons, we suggested the recent (Late Pliocene to Late glacial) activity of the Gonfolite backthrust [55,56] and the presence of the artificially confined Cosia stream and its delta. We tested our hypotheses by separately including these variables (i.e., distances from the backthrust and Cosia stream) in the regression models. Although the overall performance was improved, the models still showed some uncertainties in the NW sector of the Como basin, especially in years ranging from 1992 to 2000; these results may be due to other geological events that occurred in those years and that were constrained to that sector of the basin, for which no records are available. Further investigation is needed in order to fully verify these hypotheses, which is beyond the scope of the present study. However, the millimetric precision of InSAR data and their availability over the entire Como basin allowed the detection of new factors that may control the ground subsidence in the study area for the first time. This is a novel and remarkable result that may also influence the urban land management and risk assessment in Como Municipality.
Indeed, the possibility to combine local hydromechanical properties with InSAR measurements collected over 18 years leads to a more effective interpretation of the spatial distribution of subsidence in the Como basin. However, further work is needed to refine the statistical explanatory model and move towards subsidence forecasting models using InSAR measurements as alternatives to ground-based technologies [30].
Since the city of Como is affected by the coupling of natural and anthropogenic subsidence that also exposes the historical city center to an increasing risk of lake inundation [3,23,37,38], future improvements may also concern the modelling of the anthropic influence on the local groundwater regime due to the construction of the movable bulkheads on the lakefront. This investigation, although being spatially localized and temporally constrained, might be useful for decoupling natural subsidence causes from anthropic subsidence causes in the Como basin. In this section, details about the interpolation methods are summarized. The vertical velocity values of ERS 1 and 2 (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) and Envisat (2003-2010) were interpolated via ordinary kriging, which was modeled using a spherical semivariogram obtained from [50]: The empirical Bayesian kriging was used to interpolate the thicknesses of the reworked material (unit 1-RM) and organic silts (unit 3-OS), the overburden stress (σ v ), and the piezometric level. These interpolations were modeled using a power semivariogram according to the following equation [51]: Figure A1 shows the semivariograms obtained from the above-described interpolations.
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 21 16 These interpolations were modeled using a power semivariogram according to the following equation [51]: Figure A1 shows the semivariograms obtained from the above-described interpolations.

Appendix B
In this section, all tested regression models are summarized. Table A1 reports the correlations among all considered variables. The distance from the backthrust and the distance from the Cosia stream are strongly correlated; therefore, they induce collinearity if considered in the same model.

Appendix B
In this section, all tested regression models are summarized. Table A1 reports the correlations among all considered variables. The distance from the backthrust and the distance from the Cosia stream are strongly correlated; therefore, they induce collinearity if considered in the same model. Results from the generalized linear model (GLM) and generalized additive model (GAM) are reported in Tables A2 and A3, respectively. Table A2. Summary of all tested linear regression models (M) using generalized linear model (GLM). β 0 is the parametric component of the predictors (see Equations (3) and (4)). The adjusted R-squared (adj-R 2 ) and Akaike Information Criterion (AIC) values are also listed. Significance codes of p-values: *** = 0; ** = 0.001; * = 0.01; . = 0.05.  Table A3. Summary of all tested nonlinear regression models (M) using generalized additive model (GAM). β 0 is the parametric component of the predictors (see Equations (3) and (4)). Here, only the Akaike Information Criterion (AIC) is computed. Significance codes of p-values: *** = 0; . = 0.05.  Figure A2 shows the distribution of residuals obtained from GAM, introducing the distances from the backthrust and from the Cosia stream as further predictors.