Multi-Component and Multi-Source Approach for Studying Land Subsidence in Deltas

The coupled effects of climate change and land sinking make deltas and coastal areas prone to inundation and flooding, meaning that reliable estimation of land subsidence is becoming crucial. Commonly, land subsidence is monitored by accurate continuous and discrete measurements collected by terrestrial and space geodetic techniques, such as Global Navigation Satellite System (GNSS), Interferometry Synthetic Aperture Radar (InSAR), and high precision leveling. In particular, GNSS, which includes the Global Positioning System (GPS), provides geospatial positioning with global coverage, then used for deriving local displacements through time. These site-positioning time series usually exhibit a linear trend plus seasonal oscillations of annual and semi-annual periods. Although the periodic components observed in the geodetic signal affect the velocity estimate, studies dealing with the prediction and prevention of risks associated with subsidence focus mainly on the permanent component. Periodic components are simply removed from the original dataset by statistical analyses not based on the underlying physical mechanisms. Here, we propose a systematic approach for detecting the physical mechanisms that better explain the permanent and periodic components of subsidence observed in the geodetic time series. It consists of three steps involving a component recognition phase, based on statistical and spectral analyses of geodetic time series, a source selection phase, based on their comparison with data of different nature (e.g., geological, hydro-meteorological, hydrogeological records), and a source validation step, where the selected sources are validated through physically-based models. The application of the proposed procedure to the Codigoro area (Po River Delta, Northern Italy), historically affected by land subsidence, allowed for an accurate estimation of the subsidence rate over the period 2009–2017. Significant differences turn out in the retrieved subsidence velocities by using or not periodic trends obtained by physically based models.


Introduction
Deltas are relevant for the life and economic development of a country, harbouring small towns and megacities, feeding agricultural activities and fishing, and hosting important ecosystems. These complex natural systems, although modified over the centuries by human interventions, are founded on a delicate water-land balance that strongly depends on climatic variability effects (e.g., frequent and intense rainstorms, peaks in river discharge, and sea level fluctuations). Such effects, coupled with geodetic data by statistical models, but rarely apply the peering approach [21]. Such methodology allows the quantification of the influence, distribution and magnitude of individual processes and the understanding of the main mechanisms affecting the geodetic time series.
In contrast to the above-mentioned approaches, a multi-component and multi-source procedure is proposed to individuate the physical mechanisms that might better explain both the permanent and the periodic subsidence components, clearly exhibited by the geodetic time series. The approach is multi-component because it selectively investigates the linear component separately from the seasonal ones, and is multi-source as it allows us to define the key processes (sources) causing subsidence. The procedure is structured in steps based on multi-methodological and multi-disciplinary analyses, and on the modeling of multiple processes. It breaks through the one-component investigation and improves the estimate of the geodetic velocity, fixes the "analyses to modeling" procedure, thus enhancing the currently used qualitative or semi-quantitative procedures, and attributes the source relevance before the modeling phase, thus speeding up the source validation step.

Materials and Methods
The new approach for studying land subsidence in deltas is described in the following section and then applied to the southern sector of the Po River Delta (Northern Italy).

