Subsidence Evolution of the Firenze – Prato – Pistoia Plain ( Central Italy ) Combining PSI and GNSS Data

Subsidence phenomena, as well as landslides and floods, are one of the main geohazards affecting the Tuscany region (central Italy). The monitoring of related ground deformations plays a key role in their management to avoid problems for buildings and infrastructure. In this scenario, Earth observation offers a better solution in terms of costs and benefits than traditional techniques (e.g., GNSS (Global Navigation Satellite System) or levelling networks), especially for wide area applications. In this work, the subsidence-related ground motions in the Firenze–Prato–Pistoia plain were back-investigated to track the evolution of displacement from 2003 to 2017 by means of multi-interferometric analysis of ENVISAT and Sentinel-1 imagery combined with GNSS data. The resulting vertical deformation velocities are aligned to the European Terrestrial Reference System 89 (ETRS89) datum and can be considered real velocity of displacement. The vertical ground deformation maps derived by ENVISAT and Sentinel-1 data, corrected with the GNSS, show how the area affected by subsidence for the period 2003–2010 and the period 2014–2017 evolved. The differences between the two datasets in terms of the extension and velocity values were analysed and then associated with the geological setting of the basin and external factors, e.g., new greenhouses and nurseries. This analysis allowed for reconstructing the evolution of the subsidence for the area of interest showing an increment of ground deformation in the historic centre of Pistoia Town, a decrement of subsidence in the nursery area between Pistoia and Prato cities, and changes in the industrial sector close to Prato.


Introduction
Land subsidence is commonly defined as a gentle and gradual lowering or sudden sinking of the ground surface [1].It can be caused by natural or anthropogenic processes or by a combination of both effects.Furthermore, it can be due to endogenic or exogenic phenomena, which refer to geological-related motions and the removal of underground materials, respectively [2].Subsidence phenomena mainly affect urban areas [3] because of water overexploitation, with serious consequences such as damage to linear infrastructures, e.g., bridges, roads or railways, and building stability issues due to differential settlement [4][5][6].
Monitoring plays a key role in natural hazard management by providing cost-effective solutions to mitigate or minimize disaster losses.Traditionally, levelling networks [7] or GPS (Global Positioning System) techniques [8][9][10][11] were used to monitor land subsidence.In the last several decades, the SAR (Synthetic Aperture Radar) remote sensing methods have rapidly grown, and, currently, they are a common tool for natural hazard management.Earth Observation (EO) techniques can profitably support risk reduction strategies by taking advantage of their wide area coverage associated with a high cost/benefit ratio.One of these techniques, the Persistent Scatterers Interferometry (PSI) [12], has been successfully adopted for a wide range of applications in disaster management, and it has extensively proven to be a valuable tool to detect ground deformations due to landslides [13][14][15][16][17][18], subsidence [19][20][21][22][23][24][25][26][27][28], volcanoes [29,30] and earthquakes [31].
Subsidence is one of the main geological hazards in Tuscany, as well as landslides [32] and floods [33], but it is usually less studied.By combining the interferometric data from images acquired by the ENVISAT satellite, Rosi et al. [34] mapped the subsidence in the whole Tuscany region (central Italy) by a simple model [35] able to discriminate subsiding areas from other ground deformation.
In this work, the procedure presented by Rosi et al. [34], enhanced by combining SAR and GNSS (Global Navigation Satellite System), has been applied to both ENVISAT (acquired between 2003 and 2010) and Sentinel-1 (acquired between 2014 and 2017) datasets, focusing on the subsidence phenomena affecting the Firenze-Prato-Pistoia (FPP) plain.The radar images were processed according to the PSI (Persistent Scatterers Interferometry) technique in order to investigate the velocity of ground displacement of the Area of Interest (AoI).Successively, the resulting PS data were corrected according to the GNSS data to analyse absolute displacement and to therefore improve the precision of the investigation with respect to the usual multitemporal A-DInSAR (Advanced Interferometric SAR) analysis.The comparison between ENVISAT and Sentinel-1 data allows for investigating the evolution of zones continuously affected by ground lowering and highlighting areas characterized by trend changes.Three cross-sections were traced along the AoI to analyse the subsidence profiles and to investigate the temporal evolution from 2003 to 2017.
The main aim of the work is the monitoring of the spatial and temporal evolution of subsidence in the last 15 years in order to further develop previous works.The area of interest is historically affected by ground deformation due to its geological features and external factors, e.g., new greenhouses and nurseries.In this work, ENVISAT (2003ENVISAT ( -2010) ) and Sentinel-1 (2014-2017) SAR data were combined with GNSS data to correct the relative velocity of displacement in real velocity for monitoring this phenomenon.

