Monitoring of Land Subsidence in the Po River Delta (Northern Italy) Using Geodetic Networks

: The Po River Delta (PRD, Northern Italy) has been historically affected by land subsidence due to natural processes and human activities, with strong impacts on the stability of the natural ecosystems and signiﬁcant socio-economic consequences. This paper is aimed to highlight the spatial and temporal evolution of the land subsidence in the PRD area analyzing the geodetic observations acquired in the last decade. The analysis performed using a moving window approach on Continuous Global Navigation Satellite System (CGNSS) time-series indicates that the velocities, in the order of 6 mm/year, are not affected by signiﬁcant changes in the analyzed period. Furthermore, the use of non-permanent sites belonging to a new GNSS network (measured in 2016 and 2018) integrated with InSAR data (from 2014 to 2017) allowed us to improve the spatial coverage of data points in the PRD area. The results suggest that the land subsidence velocities in the easternmost part of the area of interest are characterized by values greater than the ones located in the western sectors. In particular, the sites located on the sandy beach ridge in the western sector of the study area are characterized by values greater than − 5 mm/year, while rates of about − 10 mm/year or lower have been observed at the eastern sites located in the Po river mouths. The morphological analysis indicates that the land subsidence observed in the PRD area is mainly due to the compaction of the shallow layers characterized by organic-rich clay and fresh-water peat.


Introduction
River deltas, which host large population and extensive economic activities, are among the territories most vulnerable to land subsidence [1,2]: the effects of this phenomenon are linked to environmental degradation, morphological changes of coastlines and emerged surfaces, damage to buildings, and interruption of services [3][4][5].
Land subsidence afflicts many areas of the world, in particular the ones located along transitional environments, such as coastal areas, deltas, wetlands, and lagoons, which are becoming increasingly vulnerable to flooding, storm surges, salinization, and permanent inundation [6][7][8][9]. In these areas, subsidence can be usually considered as a consequence of a complex combination of natural and anthropogenic factors: the compaction of Holocene sediments, tectonic movements, sinkholes formation, volcanism, thawing permafrost, and the Glacial Isostatic Adjustment (GIA), are generally considered as the main natural sources of land subsidence [10][11][12]; aquifer-system compaction associated with groundwater/oil/natural gas depletion and storage, drainage of organic soils, underground mining, hydro-compaction and stress given by new constructions, are the principal drivers of the anthropogenic land subsidence [13][14][15][16][17][18][19]. Moreover, the effects of climate change can dramatically increase the subsidence-related problems due to the rising of sea levels: the 2012 Intergovernmental Panel on Climate Change (IPCC, www.ipcc.ch, (accessed on 9 April 2021)) report, in fact, highlight an increasing occurrence of coastal and fluvial flooding,

The Po River Delta
The Po is the largest river in Italy. It opens in the northern sector of the Adriatic Sea with a large delta of about 400 km 2 that extends seawards for about 25 km [28] ( Figure  1). In the delta, the main river (Po di Venezia) is divided into seven branches: Po di Volano, Po di Goro, Po di Gnocca or di Donzella, Po di Tolle, Po di Pila, Po di Maistra and Po di Levante (Figure 1).
The formation of the modern delta is the result of natural processes and human interactions, such as the filling of the wetlands area and engineering endeavors [40]. The eastern part of the delta is mainly characterized by reclaimed territories, as shown in Figure 1. At present day, the reclamation in the PRD area is strongly reduced and part of the territory (about 200 km 2 ) is characterized by a large system of shallow water bodies. Most of the reclaimed territories are actually below the mean sea level and are poorly supplied by sediments because all its river branches have major artificial levees [41]. The complex-wide sandy beach ridges elongating from south to north can be considered as the natural western border of the reclaimed territories ( Figure 1). The compaction of the Plio-Quaternary alluvial deposits in the PRD area is an important driver of the natural ground