Multi-Component and Multi-Source Approach
The proposed multi-component and multi-source approach is based on comparative analyses between geodetic time series and datasets of different nature (e.g., geological, hydro-meteorological, hydrogeological, and mareographic data) to find correlations that can guide the selection of the possible multiple processes responsible for the observed ground subsidence, which are then validated through physically-based models.
Two general aspects should be introduced. Firstly, most of the studies described above associate the linear trend (permanent component) of the considered time series to the measured subsidence (e.g., [11][12][13]). This assumption fails when, in addition to inelastic deformations responsible for constant rates, elastic strains also occur, such as in the case of water withdrawal (e.g., [14]). To account for this occurrence, our approach investigates not only the permanent component of the ground vertical movement, but also the periodic ones through statistical and spectral analyses.
The second key aspect of the procedure concerns the choice of the natural physical processes responsible for the observed ground vertical movements. Given the strong dependence of such movements on the geological setting (e.g., [22]) and the climate effects (e.g., [23,24]), the choice of relevant sources is made according to the geology and geography of the study area. In the specific case of deltas, weather-dependent processes are the main sources to be considered. In our approach, such sources are grouped in water mass-and pressure-dependent processes. The former group is associated with the elastic response of the Earth's surface to the water mass change (rain, river, and sea water), while the latter is related to the elastic compaction and expansion of the fine-grained thin layers within the aquifer, caused by pore-pressure changes (groundwater flow).
The proposed methodological path develops through three main steps illustrated in Figure 1. Standard statistical methods and wavelet analysis are first applied to the geodetic time series for extracting the individual components (Step 1 or component recognition). Then, standard linear and non-linear techniques (e.g., cross wavelet transform and wavelet transform coherence analyses) are used for comparing geodetic data with datasets of different nature to find correlations between land and surface/subsurface systems and to identify possible sources (Step 2 or source selection). Finally, the latter are validated through physically based models (Step 3 or source validation). The choice of the most suitable analytical technique to perform Steps 1 and 2 ( Figure 1) depends mainly on the sampling rate of the available data, but also on the scale of the process to investigate.
For instance, wavelet analyses are suitable for signals with constant sampling rates and allow obtaining multi-scale data details by comparing simultaneously local amplitudes, frequencies, and phases of the analyzed signal patterns (e.g., [25]). Other linear statistical techniques, such as the centered moving average (CMA), can be applied to data with inhomogeneous sampling rates for investigating multi-scale dependences (e.g., [26]). The main drawback of the CMA is that the obtained trends extend over time lengths shorter than the original one. The outcomes can then be extended to the original time span using function estimates. An example of such a kind of function will be presented in Section 3, wherein a natural smoothing spline will be used for extending the CMA.
Furthermore, since non-stationary processes are supposed to occur at least on a local scale, the wavelet transforms are herein preferred with respect to other mathematical techniques that analyze the periodicities in frequency domain (e.g., Fourier analysis) by assuming that the underlying processes are stationary in time. Conversely, thanks to the expansion of the time signals into time-frequency space, the wavelet transform allows the identification of both stationary and local intermittent periodicities [25].
To estimate the periodic component of land subsidence, it is worth considering the correlation/ anti-correlation expected from comparing geodetic time series and hydro-meteorological patterns. In principle, in the context of delta areas, it is usually assumed that the occurrence of the water pressureand mass-dependent processes is mainly driven by the lithological composition of soil (topsoil and subsoil) and shallow layers of sedimentary coverage (bedrock). Figure 2 shows a simplified lithological sketch of subsurface layers (lower plot) with the associated physical processes (middle plot) and the ideal ground level rebound in response to the increase and decrease of water supply (upper plot).
In the presence of permeable soil and bedrock with interbedded fine-grained sediments (left side of Figure 2), the infiltration of rain/river water might become dominant in controlling the groundwater table variations. Within the shallow aquifer, the pressure-dependent mechanisms may be generated and the groundwater flow may be accompanied by the expansion and contraction of fine sediments. A peak of water supply might be associated with an uplift of ground level. Conversely, in the case of impermeable soil (i.e., clay-rich soil) (right side of Figure 2), water infiltration may be almost inhibited, while water runoff and loading on the Earth's surface dominate.
Furthermore, in the case of permeable topsoil, water infiltration may occur only weakly, changing the moisture weight of the topsoil and favoring the loading process above the subsoil. Both these cases may be ascribable to water mass-dependent processes, and the peak water supply may correspond to a depression in the ground level. Moreover, in the case of water supply in partially or fully impermeable soil extended homogeneously from floodplains to sea and cut by a great river, the loading effect may be due to river and sea water mass-dependent processes, where riverbanks might be influenced by river water weight oscillations, as well as littoral zones by sea level fluctuations. In these cases, a peak in the river or sea level also corresponds to a ground down lift. It is worth noting that the proposed methodology on the site-positioning time series allows us to investigates the physical sources of the ground motion, and it is quite different from the standard method based on the simultaneous estimation of linear trend, annual, and semi-annual sinusoids with constant phase and amplitude [15,27]. According to this procedure, indeed, the nature of the seasonal oscillations is not investigated, because the periodic terms (e.g., atmospheric and hydrological loading or human exploitation) are modeled through statistical approaches and filtering techniques.