Study Area
The area of interest is the Firenze-Prato-Pistoia basin located in the Tuscany region, central Italy (Figure 1).It has an extension of approximately 824 km 2 and a mean elevation of ca.50 m a.s.l.The plain is part of the 8.4% of the flat areas [36] of the Tuscany region, and it is crossed by the Ombrone creek in the Pistoia area, the Bisenzio River in the Prato province and the Arno River in the Florence area.The valley is bordered by the Northern Apennines in the north and in the east, by the Chianti and Senese hillslopes in the south, and by the Valdinievole valley and the Serravalle Pistoiese relief in the west.The plain has an elongated shape in the NW-SE direction and is approximately 35 km wide and 100 km long; three main cities are presented: Florence, Prato and Pistoia.
This plain is an intermontane valley formed during the late extensional phase of the formation of the northern Apennines and later filled with alluvial deposits, reaching a thickness of hundreds of metres.The shape is defined by normal faults in the NO flank forming a semi-graben (late Pliocene) filled by a relevant thickness of marshy-lacustrine sediments [37].The soft sediments are subsiding due to groundwater overexploitation, the tree nursery industry and the load of the urban centre [34,35,38], which has been reported since the 18th and 19th centuries [39,40].
From a geological point of view, the basin is comprised from the bottom to the top [36] (Figure 2) the Macigno formation composed of sandy flysch (Oligocene-Miocene), the Canetolo Complex clayschists (Eocene-Oligocene), Ophiolitic Complex (Lower Jurassic-Lower Cretaceous) and the Calvana Supergroup (Upper Cretaceous-middle Eocene).The stratigraphy is covered by lacustrine, locally fluvial sediments with different thicknesses based on the localization in the AoI and with respect to the boundary of the valley.Generally, the deposits are constituted by pebbles with variable dimensions from centimetres to metres and the presence of gravel, sand and clay.The percentage of these components depends on the distance with respect to the main rivers of the area, e.g., Arno, Ombrone and Bisenzio Rivers.In the past, the AoI has been subjected to some changes in ground movements (e.g., from subsidence to uplift) [32,41] due to the overexploitation of groundwater [41], the recession of the textile industries in the early 2000s and the opening of several greenhouses and nurseries.

Materials and Data
To investigate the subsidence phenomena and evolution between 2003 and 2017 in the AoI, satellite SAR data were combined with the derived information from eight GNSS permanent stations installed in the study area, which has presented continuous and accurate time series over the last several decades.These data have the relevant advantage of recurrent information with a high cost/benefit ratio over wide areas and millimetric precision.In addition, ancillary data such as a 10 m resolution DEM (Digital Elevation Model), derived maps (aspect and slope maps), and the geology of the study area were used to characterize and investigate the vertical ground displacements affecting the area.
The PSI processing is organized in clusters, and, for each cluster, a ground fiducial point assumed to be stable, the reference point, is chosen during the processing phase on the basis of the geological features of the area and the high coherence and reflectivity of the radar images.This means that the recorded PS displacements, derived from a differential technique, are calculated by assuming the stability of the reference points, even if it is not connected to any geodetic reference frame [42].
To by-pass the problem and then to represent the PS velocities in the same geodetic reference frame, GNSS permanent station data were involved to provide the datum alignment.