The Po River Delta
The Po is the largest river in Italy. It opens in the northern sector of the Adriatic Sea with a large delta of about 400 km 2 that extends seawards for about 25 km [28] (Figure 1). In the delta, the main river (Po di Venezia) is divided into seven branches: Po di Volano, Po di Goro, Po di Gnocca or di Donzella, Po di Tolle, Po di Pila, Po di Maistra and Po di Levante ( Figure 1).
The formation of the modern delta is the result of natural processes and human interactions, such as the filling of the wetlands area and engineering endeavors [40]. The eastern part of the delta is mainly characterized by reclaimed territories, as shown in Figure 1. At present day, the reclamation in the PRD area is strongly reduced and part of the territory (about 200 km 2 ) is characterized by a large system of shallow water bodies. Most of the reclaimed territories are actually below the mean sea level and are poorly supplied by sediments because all its river branches have major artificial levees [41]. The complex-wide sandy beach ridges elongating from south to north can be considered as the natural western border of the reclaimed territories ( Figure 1). The compaction of the Plio-Quaternary alluvial deposits in the PRD area is an important driver of the natural ground settlement in these territories: in particular, [30] suggests that the land subsidence rates in the delta area are significantly correlated with the age of the highly compressible Holocene deposits that compose the shallowest 30-40 m of the sedimentary sequence. Other authors [42,43] obtained similar results by analyzing the correlation between the age of the Holocene sediments and the compaction rate observed in the Southern Louisiana Mississippi delta. Their results show land subsidence rates of about 4-5 mm/year for deposit of ages and thicknesses similar to those present in the PRD.
An additional contribution to the natural land subsidence can be related to tectonic movements: some authors (e.g. [11,44]) suggest a contribution of about 1 mm/year due to a crustal flexing induced by the south-westward subduction of the Adriatic plate under the Apennine belt. This subduction model seems to be not compatible with the evidence provided by seismic surveys, which do not show an evident and well-developed slab [45], and with a maximum depth of more than 90 km of the earthquakes occurring in the Northern Apennines. An alternative model suggests that the present deformation pattern observed in the Apennine belt and the Po Plain is driven by the northward motion of the Adriatic plate, instead of the rollback of the Adriatic subducted margin [46].
The sedimentary basin of the Po Plain is also characterized by a complex system of non-confined, semi-confined, and confined aquifers [47]. The extensive exploitation of these aquifers can be considered as one of the principal causes of the anthropogenic land subsidence observed in both the plain and the delta areas. Because of the important historical and natural patrimony located in the Po Plain and the large concentration of industrial facilities as well as the intensive farming in this area, it is necessary to systematically monitor the occurrence and development of the land subsidence phenomenon.
The first leveling performed in the PRD was at the end of the XIX century (1877-1903) and repeated at the beginning of the 1950s [48]. These surveys were carried out before the Italian economic growth: analyzing these data, the vertical kinematic pattern measured by [49] shows in the eastern sector of the Po Plain values of land subsidence lower than 10 mm/year, while the western part is characterized by even less intense subsidence (up to 2/3 mm/year) and by upward motion (uplift). In detail, the leveling surveys carried out in that period provided average land subsidence rates of about 5-7 mm/year, the highest in the Po Plain; because the leveling measurements were carried out when the contribution of the industrial and agricultural activities of the XX century were negligible, the reported values can be considered representative of the natural contribution to the land subsidence and/or the compaction process due to the land reclamations.
The land subsidence rates increased during the Italian economic growth, after the end of the Second World War: in that period, the multi-aquifers system of the Po Plain was extensively exploited for anthropic uses (industrial and agricultural). The methane-water withdrawals from onshore and offshore reservoirs have also contributed to the increase of the land subsidence values. [28,30,[50][51][52] used precise leveling surveys performed in the 1950-1957 in the PRD and Po Plain to measure land subsidence rates up to 300 mm/year. Subsequently, leveling data obtained from 1962 to 1974, highlighted a progressive reduction of the rates with values of about 30 mm/year, in agreement with the progressive reduction of the anthropogenic withdrawals [50,53]. During the late 1970s and at the beginning of the 1980s, the construction of new public aqueducts exploiting surface waters significantly reduced the aquifer overdraft and yielded a general progressive head recovery with a significant decrease of the land subsidence rate [54]. Recent studies [3,[31][32][33]52,55] indicate a decrease of the land subsidence rates with values less than 10 mm/year, probably as a consequence of the successful policies applied 40 years ago to reduce the anthropogenic deformations in the PRD.

The PODELNET Network
The PODELNET network was developed in 2016 in agreement with the IGMI (Istituto Geografico Militare Italiano) and the Veneto Region (Unità di Progetto per il Sistema Informativo Territoriale e la Cartografia and Unità Organizzativa Genio Civile di Rovigo): it is made of 46 non-permanent sites monitored with the GNSS technique. This low-cost and low-impact network is the first application of the densification at 7 km of the IGM95 network in the Veneto Region and, for this reason, the points are named according to the IGMI procedures.
The PODELNET is used also as a reference frame for the activities carried out by the public authorities in the area, mainly the Veneto Region and the Po Delta reclamation consortium.
The PODELNET sites were chosen to take into account both the existing points belonging to the IGMI leveling network and other points belonging to local Authorities and Institutions: all points are located on stable and consolidated foundations; 16 new sites were also identified to obtain a uniform distribution of the network and avoid gaps. Among the 46 non-permanent GNSS sites of the PODELNET network (Figure 1), 24 were connected with the nearest leveling benchmark, providing the orthometric elevation of the points referred to the last IGMI geometric leveling measurements, made in 2005. The others 22 NPS are the benchmarks of the IGMI leveling network. Each point of the PODELNET is connected with its closest point by using a baseline of between 5 km and 7 km.
The first measurement was performed in June/July 2016, and the second survey was carried out two years later, but in the same months to reduce the potential influence of the seasonal signals in the measured velocities. The amplitudes of these seasonal signals are usually not negligible and can introduce biases in the velocity estimation [56,57]. The two surveys were also designed to reduce the impact of the measurement conditions on the estimation of 46 NPS positions. In particular, the two campaigns were carried out surveying the same baselines, when possible, with the same sampling rates, minimum time stationing, and instruments. The minimum observation time at each site was of 3 h, with a sampling rate of 15 s.

GNSS Data and Analysis
The GNSS observations acquired in 2016 and 2018 were analyzed using the Gamit software (release 10.7) [58]. The Earth Orientation Parameter (EOP) and precise ephemerides provided by the International GNSS Service (IGS) were included in the processing with tight constraints. The IERS/IGS 2003 models were adopted to reconstruct the temporal evolution of diurnal, semidiurnal, and terdiurnal solid Earth tides; the pole tide corrections were applied according to the same IERS standards [58]. The FES2004 model [59] was applied to model the effects of ocean loading. The atmospheric propagation delay was then implemented using the "global mapping function" [60]. We assigned also loose constraints to the coordinates of each station included in the processed network.
The daily loose solutions were translated into a local reference frame through a tridimensional movement, where the three translation parameters were estimated with the Globk-Glorg package software [61], using the coordinates and velocities of three external CGNSS sites located near the PRD area ( Figure 1): CGIA, CODI, and GARI. The daily solutions of each survey were also combined in a unique solution using the Globk-Glorg package software [61] to provide a unique value of position for each PODELNET site.
The data of the two CGNSS sites located inside the PRD (TGPO and PTO1, Figure 1) were also added to the processing of the PODELNET. The velocities of these two sites, obtained using only the observations acquired during the surveys of the network, were compared with those obtained using the entire time-series by means of the IGS sites as reference stations, to evaluate the uncertainties due to the use of a local reference system.
The coordinates and velocities of the five CGNSS sites were measured applying a procedure similar to that described in [21,33,62]. The present kinematic patterns in the Italian peninsula and surrounding regions described in these works were upgraded using the entire time-series available for each site from 1 January 2001 to 31 December 2019, aligned to the ITRF2014 reference frame, and in particular to the ETRF2014 realization, where the horizontal Eurasian plate motion was removed [63]. As suggested by [21,33,62], the time pattern of the north, east and vertical components of daily position can be modeled with the following equation: where j = 1, 2, 3 for the north, east, and vertical components; A j and V j are the intercept and trend/velocity of the best fitting straight-line model, respectively; the D j terms are the N instrumental or seismic discontinuities eventually occurred at the T l epochs; H is the Heaviside step function; B 1jk and B 2jk are the amplitude of the five (M = 5) more significant seasonal terms with period P jkj ; ε j (t) represents the time-dependent noise. The periods of the five more significant seasonal signals were searched in the estimated power spectrum analyzing the residual time-series by the Lomb-Scargle method [64,65]. The residual time-series were obtained removing the estimated linear trend model with only discontinuities (1). These parameters were computed in the first step of the procedure by a weighted least-squares approach (see [21,33,62] for more details).
The spectrum was analyzed finding the period P of the five statistically meaningful signals in the interval between 1 month and half of the observation period. This condition may represent a problem in the procedure because some climatic and geological/geophysical processes could introduce seasonal signals with relatively long periods (greater than half of the observation period) in the GNSS daily time series. Figure 2 shows the power spectrum of the three components related to the five CGNSS sites located in the PRD: it can be noted that some sites are characterized by seasonal signals with a period greater than half of the observation period (about 5-6 years, e.g., GARI). Hence, we have extended the search of the five statistically meaningful signals in the interval between 1 month and the entire observation period. The velocities obtained using equation 1 and associated uncertainties using the method suggested by [66] are reported in Table 1, and the values of the surrounding sites (CGIA, CODI, and GARI) are used to translate in the same reference frame the measured PODELNET positions.
The values were obtained under the hypothesis that the velocity is constant in the entire period of analysis; these periods are greater than the one between the two PODELNET surveys. With this hypothesis, the velocities obtained using the entire period are equal to the rate between the two PODELNET campaigns. This assumption cannot always be satisfied when the GNSS observations are used to monitor land subsidence because the anthropogenic contributions can cause quick changes in the lowering of the ground surface and, therefore, in the land subsidence velocity. These changes can have a limited temporal extension, for example, less than the recommended 2.5 years, which is the minimum period suggested to obtain reliable GNSS velocities [67]: for these reasons, the values shown in Table 1 could not represent the real velocities occurred between the two PODELNET surveys; moreover, the use of these data in the referencing procedure could introduce not negligible biases in the estimated land subsidence velocities for the PODELNET points. We attempted to overcome this problem by using the Velocity Moving Window Approach (VMWA), as suggested by other authors [24,68].
The daily time-series were analyzed with a moving window of 755 days (the time span between the two PODELNET campaigns) to obtain the "instantaneous" velocity (IV): however, the seasonal signals shown in Figure 2 can introduce noise and biases in the IVs, because they can represent local phenomena not representative of a common signal in the PRD area; for this reason, we analyzed the daily time-series where the meaningful significant (M) seasonal signals obtained with equation 1 were removed. Figure 3 shows the IV time series obtained by a window of 755 days, rejecting the velocities obtained with observations shorter than 300 days. The obtained IV corresponding to the PODELNET surveys is also shown in Table 2. The daily time-series were analyzed with a moving window of 755 days (the time span between the two PODELNET campaigns) to obtain the "instantaneous" velocity (IV): however, the seasonal signals shown in Figure 2 can introduce noise and biases in the IVs, because they can represent local phenomena not representative of a common signal in the PRD area; for this reason, we analyzed the daily time-series where the meaningful significant (M) seasonal signals obtained with equation 1 were removed. Figure 3 shows the IV time series obtained by a window of 755 days, rejecting the In Table 2 are shown (Campaigns columns) the velocities obtained with only the CGNSS observations acquired during the PODELNET surveys. The differences between the values obtained with VMWA and the ones obtained using the observations acquired during the PODELNET measurements are less than 2 mm/year. Table 1. Velocities in mm/year of CGNSS stations located in the PRD and surrounding areas, adopted in the processing of the network as references (CGIA, CODI, GARI) and benchmarks (PTO1, TGPO). In the first and second columns are reported the CGNSS station codes and the date (T), in decimal year, of the first observation used in the processing, respectively. V are the velocities in the ETRF2014 reference frame [63] of the North (V N ), East (V E ), and Vertical (V V ) components. The uncertainties associated with the velocities were estimated using the method suggested by [66]. The Weighted Root Means Square values (σ) in mm/year of the analyzed time-series are shown in the last three columns.

Code
Start Remote Sens. 2021, 13, x FOR PEER REVIEW 8 of 21 velocities obtained with observations shorter than 300 days. The obtained IV corresponding to the PODELNET surveys is also shown in Table 2. Table 1. Velocities in mm/year of CGNSS stations located in the PRD and surrounding areas, adopted in the processing of the network as references (CGIA, CODI, GARI) and benchmarks (PTO1, TGPO). In the first and second columns are reported the CGNSS station codes and the date (T), in decimal year, of the first observation used in the processing, respectively. V are the velocities in the ETRF2014 reference frame [63] of the North (VN), East (VE), and Vertical (VV) components. The uncertainties associated with the velocities were estimated using the method suggested by [66]. The Weighted Root Means Square values (σ) in mm/year of the analyzed time-series are shown in the last three columns.  Table 1. The purple vertical line represents the data related to the last day of the second PODELNET survey (18 July 2018). The IV associated with this data are the values corresponding to the time span between the two measurements.
In Table 2 are shown (Campaigns columns) the velocities obtained with only the CGNSS observations acquired during the PODELNET surveys. The differences between the values obtained with VMWA and the ones obtained using the observations acquired during the PODELNET measurements are less than 2 mm/year.  Table 1. The purple vertical line represents the data related to the last day of the second PODELNET survey (18 July 2018). The IV associated with this data are the values corresponding to the time span between the two measurements. Table 3 and Figure 4 shows the velocities of the PODELNET NPS sites obtained by means of the differences between the 2016 and 2018 coordinates: the uncertainties associated with these values cannot be estimated using the methods proposed in the literature (e.g., [66,69]), where the time-correlated noise in the GNSS time-series is taken into account. We adopted a conservative uncertainty of about 4 mm/year, twice the maximum difference reported in Table 2. The velocities of the CGNSS sites located in the study area shown a kinematic pattern characterized by velocities of about 1-2 mm/year in the north direction (Table 1), as also found in other studies (e.g., [21,33,46]). The uncertainties Remote Sens. 2021, 13, 1488 9 of 21 associated with the PODELNET velocities are greater than the horizontal kinematic values, which suggests that these rates cannot be used for a detailed study of the horizontal velocity field. Table 2. Velocities of the three components (North V N , East V E, and Vertical V V ), in mm/year, of the continuous GNSS stations estimated by means of the Velocity Moving Window Approach (VMWA columns) using a window of 755 days (from 21 June 2016 to 18 July 2018). The values of the CGNSS sites estimated using only the observations acquired during the two campaigns are shown in the Campaigns columns, in mm/year; in the last three columns (Differences) are reported the differences between the corresponding values. The uncertainties associated with the VMWA velocities were obtained by a weighted least square approach adopted in the estimation of the velocities.  Np is the number of InSAR points used to obtain the V S rate, and R is the minimum radius, in meters, adopted to find the number of points. The values of PTO1 and TGPO permanent sites were obtained using only the observations acquired during the two PODELNET surveys and processed with the data of these campaigns.

InSAR Data and Analysis
The data and results used in this study are the ones already published in [3]. In that work, the authors analyzed three different InSAR data sets with the Small Baseline Subset (SBAS) [70] approach over the entire PRD area. In particular, they obtained mean velocities and displacement time-series from the C-band ERS-1/2 (1992-2000), ENVISAT (2004-2010), and Sentinel-1A/B (2014-2017) satellites. The SBAS approach takes into account the socalled distributed scatterers, which are objects on the ground with similar backscattering properties within a SAR image resolution cell (ground pixel). For each of the detected ground pixels, the time-series analysis is performed on the stack of available images by following five main steps: (a) each image is connected multiple times (high redundancy) to generate a well-connected interferogram network; (b) the interferograms are generated and then co-registered to a selected reference image; (c) the phase of the pixels with coherence higher than a selected threshold is unwrapped using a Delaunay Minimum Cost Flow (MCF) method [71]; (d) the atmospheric phase components (Atmospheric Phase Screen) are estimated and removed by using low-pass and high-pass filters; (e) the velocities and displacements time-series are finally calculated and classified in an ArcGIS environment. For further details on the adopted processing workflow and processing parameters please refer to [3].
Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 21 calculated and classified in an ArcGIS environment. For further details on the adopted processing workflow and processing parameters please refer to [3].

Results of the GNSS-InSAR Integrated Monitoring
Among the available InSAR datasets, only the Sentinel-1A/B entirely overlaps with the GNSS observations obtained from the CGNSS stations and the PODELNET network. To compare the InSAR-and GNSS-derived velocities, it was necessary to align the two datasets into the same reference frame. Since the InSAR images do not entirely cover the study area and the CGIA and GARI sites are not covered with InSAR data, we have chosen the GNSS velocities of CODI to translate the InSAR velocities from LOS to the GNSS reference system. The velocities of PTO1 and TGPO stations were used to evaluate the accuracy of the translation process.
The first step of the GNSS-InSAR integrated monitoring procedure is to estimate the GNSS values taking into account only the period overlapping with the InSAR observations: we analyzed the daily time-series without discontinuities and principal seasonal signals with the MVWA method, using a window of 1128 days, which is the time span between the first and the last InSAR observations (17 November 2014-17 May 2017). In Figure 5 are reported the obtained IV time-series.
A limitation of the procedure to align the InSAR velocity field with a local GNSS reference frame is related to the estimation of the InSAR rate at the CGNSS site: a simple method is to find the InSAR point closest to the CODI site and compare the two velocities to measure the translation value. However, these points could not be exactly located at the CODI station and therefore the local movements of the InSAR points are potentially not representative of the CGNSS site, or noisy InSAR data could introduce biases in the alignment procedure. To avoid these potential sources of error, we adopted an averaging procedure, where the velocities of the InSAR points located at a distance less than a defined search radius are used to obtain the InSAR velocity. A critical point is to find the correct size of the search radius: values obtained with a small radius can be strongly influenced by outliers, while a large radius can include other processes not belonging to the CGNSS site.

Results of the GNSS-InSAR Integrated Monitoring
Among the available InSAR datasets, only the Sentinel-1A/B entirely overlaps with the GNSS observations obtained from the CGNSS stations and the PODELNET network. To compare the InSAR-and GNSS-derived velocities, it was necessary to align the two datasets into the same reference frame. Since the InSAR images do not entirely cover the study area and the CGIA and GARI sites are not covered with InSAR data, we have chosen the GNSS velocities of CODI to translate the InSAR velocities from LOS to the GNSS reference system. The velocities of PTO1 and TGPO stations were used to evaluate the accuracy of the translation process.
The first step of the GNSS-InSAR integrated monitoring procedure is to estimate the GNSS values taking into account only the period overlapping with the InSAR observations: we analyzed the daily time-series without discontinuities and principal seasonal signals with the MVWA method, using a window of 1128 days, which is the time span between the first and the last InSAR observations (17 November 2014-17 May 2017). In Figure 5 are reported the obtained IV time-series.
A limitation of the procedure to align the InSAR velocity field with a local GNSS reference frame is related to the estimation of the InSAR rate at the CGNSS site: a simple method is to find the InSAR point closest to the CODI site and compare the two velocities to measure the translation value. However, these points could not be exactly located at the CODI station and therefore the local movements of the InSAR points are potentially not representative of the CGNSS site, or noisy InSAR data could introduce biases in the alignment procedure. To avoid these potential sources of error, we adopted an averaging procedure, where the velocities of the InSAR points located at a distance less than a defined search radius are used to obtain the InSAR velocity. A critical point is to find the correct size of the search radius: values obtained with a small radius can be strongly influenced by outliers, while a large radius can include other processes not belonging to the CGNSS site. by a native ground resolution of about 20 m × 5 m while the SBAS technique used to process the data takes into account the averaged contribution of many independent scatterers within a larger ground resolution cell of about 20 m × 20 m [3]. Thus, we tested the stability of the averaged values at different radiuses, starting from 20 m (Figure 6a): it can be noted that increasing the radius size, the InSAR velocity increases, even if the differences are less than the weighted root mean square values adopted as uncertainties associated with the average values ( Figure 6a). As shown in Figure 6b, the number of InSAR points in the search domain increases significantly with the size of the radius, but the average distance between the points and the CGNSS site is less than the search radius (Figure 6c). The relatively high number of InSAR points located in the surroundings of the CGNSS site (Figure 6b) can be explained by the presence of the urban area where the CGNSS station is located. Figure 5. IV time-series of the CGNSS sites located in the PRD area obtained using the VMWA with a window of 1128 days (the time span between the first and last InSAR images processed). The windows containing less than 500 days were rejected. We have analyzed the time series using the model described in equation 1 without discontinuities and significant seasonal signals and estimated by a weighted least square method. The black bars represent the uncertainties associated with the IV obtained with the weighted least square method used to compute the rates. The horizontal red lines represent the velocities obtained with the entire time series and are reported in Table 1. The purple vertical line represents the data related to the last day of the InSAR observations (17 May 2017). The IV associated with this data are the values corresponding to the time span between the InSAR images.
We found that the 30 m radius was the one that better represents the CODI site with the InSAR velocities considering that the differences between the different radiuses are less than the weighted root mean square values assumed as uncertainties associated with the average values (of around 0.6 mm/year) (Figure 6a). Figure 5. IV time-series of the CGNSS sites located in the PRD area obtained using the VMWA with a window of 1128 days (the time span between the first and last InSAR images processed). The windows containing less than 500 days were rejected. We have analyzed the time series using the model described in equation 1 without discontinuities and significant seasonal signals and estimated by a weighted least square method. The black bars represent the uncertainties associated with the IV obtained with the weighted least square method used to compute the rates. The horizontal red lines represent the velocities obtained with the entire time series and are reported in Table 1. The purple vertical line represents the data related to the last day of the InSAR observations (17 May 2017). The IV associated with this data are the values corresponding to the time span between the InSAR images.
The InSAR observations acquired by the Sentinel-1A/B satellites are characterized by a native ground resolution of about 20 m × 5 m while the SBAS technique used to process the data takes into account the averaged contribution of many independent scatterers within a larger ground resolution cell of about 20 m × 20 m [3]. Thus, we tested the stability of the averaged values at different radiuses, starting from 20 m (Figure 6a): it can be noted that increasing the radius size, the InSAR velocity increases, even if the differences are less than the weighted root mean square values adopted as uncertainties associated with the average values ( Figure 6a). As shown in Figure 6b, the number of InSAR points in the search domain increases significantly with the size of the radius, but the average distance between the points and the CGNSS site is less than the search radius (Figure 6c). The relatively high number of InSAR points located in the surroundings of the CGNSS site (Figure 6b After the correction of the InSAR kinematic pattern in the local GNSS reference system (at CODI site), we have compared the InSAR and the GNSS velocities of the two benchmark sites (PTO1 and TGPO) using the selected 30 m radius, which allowed us to find at least 6 InSAR points. The obtained InSAR velocities in correspondence of the benchmark sites are shown in Table 4 together with the GNSS velocities along the LOS direction obtained using the formulas suggested by [73]: it can be noted that even if the differences with the GNSS values are less than 1 mm/year, we have adopted a conservative threshold of 2 mm/year as the uncertainty associated with the InSAR velocities translated in the GNSS reference frame.
We obtained the average InSAR velocities at each GNSS PODELNET station by selecting at least 6 InSAR points enclosed in a search circle centered on the NPS sites with a radius from 20 m to 300 m and with a step of 10 m. The NPS velocities along the LOS and the differences between these rates and the InSAR values are shown in Table 3: it can be noted that for 7 sites the InSAR velocities are not reported because the number of InSAR points enclosed in the largest search circle is less than 6, probably as a consequence of the high vegetation coverage of these areas. We found that the 30 m radius was the one that better represents the CODI site with the InSAR velocities considering that the differences between the different radiuses are less than the weighted root mean square values assumed as uncertainties associated with the average values (of around 0.6 mm/year) (Figure 6a).
After the correction of the InSAR kinematic pattern in the local GNSS reference system (at CODI site), we have compared the InSAR and the GNSS velocities of the two benchmark sites (PTO1 and TGPO) using the selected 30 m radius, which allowed us to find at least 6 InSAR points. The obtained InSAR velocities in correspondence of the benchmark sites are shown in Table 4 together with the GNSS velocities along the LOS direction obtained using the formulas suggested by [73]: it can be noted that even if the differences with the GNSS values are less than 1 mm/year, we have adopted a conservative threshold of Table 4. Velocities of the three components (North V N , East V E , and Vertical V V ), in mm/year, of continuous GNSS stations obtained using the Velocity Moving Window Approach (VMWA columns) with a window of 1128 days. The velocities of the site CODI (Figure 1) were used in the reference translation procedure to align the InSAR velocities in the GNSS reference frame. V L , in mm/year, is the GNSS velocity along the LOS direction. V S , in mm/year, is the InSAR rate obtained averaging all the points in a 30 m radius. The differences between the InSAR rate and the LOS GNSS velocities of the benchmarks sites are also reported, in mm/year, in the V D column.

GNSS on LOS InSAR LOS Differences
We obtained the average InSAR velocities at each GNSS PODELNET station by selecting at least 6 InSAR points enclosed in a search circle centered on the NPS sites with a radius from 20 m to 300 m and with a step of 10 m. The NPS velocities along the LOS and the differences between these rates and the InSAR values are shown in Table 3: it can be noted that for 7 sites the InSAR velocities are not reported because the number of InSAR points enclosed in the largest search circle is less than 6, probably as a consequence of the high vegetation coverage of these areas.
The minimum number of 6 InSAR points included in the domain around the NPS GNSS sites is not always reached with small circles. For some of the sites, in particular for the ones located in the poorly urbanized territories of the PRD, the minimum search radius necessary to obtain an adequate number of points is greater than 150 m (Figure 7b).   Figure 1) were used in the reference translation procedure to align the InSAR velocities in the GNSS reference frame. VL, in mm/year, is the GNSS velocity along the LOS direction. VS, in mm/year, is the InSAR rate obtained averaging all the points in a 30 m radius. The differences between the InSAR rate and the LOS GNSS velocities of the benchmarks sites are also reported, in mm/year, in the VD column.

GNSS on LOS InSAR LOS Differences
The minimum number of 6 InSAR points included in the domain around the NPS GNSS sites is not always reached with small circles. For some of the sites, in particular for the ones located in the poorly urbanized territories of the PRD, the minimum search radius necessary to obtain an adequate number of points is greater than 150 m (Figure 7b). Furthermore, we have investigated if the differences in velocities found in respect to the NPS points are linked to the selected search radius or to the number of InSAR points used to calculate the InSAR velocities: Figure 7c,d shows that the differences are not related to the radius or the number of points included in the domain. Figure 8a shows the InSAR kinematic pattern in the GNSS local reference frame. Furthermore, we have investigated if the differences in velocities found in respect to the NPS points are linked to the selected search radius or to the number of InSAR points used to calculate the InSAR velocities: Figure 7c,d shows that the differences are not related to the radius or the number of points included in the domain. Figure 8a shows the InSAR kinematic pattern in the GNSS local reference frame.
The differences between the InSAR and the NPS GNSS velocities along the LOS are shown in Figure 8b. The sites characterized by high positive differences (VInSAR -VGNSS), correspond to the points where the PODELNET data observed high rates of land subsidence. These sites shown in Figure 4a with yellow and green points are probably affected by noise or bias due to some problems that occurred during the surveys. These relatively high velocities can also be related to very localized land subsidence not captured by the InSAR technique. Only future GNSS measurements can provide the information necessary to understand which of the previous two hypotheses are correct.
In the Discussion section we take into consideration only the PODELNET sites where the differences between the InSAR and the GNSS values are less than 10 mm/year (Figure 8b).  (Table 3) overlapped to the InSAR kinematic pattern using the same scale; (b) Differences between InSAR and GNSS velocities (Table 3) at the PODELNET sites (VInSAR − VGNSS); are also shown the positions of the seven Po River terminal branches.

Discussion
The vertical kinematic pattern shown in Figure 4 has highlighted several aspects of the spatial and temporal distribution of the land subsidence in the PRD area. The velocities of the analyzed CGNSS sites (Table 1), obtained using the observations acquired in the last decade, are significantly less than -10 mm/year: the values related to the vertical component of PTO1 and TGPO CGNSS stations, are in agreement with those obtained by [55]. The comparison between these values and the ones reported by other authors (e.g., [21,30,31,49,74]) using observations acquired between the past centuries and the first years of the 2000s, clearly indicates a decrease in the land subsidence rate in the PRD area.
The land subsidence rates of about 6 mm/year measured in the last decade (Table 1) agree with the values linked with the consolidation of the late Holocene sediments [30,41,42,52]. However, anthropogenic activities can also contribute to these rates providing variations with periods shorter than the ones of the available CGNSS observations. For this reason, we have investigated the possibility of having short changes of the CGNSS velocities during the entire observation period using the VMWA method.  (Table 3) overlapped to the InSAR kinematic pattern using the same scale; (b) Differences between InSAR and GNSS velocities (Table 3)  The differences between the InSAR and the NPS GNSS velocities along the LOS are shown in Figure 8b. The sites characterized by high positive differences (V InSAR -V GNSS ), correspond to the points where the PODELNET data observed high rates of land subsidence. These sites shown in Figure 4a with yellow and green points are probably affected by noise or bias due to some problems that occurred during the surveys. These relatively high velocities can also be related to very localized land subsidence not captured by the InSAR technique. Only future GNSS measurements can provide the information necessary to understand which of the previous two hypotheses are correct.
In the Discussion section we take into consideration only the PODELNET sites where the differences between the InSAR and the GNSS values are less than 10 mm/year (Figure 8b).

Discussion
The vertical kinematic pattern shown in Figure 4 has highlighted several aspects of the spatial and temporal distribution of the land subsidence in the PRD area. The velocities of the analyzed CGNSS sites (Table 1), obtained using the observations acquired in the last decade, are significantly less than -10 mm/year: the values related to the vertical component of PTO1 and TGPO CGNSS stations, are in agreement with those obtained by [55]. The comparison between these values and the ones reported by other authors (e.g., [21,30,31,49,74]) using observations acquired between the past centuries and the first years of the 2000s, clearly indicates a decrease in the land subsidence rate in the PRD area.
The land subsidence rates of about 6 mm/year measured in the last decade (Table 1) agree with the values linked with the consolidation of the late Holocene sediments [30,41,42,52]. However, anthropogenic activities can also contribute to these rates providing variations with periods shorter than the ones of the available CGNSS observations. For this reason, we have investigated the possibility of having short changes of the CGNSS velocities during the entire observation period using the VMWA method. This approach provided also the velocities necessary to align the PODELNET results and the InSAR LOS rates in the same local reference frame, using only the daily positions of the CGNSS stations acquired in the overlapping period between the NPS GNSS surveys (Table 2) and the InSAR observations (Table 4) The vertical components of the VMWA time series of the CGNSS sites show some significant variations (Figures 3 and 5), especially for the stations located outside the PRD area (CGIA, CODI, and GARI). The IV time series obtained with the shorter window ( Figure 3) show more significant changes compared to the ones obtained with the 1128 days window ( Figure 5). Climatic and/or meteorological processes and human activities are usually characterized by seasonal and/or non-periodic signals with variable amplitudes and relatively short periods (e.g., few months and/or years): these characteristics can be detected by the MVWA method providing variations in the IV time-series obtained with short windows; the analysis performed using large windows can attenuate or remove the variations related to these phenomena in the 'instantaneous' time-series. The variations observed in the PTO1 and TGPO sites seem to present these characteristics: the time-series obtained with short windows (755 days, Figure 3) show some variations in the vertical velocities, while the series obtained with large windows (1128 days, Figure 5) do not show significant changes. Similar results were obtained analyzing the time series of the CGIA site, whereas the other two stations, CODI and GARI, are characterized by significant variations also in the IV time series obtained with a large window. The signals in the GARI and CODI time-series can be linked to the locations of the sites, which can highlight climatic processes or local human activities: for example, the GARI CGNSS site is located on the quay of the Porto Garibaldi harbor and, for this reason, the changes in the vertical velocities could be connected to sea-level changes or local movements of the site related to the tide levels. The analysis of these particular displacements is one of the goals of our future researches: however, the difficulty to obtain adequate information in these sites is one of the major challenges to better understand these temporal variations.
Results of Table 1 (CGNSS sites) indicate also a small but dominant N-S movement, not detectable with the InSAR data. The preliminary horizontal PODELNET velocities reported in Table 3 could indicate that the N-S component is not always dominant. The E-W component in some sites can be possibly related to technical problems that occurred during the measurements or to local effects; therefore, in this case, the comparison with the InSAR velocities can highlight the sites affected by such problems. On the other hand, land subsidence processes can introduce not negligible local E-W movements that can provide a significant contribution to the measured InSAR LOS velocity: in this case, the information inferred from the observations acquired by a ground-based GNSS network can help to detect the areas not covered by InSAR, and/or to remove possible biases in the InSAR vertical velocity field. Figure 4a shows the preliminary vertical kinematic pattern obtained analyzing the two PODELNET surveys: the green and yellow sites, characterized by velocities greater than 15/20 mm/year, have not been considered because they are most likely outliers: particular attention must be taken for the future measurements to verify if these high velocities are linked to a particular land subsidence process or just to a technical problem occurred during the surveys. The comparison between the InSAR and GNSS data along the LOS direction shows differences lower than 10 mm/year (Figure 8b): it can be noted that the NPS sites located on the complex sandy beach ridge of the PRD area (Figures 1 and 4a) are characterized by positive velocities or low/moderate negative rates (up to −2 mm/year), whereas the ones on the land reclamation and sedimentary areas are characterized by more heterogeneous velocity patterns, with land subsidence rates lower than 5 mm/year (Figure 4a). However, the obtained kinematic pattern is in agreement with both the results provided by [55] using COSMO-SkyMed and ALOS-PALSAR InSAR data and those of [30] derived from ERS-1/2 InSAR data. Figure 4b shows the comparison between the vertical kinematic pattern and the highstand PRD lobe map: here there is a good correlation between the PODELNET land subsidence rates and the "younger" territories (0-200 B.P. and 200-350 B.P.) (Figure 4b, dark yellow and pink areas). Land reclamation activities in these areas were carried out by the Republic of Venice starting from the XVII century, with the main intervention that forced the Po river to flow southward through an artificial delta mouth to prevent the sedimentary infilling of the Venetian Lagoon [72]. These results indicate that the land subsidence is mainly due to the compaction of the depositional elements, organic-rich clay, and fresh-water peat located in the shallowest layers of the stratigraphic succession [72]. A similar relationship between sediment age and the land subsidence velocity was found by [30]: they provided land subsidence rates of about 2 mm/year on the sandy beach ridges, and higher velocities (10-15 mm/year) in the eastern deltaic area, in agreement with the results obtained in this work. These authors assumed that only vertical movements contributed to the LOS displacements measured with InSAR: however, the preliminary horizontal velocity field reported in Table 3 indicates that in some areas this hypothesis may not be verified.
The area of the Po di Goro mouth, located in the southern sector of the PRD on 'younger' territories, seems to be characterized by low values of land subsidence, as observed in three different PODELNET sites (Figure 4a,b). A similar pattern can be observed in the Po di Tolle mouth area, where two PODELNET sites present low values of the land subsidence rates (Figure 4a,b). Similar rates are also visible in other two PODELNET sites, located in the area of Albarella Island (Figure 4a,b) and in agreement with the InSAR data: these NPS are located on a sandy beach ridge (Figures 1 and 4a) that correspond to the territories reclaimed during the reorganization of the drainage system occurred in the Medieval age (Figure 4b): this result confirms the previous conclusion about the NPS located on the sandy beach ridges in the western portion of the study area, i.e., the territories reclaimed before the XV-XVI century. Figure 8a shows the PODELNET and LOS InSAR kinematic patterns: it can be noted that the points located on the sandy beach ridges (Figure 1) and in the western PRD area are characterized by LOS rates lower than the ones located on the easternmost "younger" [72] reclaimed territories, in agreement with the results provided by [30,55]; this confirms that the main contribution to the land subsidence is the compaction of the superficial layers caused by the recent reclamation processes.

Conclusions
The present land subsidence rates in the PRD area were monitored by integration of CGNSS, NPS GNSS, and InSAR data. The results obtained analyzing 5 CGNSS stations show vertical land motions of about −6 mm/year in the deltaic area and about −3 mm/year in the surrounding territories, without significant variations in the last decade. The preliminary PODELNET NPS GNSS vertical kinematic pattern, obtained analyzing the data acquired through a network of 46 sites in the 2016 and 2018 surveys, shows a heterogeneous land subsidence pattern in the PRD area: the lower vertical land motion rates (values greater than −5 mm/year) are located on the sandy beach ridge that corresponds to the Bronze/Iron Age stranded beach ridges. On the contrary, the higher land subsidence velocities are measured in the easternmost PRD territories formed during the last three centuries. Some areas located in the south-eastern sector of the PRD are characterized by low subsidence velocities (Po di Goro and Po di Tolle mouths): this could indicate more efficient management of the land subsidence processes, but further studies are needed to confirm such hypothesis. Considering the fact that in general the InSAR technology does not perform well in high vegetated areas or areas with high temporal variability, the 46 PODELNET sites represent an important improvement for the monitoring of the land subsidence in the PRD and can represent a strategic monitoring tool for the management of the PRD area in the next future. A scheduled continuous campaign and/or the transformation of some sites in semi-permanent continuous GNSS sites can provide further useful information to better understand the deformations in these territories. An integrated monitoring system, combining GNSS and InSAR data, allows to overcome the limits of the two techniques: the processed InSAR images, validated using the available CGNSS stations, can be integrated with the points of the PODELNET network, increasing the spatial and the temporal resolution together with the cost-effectiveness of both approaches. This integrated monitoring system find considerable applications in coastal management, especially for the safeguarding of natural ecosystems and human activities; the monitoring of the defense systems against flooding is crucial in areas largely below the mean sea level. Further developments of this study will be linked to the monitoring of the 350 km-long primary levees network of the eastern PRD using new high-resolution InSAR data and GNSS observations. Finally, the method presented in this work, which includes the integration between different types of data for the evaluation and monitoring of land subsidence, can be applied worldwide in other coastal areas where significant vertical land motion is expected as a consequence of human activities and increasing rates of sea-level rise due to climate changes.