Study Area
The Po River Delta ( Figure 3) is an area historically affected by subsidence due to natural and anthropogenic causes (e.g., [28][29][30]) and also influenced by global climatic changes [31]. The area is characterized by lower delta plain, with ancient coast, delta fronts and lagoons constituted by inter-distributary bay deposits, marsh deposits, beach-ridge, and aeolian dune sediments [32]. The modern Po Delta was generated 400 years ago by an artificial fluvial-mouth cut done in Porto Viro [33], and it is evolving today towards a wave-dominated and cuspate shape [34]. The delta plain overlays 10-30 meters of deltaic, littoral and marine sediments, deposed during the last Holocene prograding phase [35]. The different deltaic environments influence also the soil properties and evolution [36]. In particular, over the most recent dunes, the soils are sandy-rich and calcareous, with very low organic carbon content and scarce horizon differentiation, while in the inner dunes, the soils are more developed, with higher organic content and a better stratification. In the inter-dune lowlands terrains, the water saturation is nearly permanent, while in the inner wetlands, the soils are affected by the water regulation due to the agricultural use of the land and by the biodiversity preservation [36]. Such terrains are characterized by coarse-loamy texture and show conditions of nearly permanent saturation, since the fresh water table is close to the surface. The largest land subsidence rates in the Po Delta area have been measured during the period 1950-1957, when the Delta was subjected to an intense phase of methane extraction from the relatively-deep aquifers [37,38]. In 1960, the production wells were closed and, as a consequence, the ground lowering rates diminished. In the meanwhile, the subsidence history influenced by man-induced activity, and then suspended, shows an evident correlation with the greatest modifications of land use and cover [39]. In fact, unlike for other deltas, where the land-use practices influence subsidence processes (e.g., [40]), the Po River Delta experienced over the period 1955-1983 the greatest extension of urban fabric, industrial areas, port and sport facilities, and the largest enlargement of agricultural covering, despite the swamps and salt marshes.
Over the last few decades, the advent of space geodetic techniques for monitoring the Po River Delta and near coastal areas, such as the GNSS and DInSAR, and the integration among different tools, favored a more complete observation of the subsidence phenomenon. So far, the present day subsidence rates are found ranging from 2 to 5 mm/yr in the central part of the delta, while they increase in the eastern sector, reaching 10-12 mm/yr, corresponding to the Po River mouth, and up to 30 mm/yr in the proximity of the littoral zones (e.g., [30,41,42]).
To determine the subsidence causes, some authors studied the Plio-Pleistocene sedimentary sequence of the entire Po Basin, discriminating the natural short-and long-term components of subsidence [29,43,44]. They give particular emphasis to the role of the flexure of the Adriatic Plate beneath the Tyrrhenian Plate (tectonic process), and to the sediments-related processes, such as compaction of Plio-Pleistocene strata and sediment isostasy. In addition, these authors quantitatively described the contribution due to the isostatic rebound in relation to the Alpine glaciers load after the Last Glacial Maximum. Other researches focused their attention on the natural short-term component of the subsidence on the basis of the observed correlation between the subsidence rate and the age of the shallower sediments of the Po River Delta [30]. In particular, Teatini et al. (2011) [30] indicate the consolidation of the Late Holocene muddy sediments, which rapidly prograded eastwards in recent centuries, as the main cause of the present-day subsidence. Finally, Fiaschi et al. (2018) [42] suggest that over the Po River Delta, the recent-constructed buildings may locally accelerate the natural subsidence processes causing the settlement of the foundations.

Datasets
All data used in the present case study are open source or made available upon request by the handling institutions. The data are grouped in geodetic and hydro-meteorological datasets.