SAR Data
PSI data from ENVISAT satellite (2003-2010) are made available, in both ascending and descending orbits, from the Portale Cartografico Nazionale (PCN) of the Italian Ministry of Environment and Territory of the Sea (METS) [43] (http://www.pcn.minambiente.it/),while Sentinel-1 radar images can be freely downloaded from Sentinel Scientific Data Hub [44] (https://scihub.copernicus.eu).The Sentinel-1 constellation is composed of Sentinel-1A and Sentinel-1B, which have the same reference orbit with a 180-degree orbital phasing difference, allowing a short revisiting time [32]; in addition, it is a continuation of the ENVISAT mission.Both satellites carry sensors in C-band (wavelength 5.6 cm), so they can be easily compared, even if they have different revisit times: Sentinel satellite has a revisiting period of 12 days from April 2014 and of six days from April 2016 when the twin satellite was launched, while the ENVISAT constellation had a revisiting period of 35 days.Details on the ENVISAT and Sentinel-1 data available and used for the investigation of the FPP plain are reported in Table 1.The ENVISAT data have been processed by the PSP-DIFSAR (Persistent Scatterers Pairs Differential InSAR) [45] and PSInSAR [13] approaches by e-GEOS and TRE-Altamira companies, respectively [34], while the Sentinel-1 images were elaborated by TRE-Altamira company by the SqueeSAR technique [46].Both ENVISAT and Sentinel-1 PSI products cover the entire area with a good point density due to the high amount of reflecting elements, e.g., urban areas, and the advantageous morphology.The difference in the PS number between the ENVISAT and the Sentinel-1 data is due to the algorithm used for the processing of radar images, as described before, and the different coherence thresholds applied in the processing phase (0.6 for the ENVISAT products and 0.7 for the Sentinel-1 data).

GNSS Data
Eight GNSS permanent stations (Figure 1), which have been continuously collecting data for several years, are present inside and close to the study area.Two of them, Florence (IGMI) and Prato (PRAT), belong to the European GNSS Permanent Network (EPN) and to the Italian GNSS Geodetic Reference System named Rete Dinamica Nazionale (RDN).Three of these GNSS stations are located outside of the FPP basin (EMPO, CRMI and LMPR), one in the city centre (FIPR) and one in the periphery (IGMI) of Florence.The other three GNSS stations are located in the town of Calenzano (CALA) near Florence, in the city of Prato (PRAT) and in Seano village (SEAN), on the east side of Prato city (Table 2).The data availability of the GNSS stations start in 1998 for PRAT station, 2006 for IGMI and CRMI, and 2010 for the remaining sites.Note that the GNSS and PSI datasets overlap for more than seven years, and the hypothesis of linear motions has been assumed to extrapolate the GNSS velocities.The reason is that all the GNSS sites located in study area present velocities determined by GNSS and PSI techniques that are characterized by a linear and almost constant pattern across time in horizontal and vertical components but with different rates.The GNSS data processing was performed using Bernese GPS Software (Version 5.0.2007, Astronomical Institute of University of Bern (AIUB); Bern, Switzerland) following the strategy processing of the International Earth Rotation and Reference Systems Service Conventions 2010 and EUREF (European Reference Frame) guideline.Daily double differenced GNSS phase observations were used to estimate station coordinates and the velocity uncertainties were estimated with the combination of flicker and white noise.Daily free normal equation solutions were stacked to produce a single weekly normal equation.Time series weekly positioning solutions were analysed for determining the velocity by checking for the possible outliers and estimating the annual and semi-annual seasonal effects, correcting the eventual discontinuities and aligning with the ETRS89 reference frame [47,48].

Methodology
The subsidence mapping by ENVISAT and Sentinel-1 data is based on two steps.The first step aims to correct the SAR velocities since they are calculated with respect to a set of unknown reference points defined during the processing phase as motionless stable points according to the geological data and high coherence.This hypothesis that reference points are motionless causes incorrect determination of PS velocities in the case that reference points are not stable.In detail, the displacements and velocities determined with SAR interferometry can represent only relative motions, which is, in the case of PS Interferometry, the relative motion between the PS points.In order to obtain absolute measurements, it is necessary to determine displacements and velocities with respect to a geodetic datum (ETRS89), a step that was performed by using the GNSS technique.In particular, a set of GNSS permanent stations located in the study area were involved in the calibration of SAR images.GNSS permanent stations collect data every 30 s, 24 h a day, seven days a week and these conditions permit determining their positions with millimetric accuracy.The determination of the horizontal and vertical components of velocity of each GNSS site is calculated by the analysis of time series positions after removing periodic fluctuations and eventual discontinuities.
To compare PS displacements and velocities with GNSS, SAR measurement vectors that are recorded along the satellite LOS (Line Of Sight) direction must be projected into their horizontal and normal components with respect to the international ellipsoid (GR80) such as the GNSS dataset.The velocity v LOS measured by the SAR satellite for a ground movement v = [v N , v E , v V ] can be represented by the direction cosine S = [S N , S E , S V ] that identifies the unit vector in the direction from the ground to the satellite (LOS direction) (Figure 3): (1)

Methodology
The subsidence mapping by ENVISAT and Sentinel-1 data is based on two steps.The first step aims to correct the SAR velocities since they are calculated with respect to a set of unknown reference points defined during the processing phase as motionless stable points according to the geological data and high coherence.This hypothesis that reference points are motionless causes incorrect determination of PS velocities in the case that reference points are not stable.In detail, the displacements and velocities determined with SAR interferometry can represent only relative motions, which is, in the case of PS Interferometry, the relative motion between the PS points.In order to obtain absolute measurements, it is necessary to determine displacements and velocities with respect to a geodetic datum (ETRS89), a step that was performed by using the GNSS technique.In particular, a set of GNSS permanent stations located in the study area were involved in the calibration of SAR images.GNSS permanent stations collect data every 30 s, 24 h a day, seven days a week and these conditions permit determining their positions with millimetric accuracy.The determination of the horizontal and vertical components of velocity of each GNSS site is calculated by the analysis of time series positions after removing periodic fluctuations and eventual discontinuities.
To compare PS displacements and velocities with GNSS, SAR measurement vectors that are recorded along the satellite LOS (Line Of Sight) direction must be projected into their horizontal and normal components with respect to the international ellipsoid (GR80) such as the GNSS dataset.The velocity measured by the SAR satellite for a ground movement = [ , , ] can be represented by the direction cosine = [ , , ] that identifies the unit vector in the direction from the ground to the satellite (LOS direction) (Figure 3): Writing Equation (1) for the ascending (VA) and descending (VD) velocities by using the direction cosine gives the following: Equation ( 3) is a system of two equations with three variables and is thus not solvable.Displacements occurring along the north-south direction, almost parallel to the satellite orbit, cannot be measured accurately, as their projection along the LOS is negligible for both the ascending and The component of unit vector S for ascending (S A ) and descending (S D ) orbits are derived by the inclination of the orbits with respect to the equator, δ, and the mean off-nadir, θ: Writing Equation (1) for the ascending (V A ) and descending (V D ) velocities by using the direction cosine gives the following: Equation ( 3) is a system of two equations with three variables and is thus not solvable.Displacements occurring along the north-south direction, almost parallel to the satellite orbit, cannot be measured accurately, as their projection along the LOS is negligible for both the ascending and descending orbits.Assuming that the projection of the north component of the velocity along the LOS is negligible in both satellites, Equation (3) can be approximated as: where the two components of velocity v can be derived as follows: The parameters of ENVISAT and Sentinel-1 constellations, i.e., the orbit inclinations with respect to the equator and the mean off-nadir angle and the modules of the components, are presented in Table 3.The GNSS velocities derived by the permanent station located in the study area were used for the data alignment of the PS velocities to the ETRS89 reference frame.The correction value to be applied to the PS dataset is determined in each GNSS site by the difference of the velocity measured with the GNSS and PS techniques.The difference of the velocity derived by these techniques for each GNSS site consists of a set of sparse correction velocities.According to the ground deformation maps made by the vertical velocity corrected by the GNSS data and the approach described by Crisci et al. [33] to identify the areas affected by vertical displacement, the available PS velocities were interpolated by Inverse Distance Weight (IDW) [49,50] applying a weight calculated as the inverse of the square of the Euclidean distance.After the correction, the resulting PS velocities were aligned to the European geodetic reference frame ETRS89.These absolute (not the mathematical meaning) values of vertical velocity were then used to map the subsidence in the AoI from 2003 to 2017.
To calculate the subsidence rate, the corrected data from both ENVISAT and Sentinel-1 were analysed using the procedure described by Rosi et al. [35].This procedure was developed to map the subsidence at regional scale, taking into account the variations in the morphology of a certain study area.This procedure, starting from the vertical and horizontal components of displacements measured by the satellite, is based on the calculation of the direction of the displacement; then, the horizontal and vertical displacements are divided, and the ground movements that are mainly horizontal are filtered out since subsidence (or uplift) can be considered as mainly vertical movements (Figure 4).

Results
The analysis of the GNSS data demonstrated that the areas considered as stable are instead affected by the continuous displacement of the surface.In the study area, the GNSS time series exhibit fluctuant and cyclical ground movement by means of a millimetric accuracy and daily values that generally follow a linear velocity trend.The time series of the vertical component of PS data around each GNSS station do not show any cyclical short-term oscillations (Figure 5) but confirm the general trend of deformation, which is comparable with the outcomes of two recent studies [51,52] about the geodynamics of the Italian peninsula.
The absolute values of the vertical ground displacement of the FPP basin were assessed for both the ENVISAT and Sentinel-1 SAR data.The differences in the vertical displacements detected by the ENVISAT satellites before and after the GNSS correction range between +1 mm/year in the southwest boundary of the AoI and −2 mm/year close to the central part of the AoI (Figure 6).The southeastern portion of the plain has a

Results
The analysis of the GNSS data demonstrated that the areas considered as stable are instead affected by the continuous displacement of the surface.In the study area, the GNSS time series exhibit fluctuant and cyclical ground movement by means of a millimetric accuracy and daily values that generally follow a linear velocity trend.The time series of the vertical component of PS data around each GNSS station do not show any cyclical short-term oscillations (Figure 5) but confirm the general trend of deformation, which is comparable with the outcomes of two recent studies [51,52] about the geodynamics of the Italian peninsula.
The absolute values of the vertical ground displacement of the FPP basin were assessed for both the ENVISAT and Sentinel-1 SAR data.

Results
The analysis of the GNSS data demonstrated that the areas considered as stable are instead affected by the continuous displacement of the surface.In the study area, the GNSS time series exhibit fluctuant and cyclical ground movement by means of a millimetric accuracy and daily values that generally follow a linear velocity trend.The time series of the vertical component of PS data around each GNSS station do not show any cyclical short-term oscillations (Figure 5) but confirm the general trend of deformation, which is comparable with the outcomes of two recent studies [51,52] about the geodynamics of the Italian peninsula.
The absolute values of the vertical ground displacement of the FPP basin were assessed for both the ENVISAT and Sentinel-1 SAR data.The differences in the vertical displacements detected by the ENVISAT satellites before and after the GNSS correction range between +1 mm/year in the southwest boundary of the AoI and −2 mm/year close to the central part of the AoI (Figure 6).The southeastern portion of the plain has a good correlation with the geodetic data, and the correction is close to zero, whereas the correlation is higher in the southwestern sector of the AoI with values of approximately −2 mm/year.The corrected The differences in the vertical displacements detected by the ENVISAT satellites before and after the GNSS correction range between +1 mm/year in the southwest boundary of the AoI and −2 mm/year close to the central part of the AoI (Figure 6).The southeastern portion of the plain has a good correlation with the geodetic data, and the correction is close to zero, whereas the correlation is higher in the southwestern sector of the AoI with values of approximately −2 mm/year.The corrected vertical velocities show relevant ground subsidence in the southern area of Pistoia Province, with values higher than −20 mm/year.The involved area is approximatively 80 km 2 wide and has a NW-SE elongated shape.In addition, some localized subsidence areas can be identified close to Campi Bisenzio town, with a peak of −20 mm/year of subsidence, and at two zones in the southern portion of Prato city and in the industrial area between Calenzano and Sesto Fiorentino towns, with displacement velocities of approximately −13 mm/year and −10 mm/year, respectively.Except for these two localized cases, the Firenze Province results in being stable and without evident vertical ground deformations from 2003 to 2010.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 19 vertical velocities show relevant ground subsidence in the southern area of Pistoia Province, with values higher than −20 mm/year.The involved area is approximatively 80 km 2 wide and has a NW-SE elongated shape.In addition, some localized subsidence areas can be identified close to Campi Bisenzio town, with a peak of −20 mm/year of subsidence, and at two zones in the southern portion of Prato city and in the industrial area between Calenzano and Sesto Fiorentino towns, with displacement velocities of approximately −13 mm/year and −10 mm/year, respectively.Except for these two localized cases, the Firenze Province results in being stable and without evident vertical ground deformations from 2003 to 2010.The differences between the Sentinel data before and after the GNSS correction are lower than those for the ENVISAT data, ranging from less than 0.5 mm/year in the southeastern end of the AoI, i.e., Firenze Province, to approximately 1.0 mm/year in Pistoia Province, with a peak of 1.5 mm/year close to Prato Province (Figure 7).The period monitored by the Sentinel-1 constellation shows a general reduction in the ground movements in the AoI.The two areas with higher ground deformations (Pistoia and Campi Bisenzio) in the period monitored with the ENVISAT satellites show different vertical displacements from 2014 to 2017.In the Pistoia Province, a smaller subsiding area is recognizable in the south end of the city, and the elongated shape extending in the SE direction disappeared, i.e., the area subsided over time by approximatively 50 km 2 .On the other hand, two new subsidence areas are recognizable with respect to the previous dataset: (i) the historical centre of Pistoia, characterized by a ground subsidence of ca. 13 mm/year and (ii) the SSW part of Pistoia, a recently urbanized territory where the ground subsidence is higher than 20 mm/year.Moreover, a small region of uplift, with values of approximately 5 mm/year, is recognizable in the northwest part of the AoI between the cities of Prato and Firenze.
The area of Campi Bisenzio town and the region between Campi Bisenzio and Calenzano towns The differences between the Sentinel data before and after the GNSS correction are lower than those for the ENVISAT data, ranging from less than 0.5 mm/year in the southeastern end of the AoI, i.e., Firenze Province, to approximately 1.0 mm/year in Pistoia Province, with a peak of 1.5 mm/year close to Prato Province (Figure 7).The period monitored by the Sentinel-1 constellation shows a general reduction in the ground movements in the AoI.The two areas with higher ground deformations (Pistoia and Campi Bisenzio) in the period monitored with the ENVISAT satellites show different vertical displacements from 2014 to 2017.In the Pistoia Province, a smaller subsiding area is recognizable in the south end of the city, and the elongated shape extending in the SE direction disappeared, i.e., the area subsided over time by approximatively 50 km 2 .On the other hand, two new subsidence areas are recognizable with respect to the previous dataset: (i) the historical centre of Pistoia, characterized by a ground subsidence of ca. 13 mm/year and (ii) the SSW part of Pistoia, a recently urbanized territory where the ground subsidence is higher than 20 mm/year.Moreover, a small region of uplift, with values of approximately 5 mm/year, is recognizable in the northwest part of the AoI between the cities of Prato and Firenze.
The area of Campi Bisenzio town and the region between Campi Bisenzio and Calenzano towns show values of approximately −5.8 mm/year and −3.5 mm/year, respectively, which is lower than the previous deformation map made from the ENVISAT data.The area of Sesto Fiorentino exhibits an uplift with values of approximately 4 mm/year, contrary to the subsidence recorded in the same area during the period of 2003-2010 by the ENVISAT constellation.
The stability of the southern portion of the AoI (Firenze Province) is confirmed since areas with significant uplift or subsidence were not detected; these results for Firenze are in agreement with the previous deformation map covered by the ENVISAT data.Two time series for the main subsiding area (Pistoia area), the first time series in the Pistoia centre and the second time series in the southern part of the town, were realized combining the ENVISAT (2003-2010) and the Sentinel-1 (2014-2017) information in order to visualize the entire subsidence evolution of the two areas.To stich the two time series, the value of the cumulative displacement of the first acquisition of Sentinel-1 was predicted on the basis of the ENVISAT data assuming a linear trend of deformation.
The cumulative time series shows different evolutions for the two sites (Figure 8a), as already recognized by the visual comparison of the deformation maps from ENVISAT and Sentinel-1.The historic centre of Pistoia shows a stable trend, even if partially affected by signal noise, until December 2007 in the period monitored by the ENVISAT constellation, while, successively until July 2010, the area slowly subsided.An important increase in the subsidence rate was recorded starting at the end of 2014 by the Sentinel-1 sensors.
The southern area of Pistoia Province shows, instead, an opposite behaviour.In the period monitored by the ENVISAT constellation, an important and constant trend of subsidence is recognizable, while the displacement recorded by the Sentinel-1 satellites show a lower rate of subsidence, with a relatively stable trend in the last period of 2017 (Figure 8b).Two time series for the main subsiding area (Pistoia area), the first time series in the Pistoia centre and the second time series in the southern part of the town, were realized combining the ENVISAT (2003-2010) and the Sentinel-1 (2014-2017) information in order to visualize the entire subsidence evolution of the two areas.To stich the two time series, the value of the cumulative displacement of the first acquisition of Sentinel-1 was predicted on the basis of the ENVISAT data assuming a linear trend of deformation.
The cumulative time series shows different evolutions for the two sites (Figure 8a), as already recognized by the visual comparison of the deformation maps from ENVISAT and Sentinel-1.The historic centre of Pistoia shows a stable trend, even if partially affected by signal noise, until December 2007 in the period monitored by the ENVISAT constellation, while, successively until July 2010, the area slowly subsided.An important increase in the subsidence rate was recorded starting at the end of 2014 by the Sentinel-1 sensors.
The southern area of Pistoia Province shows, instead, an opposite behaviour.In the period monitored by the ENVISAT constellation, an important and constant trend of subsidence is recognizable, while the displacement recorded by the Sentinel-1 satellites show a lower rate of subsidence, with a relatively stable trend in the last period of 2017 (Figure 8b).A successive step was the creation of a contour map of the areas affected by the subsidence in the periods 2003-2010 and 2014-2017 in order to analyse the evolution of the subsidence bowls.According to the ground deformation maps made by the vertical velocity corrected by the GNSS data and the approach described by Crisci et al. [33] to identify areas affected by vertical displacement, the available PS data were interpolated by means of the inverse IDW method [49].The result is a raster map from which it is possible to extract displacement isolines for both the ENVISAT (Figure 9a) and Sentinel-1 data (Figure 9b).
To better highlight the spatial subsidence evolution, three sections, one longitudinal (NW-SE) and two transversals (NE-SW), were traced for both Pistoia (Figure 10a) and Prato (Figure 11a) cities.
The comparison of the profiles extracted for the ENVISAT (orange lines) and the Sentinel-1 (blue lines) confirm the interpretations made by the ground deformations and by the time series.In profile A-A' (Figure 10b), it is evident that the southern area of Pistoia is affected by the ground subsidence, mainly in the period 2003-2010, while the centre of Pistoia (2000 m in the x-axis of Figure 10b) in the period 2014-2017 was affected by an important subsidence event that was not detected before.The same phenomenon is identifiable on the left side (corresponding to the northwestern part of Pistoia) of transversal section B-B' (Figure 10c), while the right side adequately represents the subsidence affecting the recently built district.The third profile C-C' (Figure 10d) shows that the evolution of the subsidence in the southern area close to Pistoia is still substantial.
The same analysis was conducted for Prato city (Figure 11a).A successive step was the creation of a contour map of the areas affected by the subsidence in the periods 2003-2010 and 2014-2017 in order to analyse the evolution of the subsidence bowls.According to the ground deformation maps made by the vertical velocity corrected by the GNSS data and the approach described by Crisci et al. [33] to identify areas affected by vertical displacement, the available PS data were interpolated by means of the inverse IDW method [49].The result is a raster map from which it is possible to extract displacement isolines for both the ENVISAT (Figure 9a) and Sentinel-1 data (Figure 9b).
To better highlight the spatial subsidence evolution, three sections, one longitudinal (NW-SE) and two transversals (NE-SW), were traced for both Pistoia (Figure 10a) and Prato (Figure 11a) cities.
The comparison of the profiles extracted for the ENVISAT (orange lines) and the Sentinel-1 (blue lines) confirm the interpretations made by the ground deformations and by the time series.In profile A-A' (Figure 10b), it is evident that the southern area of Pistoia is affected by the ground subsidence, mainly in the period 2003-2010, while the centre of Pistoia (2000 m in the x-axis of Figure 10b) in the period 2014-2017 was affected by an important subsidence event that was not detected before.The same phenomenon is identifiable on the left side (corresponding to the northwestern part of Pistoia) of transversal section B-B' (Figure 10c), while the right side adequately represents the subsidence affecting the recently built district.The third profile C-C' (Figure 10d) shows that the evolution of the subsidence in the southern area close to Pistoia is still substantial.
The same analysis was conducted for Prato city (Figure 11a).Profile D-D' (Figure 11b) shows that the subsidence affecting the southeastern region of Prato increased in the recent period monitored by the Sentinel-1 constellation (blue line) with respect to the previous period monitored by the ENVISAT satellite (orange line).An apparent inversion of the displacement trend is shown in profile E-E' (Figure 11c), but the velocity values are within one standard deviation of the data (±2.0 mm/year), which is considered the stability range.Profile F-F' (Figure 11d) confirms the detected reduction in the involved area affected by ground subsidence.The same decrease in the subsidence during the recent period of 2014-2017 was also recorded by the deformation maps and highlighted by the contours of the subsidence areas during 2003-2010 in the south of Prato, the industrial district between Sesto Fiorentino and Calenzano towns.Profile D-D' (Figure 11b) shows that the subsidence affecting the southeastern region of Prato increased in the recent period monitored by the Sentinel-1 constellation (blue line) with respect to the previous period monitored by the ENVISAT satellite (orange line).An apparent inversion of the displacement trend is shown in profile E-E' (Figure 11c), but the velocity values are within one standard deviation of the data (±2.0 mm/year), which is considered the stability range.Profile F-F' (Figure 11d) confirms the detected reduction in the involved area affected by ground subsidence.The same decrease in the subsidence during the recent period of 2014-2017 was also recorded by the deformation maps and highlighted by the contours of the subsidence areas during 2003-2010 in the south of Prato, the industrial district between Sesto Fiorentino and Calenzano towns.

Discussion
The analysis and monitoring of the subsidence phenomena affecting the Firenze-Prato-Pistoia plain is important to control the evolution of this natural hazard to avoid the occurrences of issues for buildings and infrastructure.In this work, ancillary, geodetic and SAR data, i.e., ENVISAT (2003ENVISAT ( -2010) ) and Sentinel-1 (2014-2017) were combined and compared to monitor the evolution of the subsidence in the AoI from 2003 to 2017.
The PS data velocities are relative measurements reported to reference points assumed stable according to geological and signal analysis.To have absolute values of the velocity displacement, the PS data for both the ENVISAT and Sentinel-1 constellations were projected into vertical and

Discussion
The analysis and monitoring of the subsidence phenomena affecting the Firenze-Prato-Pistoia plain is important to control the evolution of this natural hazard to avoid the occurrences of issues for buildings and infrastructure.In this work, ancillary, geodetic and SAR data, i.e., ENVISAT (2003ENVISAT ( -2010) ) and Sentinel-1 (2014-2017) were combined and compared to monitor the evolution of the subsidence in the AoI from 2003 to 2017.
The PS data velocities are relative measurements reported to reference points assumed stable according to geological and signal analysis.To have absolute values of the velocity displacement, the PS data for both the ENVISAT and Sentinel-1 constellations were projected into vertical and horizontal components, i.e., east-west and zenith-nadir, and then the vertical component was corrected by the GNSS network to estimate the subsidence rate in the AoI, which is historically known for subsidence phenomena.
This approach allowed the identification of different areas affected by subsidence during both monitored periods: 2003-2010 with ENVISAT and 2014-2017 by Sentinel-1 satellites.The comparison of these data was used to investigate the evolution of the subsidence in the AoI according to the vertical absolute displacement and to trace the limits of the subsiding areas for both investigated periods.
Pistoia town and the surrounding area of Prato city show the most relevant ground deformation evolutions, while the Florence Province results are stable in both datasets, except for in three localized zones.
During the period monitored by the ENVISAT satellite, an important subsidence event affecting the southern region of Pistoia Province was identified, with an elongated shape in the SE direction (Figure 9a).In addition, three subsiding areas were identified, i.e., Campi Bisenzio town, the region between Campi Bisenzio and Calenzano towns, and Sesto Fiorentino city.During the period monitored by the Sentinel-1 satellites, the area of Pistoia town affected by the lowering of the ground had lower subsidence rates than before, but two subsidence bowls were identified for the first time: one in the historical city centre and the second in a recently built-up area to the northwest with strategic infrastructure, i.e., the hospital (Figure 8b).Furthermore, some localized subsidence events were identified in the Eastern Pistoia Province, even if with lower values and smaller extensions.
A comparison between the subsidence profile (Figure 10b) and a geological cross section (Figure 12) [24], traced close section A-A', allowed the proposal of a theory about the subsidence causes and evolution.The main subsidence area is located ca. in the central part of the sedimentary valley (between Brama stream to the NW and Ombrone River to the SE in Figure 12), where very thick sediments occur.It is also possible to notice that, between 2003 and 2010, the subsidence gradually decreased from SE to NW in the same direction that the sediment thickness decreased.From 2014, however, it is possible to notice a subsiding area in the historical centre of Pistoia with velocities up to ca. 20 mm/year; in the cross section in Figure 12, this subsiding area is located between the Ombrone River to the NW and Brama stream to the SE, where the sediment layers are thinner.
Furthermore, it is possible to highlight the presence of a sort of step corresponding to Brama stream that separates the main valley to the SE and a secondary valley to the NW. Figure 10b shows that the subsidence in the SE part of Pistoia has recently decreased, probably because of the diffusion of several industries in this area.
At the beginning of the 21st century, several greenhouses and nurseries were opened, causing substantial subsidence because of the overexploitation of groundwater.The high request of this resource from these companies could have caused a rapid reduction in the groundwater level, with the relevant consequence of ground subsidence, which was probably due to the compaction of the soft layers in the sedimentary deposits.Furthermore, the initial widening of the cone of depression around the exploitation area could lead to a progressive expansion of the area affected by sediment compaction, which can explain the linearly decreasing subsidence rate from the SE (where well fields are located) towards the NW.In the following years (2014-2017), the area was influenced by the decreasing subsidence, probably due to the stabilization of the ground water level and the consequent reduction in the compaction rate of the soft layers of sediment.
The cause of the recent evolution of subsidence in the historical centre of Pistoia could be ascribed to the ongoing widening of the depression cone around the groundwater exploitation area in the SE part of Pistoia.This widening could be caused by the reduction in water pressure in the alluvial/lacustrine deposits of the NW basin, with a consequent compaction and settlement of the soft layers.
The zone characterized by a lower rate of subsidence between the two main subsiding areas corresponds to the step under the Brama stream.Here, the thickness of the deposits is lower than that in the other two valleys, while their granulometry is generally higher; this could lead to a lower compaction of the sediments and subsequently to a lower subsidence rate.It is obvious that this is only a theoretical scenario, since direct groundwater level measurements or characterizations of the deposits are not available.
only a theoretical scenario, since direct groundwater level measurements or characterizations of the deposits are not available.The area surrounding Prato shows a reduction in the subsidence due to the continuous decrease in the textile industries in this area.Several companies were moved to other countries to reduce the manpower costs; this means lower groundwater exploitation and, therefore, a rise in the groundwater level.

Conclusions
The evolution of ground deformation in the Firenze-Prato-Pistoia basin, historically known for subsidence phenomena, was investigated for the periods 2003-2010 and 2014-2017.The analysis was conducted by means of a combined use of ENVISAT and Sentinel-1 satellite data and GNSS data to investigate the absolute velocity values of displacements affecting the AoI.In the period monitored with the ENVISAT satellite, an area close to Pistoia of approximately 80 km 2 that was affected by subsidence was identified, with velocities higher than 20 mm/year.In addition, three more subsiding areas in Prato, Campi Bisenzio and in the territory between Calenzano and Sesto Fiorentino were identified and mapped.The Sentinel-1 data showed that the subsiding area close to Pistoia is reducing, from 80 km 2 to ca.50 km 2 , but with similar subsidence rates.The work allowed for highlighting a general reduction of both subsidence rates and extension in the whole plain, but two new areas affected by ground deformation after 2014 were identified in the NW part of the plain: (i) the historical centre of Pistoia and (ii) the recently urbanized area in the SSW end of Pistoia.In the same period, the other areas showed a lower subsidence rate or uplift than that in the Sesto Fiorentino area.As a result of the investigation of the ground deformation over the whole plain, the isolines of the subsidence rates for both monitoring periods were created.The investigation of the differences between the ENVISAT and Sentinel-1 data allowed the evaluation of the evolution of the deformation in the AoI; for this purpose, several sections in the main subsiding areas were traced to analyse the subsidence profiles.These sections were later compared with a geological section to define a plausible hypothesis about the causes and the geological mechanism of the subsidence of Pistoia.The area surrounding Prato shows a reduction in the subsidence due to the continuous decrease in the textile industries in this area.Several companies were moved to other countries to reduce the manpower costs; this means lower groundwater exploitation and, therefore, a rise in the groundwater level.

Conclusions
The evolution of ground deformation in the Firenze-Prato-Pistoia basin, historically known for subsidence phenomena, was investigated for the periods 2003-2010 and 2014-2017.The analysis was conducted by means of a combined use of ENVISAT and Sentinel-1 satellite data and GNSS data to investigate the absolute velocity values of displacements affecting the AoI.In the period monitored with the ENVISAT satellite, an area close to Pistoia of approximately 80 km 2 that was affected by subsidence was identified, with velocities higher than 20 mm/year.In addition, three more subsiding areas in Prato, Campi Bisenzio and in the territory between Calenzano and Sesto Fiorentino were identified and mapped.The Sentinel-1 data showed that the subsiding area close to Pistoia is reducing, from 80 km 2 to ca.50 km 2 , but with similar subsidence rates.The work allowed for highlighting a general reduction of both subsidence rates and extension in the whole plain, but two new areas affected by ground deformation after 2014 were identified in the NW part of the plain: (i) the historical centre of Pistoia and (ii) the recently urbanized area in the SSW end of Pistoia.In the same period, the other areas showed a lower subsidence rate or uplift than that in the Sesto Fiorentino area.As a result of the investigation of the ground deformation over the whole plain, the isolines of the subsidence rates for both monitoring periods were created.The investigation of the differences between the ENVISAT and Sentinel-1 data allowed the evaluation of the evolution of the deformation in the AoI; for this purpose, several sections in the main subsiding areas were traced to analyse the subsidence profiles.These sections were later compared with a geological section to define a plausible hypothesis about the causes and the geological mechanism of the subsidence of Pistoia.

19 Figure 1 .
Figure 1.Locations of the Area of Interest (AoI) and the eight GNSS stations used for the geodetic correction: five stations are inside the AoI and three are outside the AoI but close to the FPP plain.

Figure 2 .
Figure 2. Sample of a stratigraphic column of the AoI (modified from[36]).

Figure 1 . 19 Figure 1 .
Figure 1.Locations of the Area of Interest (AoI) and the eight GNSS stations used for the geodetic correction: five stations are inside the AoI and three are outside the AoI but close to the FPP plain.

Figure 2 .
Figure 2. Sample of a stratigraphic column of the AoI (modified from[36]).

Figure 2 .
Figure 2. Sample of a stratigraphic column of the AoI (modified from[36]).

Figure 3 .
Figure 3. Representation of the cosine directors in both ascending and descending geometries for decomposing the LOS velocity in Vertical and Horizontal components.

Figure 3 .
Figure 3. Representation of the cosine directors in both ascending and descending geometries for decomposing the LOS velocity in Vertical and Horizontal components.

Figure 5 .
Figure 5.Comparison of the GNSS data recorded by the IGMI station and the nearby ENVISAT and Sentinel-1 PS data.

Figure 5 .
Figure 5.Comparison of the GNSS data recorded by the IGMI station and the nearby ENVISAT and Sentinel-1 PS data.

Figure 5 .
Figure 5.Comparison of the GNSS data recorded by the IGMI station and the nearby ENVISAT and Sentinel-1 PS data.

Figure 6 .
Figure 6.Corrected ENVISAT data (vertical velocity) and isolines of the correction values adopted for the vertical velocities.

Figure 6 .
Figure 6.Corrected ENVISAT data (vertical velocity) and isolines of the correction values adopted for the vertical velocities.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 19 previous deformation map made from the ENVISAT data.The area of Sesto Fiorentino exhibits an uplift with values of approximately 4 mm/year, contrary to the subsidence recorded in the same area during the period of 2003-2010 by the ENVISAT constellation.The stability of the southern portion of the AoI (Firenze Province) is confirmed since areas with significant uplift or subsidence were not detected; these results for Firenze are in agreement with the previous deformation map covered by the ENVISAT data.

Figure 7 .
Figure 7. Corrected Sentinel-1 data (vertical velocity) and isolines of the correction values adopted for the vertical velocities.

Figure 7 .
Figure 7. Corrected Sentinel-1 data (vertical velocity) and isolines of the correction values adopted for the vertical velocities.

19 Figure 8 .
Figure 8.Time series of the two PS sites in the two areas affected by subsidence in Pistoia Province: city centre of Pistoia (a) and southern area of Pistoia (b).

Figure 8 .
Figure 8.Time series of the two PS sites in the two areas affected by subsidence in Pistoia Province: city centre of Pistoia (a) and southern area of Pistoia (b).

19 Figure 9 .
Figure 9. Boundary of the subsidence areas recorded by the ENVISAT (2003-2010, (a)) and Sentinel-1 (2014-2017, (b)) constellations.Different subsidence rates between the two monitored periods may be due to changes in groundwater exploitation rates.

Figure 9 .
Figure 9. Boundary of the subsidence areas recorded by the ENVISAT (2003-2010, (a)) and Sentinel-1 (2014-2017, (b)) constellations.Different subsidence rates between the two monitored periods may be due to changes in groundwater exploitation rates.

Figure 12 .
Figure 12.Geological section of Pistoia town (modified from [37]) across to nearly section A-A' in Figure 10.

Figure 12 .
Figure 12.Geological section of Pistoia town (modified from [37]) across to nearly section A-A' in Figure 10.

Table 1 .
Details of the available datasets of ENVISAT and Sentinel-1A and Sentinel-1B used to investigate the ground deformation in the Firenze-Prato-Pistoia plain.

Table 3 .
Parameters of the Sentinel and ENVISAT images used for this work.