Site-Positioning Time Series (Up Component)
The geodetic data comprise three site-position time series, whose daily solutions have been produced by the Nevada Geodetic Laboratory (University of Nevada, Reno). The three continuous GPS (CGPS) stations [45] are located at the Taglio di Po (TGPO), Porto Tolle (PTO1), and Codigoro (CODI) villages ( Figure 3). They are placed, respectively, on the roof of a barn (TGPO), a three-floors (PTO1) and one-floor (CODI) buildings, whose foundations do not exceed 2 m in depth and lay on sandy and clayey sediments. TGPO and PTO1 belong to the Veneto Region GPS network, while CODI belongs to the Emilia Romagna GPS network. The Analysis Centre of Nevada Geodetic Laboratory processed the primary RINEX (Receiver INdependent EXchange format) data by using GIPSY/OASIS-II software (developed by Jet Propulsion Laboratory-NASA) [46] and obtained station 3D positions in the International Terrestrial Reference Frame 2008 [47]. The processing provides cleaned data, which consist in time series whose outliers are removed, plus information regarding earthquake and equipment events ( Figure 4). In the case of large earthquakes (Mw > 6.9) occurring near to the CGPS stations, the background trends of the American datasets are corrected during the time series analysis phase. Moreover, the series clearly exhibit non-tidal loading effects (i.e., atmospheric pressure, ocean bottom pressure, and hydrological effects). Finally, North, East, and Up CGPS components were analyzed; since the two horizontal ones displayed negligible rates after removal of the Eurasia plate movement, we considered only the vertical (up) component of ground deformation, spanning the time interval 2009-2017.

Hydro-Meteorological Data Series
The hydro-meteorological data used in this work have been collected from 20 stations located in an area of 70 × 40 km 2 extending, from Pontelagoscuro to Pila villages ( Figure 3). The data refer to three parameters: pluviometric, river hydrometric, and piezometric levels ( Table 1). In particular, the daily river hydrometric stages are collected from five stations monitored by the inter-regional agency for the Po River. The daily rainfall series, originally acquired by several institutions, such as the Air Force of Italy and Regional Environmental Protection Agency (ARPA), are available on the national system of the Italian National Institute for Environmental Protection and Research (ISPRA) and on the regional system of Emilia-Romagna ARPA. The pluviometric data come from 11 stations located in the south-central and eastern sectors of the Po River Delta area. Finally, the piezometric levels have been collected by a regional monitoring network (FaldaNET; see Table 1). They refer to four boreholes that reach the shallower unconfined aquifer (phreatic water) and the soil water (vadose water) within a few meters below the ground level (b.g.l.). Even if the groundwater data are inhomogeneous in terms of temporal sampling, the performed statistical analyses (i.e., monthly or annual means and moving averages) allowed us to define the main trends and to identify significant correlations (see Section 2.1 for details).

Results
In this section, the main results obtained by the new approach applied to the Codigoro area are illustrated following the conceptual scheme depicted in Figure 1.

Component Recognition (Step 1)
In order to quantify the permanent component of land subsidence in the Codigoro area, the geodetic data analysis has been performed not only on the CODI station, but also on the TGPO and PTO1 sites, located about 21 km and 22 km, respectively, away from CODI (see Figure 3). The permanent component has been preliminarily estimated through a linear least-square fit over a 4-year period of common time interval and inter-seismic deformation (October 2012-October 2016). The obtained linear rate is found between −1.7 and −3.4 ± 0.1 mm/yr, with a minimum in the southern part of the delta (CODI) and a maximum in the central one (TGPO, PTO1). Moreover, the geodetic velocities calculated at the TGPO (−3.4 ± 0.1 mm/yr) and PTO1 (−3.3 ± 0.1 mm/yr) stations are quite similar in the analyzed 4-year period, while PTO1 shows higher rates compared to TGPO during different data time spans [48,49]. Finally, the linear rates calculated on the CODI station, which are characterized by the longest data time span (available from January 2009), change in relation with the analyzed time interval (Table 2), thus indicating the occurrence of non-uniform processes controlling the ground motion. As far as the periodic components, they have been recognised in CODI CGPS data by using two analytical tools: the smoothed CMA and the wavelet analysis, with the latter performed by a MATLAB package software [25]. The 6-month smoothed CMA applied on the whole data set (January 2009-December 2017) reveals a narrow oscillation with a maximum amplitude of 10 mm, characterized by biennial and annual period oscillations (Figure 5a). The continuous wavelet power spectrum of CODI CGPS time series confirms the presence of these two stronger power peaks at 2.1-year (period = 770 days) and 1-year (period = 365 days) scales (Figure 5b). The thick black contour depicted in Figure 5b designates the 95% significance level against red noise, while the black contoured shading is the cone of influence where edge effects occur [25]. Significant regions are identified using the Monte Carlo method. The analyzed record also exhibits a weaker intra-annual oscillation at ã 6-month scale (period = 190 days) and scattered oscillations at a~3-month scale (period = 90 days).

Source Selection (Step 2)
Multi-disciplinary and multi-methodological comparative analyses have been performed to retrieve useful indications for getting insight into the main physical mechanisms active in the Codigoro area.
In the way to disclose the mechanism responsible for the permanent component of subsidence, the geodetic vertical velocities computed on the three CGPS stations have been compared with the thicknesses of the shallower and deeper sedimentary sequences, following a procedure previously used by other authors on the Venice coast and inland (e.g., [50]). The considered thicknesses correspond to the Plio-Pleistocene, Late Pleistocene-Holocene, and Late Holocene sequences, obtained from stratigraphic information available in the literature (e.g., [51][52][53]) and from the geological map of Italy (1:50,000) F. 187-Codigoro [35]. The comparative analysis, carried out by a linear fitting, reveals that the best correlation between velocity and thickness is obtained with the prograding sequences of the Late Holocene (see R 2 in Table 3). As suggested by other authors (e.g., [30]), this result seems to enhance the compaction process of the shallow sedimentary layers. In an attempt to fathom the causative mechanism of the periodic components detected in CODI CGPS time series, the periodicities retrieved at Step 1 are searched by analysing hydro-meteorological data both at local (rainfalls of Po River Delta) and regional scale (river hydrometric level drained from Po River catchment). In addition, the time evolution of hydro-meteorological series has been also compared with the shallow groundwater oscillations occurring within a range of a few meters of depth b.g.l., because the analyzed CGPS data refer to antennas monumented on buildings characterized by shallow foundations (up to 2 m b.g.l.), which might be affected by groundwater oscillations and relative capillary fringes [54]. Related to the four piezometers considered for the comparison, two holes are located in the clay-silty lower delta plain deposits (31FE and 35FE), while the other two wells (26FE and 23FE) reach the sandy delta front deposits of the unconfined aquifer (aquifer A0, according to Emilia-Romagna Region classification) (Figure 6a). The piezometric level of the 31FE borehole better correlates with the rain pattern of the Jolanda di Savoia station, as shown by the smoothed 6-and 9-months CMA during the 2009-2016 period (Figure 6b).
In addition to the rainfall correlation, the histogram analysis of mean annual and cumulative annual values, applied on both rainfalls and river-water level during the 2010-2017 period, reveals that the piezometric level trend correlates both with river-and rain-water patterns. Indeed, as shown in Figure 6d, the main difference between river-and rain-water trends occurs during 2015-2017, wherein the river stage decreases, meanwhile the rainfall profile increases or keeps high and quite constant values. During 2015-2017, the piezometric level of borehole 26FE, located close to the Po di Goro River, better correlates with the river level measured at Ariano SIAP hydrometric station (Figure 6c,d). Conversely, in the same period, the piezometric level of the borehole 23FE, located far from the main Po River course and the Po di Goro River, shows quite constant high values such as the rainfall pattern observed at the Goro pluviometric station (Figure 6c,d). The only exception is represented by the borehole 35FE, which presents a piezometric level that does not correlate with the rainfall trend of the near Volano pluviometric station.
In order to highlight possible common features between daily hydro-meteorological data and CODI CGPS time series, frequency analyses based on cross wavelet transform (XWT) and wavelet transform coherence (WTC) have been also performed. In particular, in time-frequency space, the XWT analysis finds regions where two simultaneous time series show high common power, while the WTC analysis finds regions where time-series covary but without having necessarily high power [25]. The analyses have been carried out between daily CGPS data of the CODI station and Po di Goro River stages measured at the Ariano SIAP hydrometric station, within the time interval from January 2009 to December 2017. The red areas in the diagram of Figure 7a reveal that the geodetic trend compared to the hydrological pattern shares strong power peaks at time scales of about 2 years, 1 year, and 6 months, while coherence is significant for the annual scale (period = 365 days) (Figure 7b). The opposition of phase indicates that when the river water level is higher, the ground level is lower and vice-versa. Meanwhile, the out of phase by 90 • between the two records indicates that when the river water level is higher, the ground level at Codigoro rises with a response delayed by~3 months. The highest cross-correlation and coherence between the daily CODI CGPS time series and the mean daily rainfall data collected at Copparo and Volano stations are found at the annual scale (period = 365 days) for both pluviometric stations (Figures 8 and 9, respectively). During the 2009-2013 and 2015-2017 periods, the CODI CGPS time series results in opposition of phase with the rainfall time series observed at Copparo (Figure 8) and Volano (Figure 9), thus indicating that when the pluviometric level is higher, the ground level is lower and vice-versa. Whilst, in 2014, the geodetic time series of CODI CGPS is out of phase of 135 • with the two rainfall datasets of Copparo ( Figure 8) and Volano (Figure 9), thus indicating that when the rainfall amount is higher, the ground level at Codigoro lowers with a response delayed by~3 months.

Source Validation (Step 3)
The results of the previous step led us to preliminarily investigate the hydrological loading mechanism depending on rainfall, soil moisture, and canopy water.

Hydrological Modeling
Continental hydrology loading is the largest contributor to surface deformation at seasonal and longer-term timescales [55,56]. To this end we estimated ground displacement due to the hydrological loading over the entire Po River Delta using physically based hydrological models. Specifically, the Modern-Era Retrospective analysis for Research and Applications (MERRA2) [57] and the Global Land Data Assimilation System (GLDAS) [58] models have been used, which allow the calculation of the elastic displacements at specified station sites [56,59].
Even if these physically-based models do not account for surface water (e.g., rivers) and deep groundwater processes, they simulate the mass transfer at the Earth's surface and other environmental parameters by integrating satellite and ground-based observations into land surface models through data assimilation techniques. They are forced by precipitation, near-surface atmospheric conditions (humidity, wind speed, pressure, etc.), as well as heat fluxes (downward shortwave and longwave radiation), and provide land surface states, including soil-moisture, snow and, eventually, canopy water (only GLDAS here) using water and energy balance equations [60].
Loading computations are performed by convolving continental water storage (soil moisture, snow, and eventually, canopy water) with appropriate elastic Green's functions [61], describing the Earth response to loading processes. Permanently ice-covered regions have been masked out, as these models are inaccurate over these areas, as ice flow dynamics is not taken into account. We also enforced the total water mass conversation by compensating any lack/excess of water over land by removing/adding a uniform layer in the oceans. Main input parameters-such as soil properties (texture, hydraulic parameters, pedo-transfer functions), bedrock, surface and root zone soil moisture thicknesses, vertical decay factor for saturated hydraulic conductivity, etc.-may be found for instance in [62].

Modeling Results
Computation of the ground displacement for all the GPS stations used in this study is available at the EOST loading service [63]. The results from MERRA2, GLDAS/Noah v1.0 and GLDAS/Noah v2.1 models slightly differ due to some modeling basic conditions. For instance, GLDAS uses the Noah v1.0/v2.1 land surface models for calculating the surface loads due to soil moisture, snow, and vegetation water content on a global mesh of 0.25 • [64] and 3-hourly samples. Conversely, the land surface model used in MERRA2 is the catchment model [60], which downscales the soil moisture variability and its effect on the runoff and the evaporation amounts. The output grid of this model has a spatial resolution of 0.5 • × 0.625 • latitude by longitude [57] and an hourly temporal resolution. Since the Po River Delta extends for an area of 40 × 40 km 2 , a single grid cell of GLDASs and MERRA2 models covers the entire delta. Except for the period 2014-2015, the annual phase and the amplitude of 6-10 mm of the elastic displacements provided by the global hydrological loading models (MERRA2, GLDAS, and GLDAS2) are in good agreement with the de-trended geodetic time series of the CODI CGPS station, as shown in Figure 10a. The linear fit preliminarily applied on original CODI CGPS time series from June 2012 to October 2016 (inter-seismic deformation time span) indicates a mean geodetic velocity of −1.9 ± 0.1 mm/yr (Figure 10b). After periodic trend removal with the best fitting model (GLDAS2), the most reliable mean subsidence rate is 2.3 ± 0.1 mm/yr (Figure 10c).

Discussion
The proposed methodological path introduces three main steps: analyses of permanent and periodic components of geodetic signals (Step 1); multi-disciplinary and multi-methodological comparative analyses, which is a key aspect to identify relevant sources (Step 2); physically-based modeling to validate the processes that cause ground deformations (Step 3).
In the Po River Delta area, Step 1 highlighted that the geodetic signals are characterized by periodicities at biennial, annual, and intra-annual (6 and 3 months) scale. These periodicities are typically found in the rainfall patterns [65] and could be related to other weather-dependent global processes, such as the Quasi Biennial Oscillation (QBO) (e.g., [66]) and the fast component of the El Niño Southern Oscillation (ENSO), the Atlantic Multidecadal Oscillation (AMO), and the North Atlantic Oscillation (NOA) [67]. The correlations found between the analyzed geodetic datasets and the hydro-meteorological data both at local (delta rainfalls) and regional (Po river hydrometric level) scale (Figure 7), allow us to preliminarily link the retrieved cyclicities to the elastic loading effects as relevant sources.
Indeed, focusing on the significant regions of the XWT and WTC plots in Figure 7 (out of the shaded zones), an overall opposition of phase in the annual signals between CODI CGPS ground displacement and river stage (time lag of 180 degree or 6 months) is observed. Such behavior preliminarily suggests a loading effect. It means that the ground level oscillation at CODI station could be triggered by changes of the stream volume of Volano River, which can be approximated by Po di Goro hydrometric level trend, and that the water river loading could affect the southern portion of Delta. The observed delay (time lag of 90 degree or 3 months) between the river stage and the ground rebound, is likely due to the presence of clayey terrains [54]. In fact, after the wet period occurred from 2013 to 2015, the fine-grained deposits could have undergone swelling effects due to water infiltration [68]. In principle, the swelling might be able to contrast and delay the up-and down-lift of the ground induced by the surface loading process. These interactions have been already unveiled around CODI CGPS site [54,69] and at TGPO and PTO1 CGPS stations, near the main Po River course [49].
In addition, the opposition of phase between the CODI CGPS time series and the rainfall mean data collected at Copparo and Volano pluviometric stations (see Figures 8 and 9, respectively), suggests the occurrence of a rain water mass-dependent process in controlling land vertical movement at the CODI site, and unveils the occurrence of an elastic loading effect involving rain water and soil. Moreover, the correlations between piezometric level and rainfall patterns ( Figure 6) reveal the occurrence of an infiltration process in the shallower part of the sedimentary coverage and confirm that soil canopy can be influenced by infiltrated rainwater supply. In particular, the abundant rainfall occurred in 2013 and 2014 ( Figure 6) could have caused the saturation of the soil canopy inducing a delay of water release in presence of clayey terrains (time lag of 135 degrees or 9 months).
Finally, the largest extent of the 95% statistical confidence areas found in the cross-correlation power and coherence index (Figure 7) indicates that the ground oscillations at CODI CGPS is mainly controlled by the rainwater loading instead of by river water loading. The ground displacement estimated over the entire Po River Delta by physically based modeling based on hydrological loading effects shows discrepancy with the observed one, especially during 2014 ( Figure 10). This difference could be due to the rainfall input, which is considered by the models as homogeneous over the entire delta because of their low resolution. Really, during the wet stage 2013-2014, some pluviometric stations located in the north-western part of the Delta show an evident peak of rainfall in 2014 [49], while other stations located in the eastern (e.g., Pradon Porto Tolle in Figure 6) and southern part (e.g., Volano in Figure 6), indicate 2013 as the wettest year. Therefore, during 2014 the models could have underestimated the actual rainwater supply.
Concluding, we remark that significant differences in the geodetic vertical velocity estimate are observed by using statistical or hydrological models [69]. In particular, it is worth noting that the mean subsidence rate of 2.3 ± 0.1 mm/yr, obtained after periodic trend removal with the best fitting model (GLDAS2), is consistent with the value 2.4 mm/yr calculated over 9-year period on the original data (January 2009-December 2017 in Table 2). The absolute values of the preliminary linear rates, calculated on datasets that include both permanent and periodic components, are smaller when time intervals get longer ( Table 2). Such evidence confirms that the longer the time span is, the lesser the effect of the periodic terms on the linear rate is [16].

Conclusions
This research introduces a methodological path aimed at dealing with the complexity of the land subsidence phenomenon. In our opinion, the choice of the relevant sources that mainly control the periodic component of land subsidence is crucial for providing accurate estimates of subsidence rates, and it needs comparative analyses with data of different nature accounting for the relevant processes related to geography and geology of the study area.
With respect to past and current studies, the developed procedure enables us to: i. overcome the one-component investigation (permanent component), thereby improving the accuracy in the geodetic velocity estimate, even when short time series span is available; ii. drive the selection of the sources by an analysis to modeling procedure, enhancing qualitative or semi-quantitative classical approaches; iii. quicken the source validation by highlighting the relevance of the source on the basis of the comparative analyses and before the modeling, differently from the peering approach, which validates the source on the basis of the model findings.
The proposed procedure has been applied to the Codigoro area, located in the southern part of the Po River Delta (Northern Italy), which is an area historically affected by land subsidence and recently interested by accurate continuous geodetic monitoring through GNSS and GPS stations. The study has focused on the analysis of three daily-CGPS time series spanning the time interval 2009-2017, several meteo-hydro-parameters collected from twenty stations and stratigraphic-geological information. The main findings suggest that a rainwater mass-dependent process might be active on the floodplain of Codigoro at the regional scale and it could explain the seasonal annual component detected in the CODI CGPS station. Since this mechanism is mass-dependent, the ground level up-or-down lift occurs according to the elastic rebound due to the soil moisture mass change, which appears to be controlled by shallow hydrological mechanisms (rain water infiltration in the topsoil, runoff, and evapo-transpiration). Moreover, the use of the most accurate physically based global hydrological models-i.e., MERRA2 and GLDAS2-enables to account for the loading effects due to rainfalls and soil moisture. Indeed, significant differences turned out by using or not the hydrological models for retrieving the subsidence velocity, which passes from 1.9 ± 0.1 mm/yr to 2.3 ± 0.1 mm/yr over the period June 2012-October 2016 at CODI station. Finally, the low subsidence rate found at CODI site compared with the higher values observed in the central part of delta (TGPO and PTO1 sites) better correlates with the thickness of the last Holocene prograding sequences. This evidence claims for the need of a future work on the modeling of the compaction process affecting this sequence. It might help to discriminate the contribution to the present-day subsidence of the Po River Delta, due to shallow sediment compaction and tectonics-driven mechanisms.
On the wake of the methodological advancement, two aspects of the new approach will be implemented through further research. Firstly, a suitable noise analysis in the component recognition phase will be introduced, in order to improve the accuracy of the subsidence rate estimates. Secondly, probabilistic modeling will be applied at the aim to better manage the uncertainty on geological data for obtaining more realistic simulations of the subsidence phenomena.