S Zerbini

F Raicich

B Richter

V Gorini

M Errico

Hydrological Signals in Height and Gravity in Northeastern Italy inferred from Principal Components Analysis

Keywords: GPS height, gravity, hydrology, subsidence, Empirical Orthogonal Functions, Singular Value Decomposition method

This work describes a study of GPS heights, gravity and hydrological time series collected by stations located in northeastern Italy. During the last 12 years, changes in the long-term behaviors of the GPS heights and gravity time series are observed. In

particular, starting in 2004-2005, a height increase is observed over the whole area. The temporal and spatial variability of these parameters has been studied as well as those of key hydrological variables, namely precipitation, hydrological balance and water table by using the Empirical Orthogonal Functions (EOF) analysis. The coupled variability between the GPS heights and the hydrological balance and precipitation data has been investigated by means of the Singular Value Decomposition (SVD) approach. Significant common patterns in the spatial and temporal variability of these parameters have been recognized. In particular, hydrology-induced variations are clearly observable starting in [2002][2003] in the southern part of the Po Plain for the longest time series, and from 2004-2005 over the whole area. These findings, obtained by means of purely mathematical approaches, are supported by sound physical interpretation suggesting that the climaterelated fluctuations in the regional/local hydrological regime are one of the main

Introduction

The long (multi year) and short (seasonal to a few years) period oscillations observed in the continuous time series of GPS heights and gravity result from the superimposition of various phenomena of different physical nature. Among them, tectonics and the mass transport in the Earth's system which induces deformation of the Earth's crust in response to variations of the load due to hydrology, air pressure and nontidal oceanic effects. Continuous GPS time series exceeding a decade are now available which can provide reliable information on the long-term height evolution (linear and nonlinear) as well as on the short-period seasonal fluctuations. However, it is complicated to unravel the contribution of the different components in the observed long-term height behavior because, at global scale, the Glacial Isostatic Adjustment (GIA) is the only coherent geological contribution to height variation for which understanding of the physical process has been achieved. Superconducting gravimeter and/or absolute gravity A c c e p t e d M a n u s c r i p t 3 measurements are not as widely spread as those of GPS, but where both types of measurements are available, they offer a unique means to detect and understand mass contributions.

It was shown that hydrological mass variations play a major role in the seasonal height and gravity variability [START_REF] Van Dam | Crustal displacements due to continental water loading[END_REF][START_REF] Zerbini | Combination of space and terrestrial geodetic techniques to monitor land subsidence: Case study, the southeastern Po Plain, Italy[END_REF]. On long time scales, climate-related variations of GPS heights and gravity have not yet been clearly identified mainly because of the limited temporal extent of most of the continuous series and for the lack of information on the spatial and temporal variability of groundwater storage.

We have studied the variations observed during the last decade in the GPS heights, gravity and the hydrological time series in northeastern Italy (Fig. 1), as well as their mutual relationships. The investigated area is affected by both natural and humaninduced subsidence. The Po Plain is a subsiding sedimentary basin encompassed by the Alps, Apennines and Dinarides Chains. Two structural environments occur within the basin: the north-verging Apennine fold-and-thrust belt system that is buried under the Plio-Quaternary cover and a platform gently dipping from the Alps into the basin [START_REF] Carminati | Subsidence rates in the Po Plain, northern Italy: the relative impact of natural and anthropogenic causation[END_REF]. Sedimentation has filled the basin with alternate stratigraphic sequences of sands, silts and clays variably interbedded and normally consolidated and containing water in the form of impregnating water and groundwater [START_REF] Zerbini | Vertical crustal motions and mean sea level: an experiment in the Eastern Po Plain[END_REF]2002). These conditions are most favorable for the development of both natural and anthropogenic subsidence. Long-term natural subsidence rates range between 1-to-2 mm/yr [START_REF] Pignone | Pianura bolognese e costa le aree più critiche[END_REF], with maximum rates of 2.5 mm/yr evaluated from the analysis of borehole stratigraphies and from available seismic sections A c c e p t e d M a n u s c r i p t 4 [START_REF] Carminati | Apennines subduction related subsidence of Venice (Italy)[END_REF] with the largest rates occurring in the southern part of the Po Plain and in the Po Delta. As discussed by [START_REF] Carminati | Separating natural and anthropogenic vertical movements in fast subsiding areas: the Po Plain (N. Italy) case[END_REF], backstripping analysis suggests that tectonics accounts for about 50% of the long-term natural subsidence, whereas compaction and sediment load account for about 30 and 20%, respectively. However, the present-day geodetic observations show subsidence rates with peaks up to 40 mm/yr resulting from an anthropogenic component superimposed to the natural subsidence. The anthropogenic contribution has been generated by systematic exploitation of groundwater for both industrial and civil use for about 50 years starting during 1950s. Groundwater control policies have been adopted since the beginning of the 1980s with the consequent reduction of fluids withdrawal.

In our study, we have adopted the Empirical Orthogonal Functions (EOF) and the Singular Value Decomposition (SVD) analyses. Significant common patterns in the spatial and temporal variability of GPS heights, gravity and hydrological parameters have been identified. In particular, hydrology-induced variations are clearly observable starting in 2002-2003 in the southern part of the Po Plain for the longest GPS time series as well as for gravity, and from 2005 for the GPS heights over the whole area. Observing and modeling long-and short-period signals in the height and gravity time series allows to better understand and quantify subsidence. This knowledge is important because of the relevant societal impacts of this phenomenon.

Data

A c c e p t e d M a n u s c r i p t 5 In northeastern Italy, several data sets are available that can be used for a spatial/temporal analysis of potentially common signals. These are time series of GPS heights, hydrological information, such as precipitation, simplified hydrological balance (precipitation minus evapotranspiration) and water table data, and series of absolute gravity and superconducting gravimeter measurements. Four of the GPS time series are rather long (about a decade or longer) and continuous.

GPS

We have used the continuous GPS data of a network located in northeastern Italy (Fig. 1 and Table 1). The observations were analyzed by means of the Bernese software package version 5.0 [START_REF] Dach | Bernese GPS Software version 5.0[END_REF]. The data were processed by using high-accuracy International GNSS Service (IGS) products including the ionospheric files, the satellite orbits and the Earth rotation parameters. The ITRF2005 coordinates and velocity field [START_REF] Altamimi | ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters[END_REF] were adopted for five IGS stations in Europe (ZIMM, GRAZ, MEDI, MATE and CAGL) used as fiducial sites in the network adjustment procedure. In Medicina, there are two GPS receivers, MEDI and MSEL, installed at a distance of about 30 m from each other. Both MEDI and MSEL are EUREF EPN stations. One important aspect is that the antennas/receivers of the stations, except for MEDI, were never substituted during the time period of the present analysis. This ensures a desirable homogeneity of the height series and the absence of jumps due to a change of instrumentation. The MEDI height time series presented in this paper for comparison purposes is a EUREF solution kindly provided by K. Szafranek and A. Kenyeres (2009).

A c c e p t e d M a n u s c r i p t 6

Hydrology

For all stations, we have estimated a simplified daily hydrological balance. The data used were the daily precipitation and maximum and minimum air temperature to derive estimates of the potential evapotranspiration according to the Hargreaves equation (Hargreaves andSamani, 1982, 1985). In addition, at Medicina, continuous measurements of the surficial water table are available.

In order to compare the hydrological balance data with those of GPS heights and gravity, we have built hydrological time series in which the n-th element is the time integral of the daily values from 1 (first epoch) to n. The linear trends have been removed from these time series in order to highlight the seasonal to interannual variability of the local hydrology. The residuals turn out to be useful proxies for the surficial water table fluctuations. Precipitation, treated in the same way as the hydrological balance, also appears to be a good proxy for the multi-year behavior of the water table. The identification of reliable proxies is necessary because water table data are not, in general, readily available. Figure 2 illustrates, as an example for the Medicina station, a comparison between hydrological balance, precipitation and water table data residuals.

The three series are in good agreement as regards the long-term behavior and the phase of the seasonal cycle. However, the amplitude of the oscillations of the water table is not always well represented by the hydrological balance and precipitation.

Gravity

A c c e p t e d M a n u s c r i p t 7 Superconducting gravimeter (SG) and absolute gravity (AG) observations are performed at Medicina, while in Loiano, Bologna, Voltabarozzo, Marghera, Cavallino and Treviso only AG measurements are carried out (Table 1). The continuous SG measurements at Medicina start at the beginning of 1998. An instrumental linear drift (18.6 nms -2 /yr) was estimated by comparison with a series of AG measurements performed by means of FG-5 instruments. The gravity data have been analyzed according to the description given in [START_REF] Zerbini | Combination of space and terrestrial geodetic techniques to monitor land subsidence: Case study, the southeastern Po Plain, Italy[END_REF].

Methods: EOF and SVD analyses

Our objective is twofold: we aim at recognizing both spatial and temporal variability patterns. We will study the behaviour of individual variables, namely GPS height, gravity, precipitation, hydrological balance and water table and of variable pairs, for example, GPS height and precipitation. The analysis is performed using a principal component approach.

In particular, following the terminology proposed by [START_REF] Björnsson | A Manual for EOF and SVD Analyses of Climate Data[END_REF] and [START_REF] Venegas | Statistical Methods for Signal Detection in Climate[END_REF], the EOF method is used to analyse individual variables and the SVD method for variable pairs. It should be mentioned that the terminology is not univocal in the literature. The EOF and SVD techniques are widely used in geophysics to analyse temporally and spatially varying fields allowing transforming the data into a different set with some desirable properties (for more details, see [START_REF] Bretherton | An intercomparison of methods for finding coupled patterns in climate data[END_REF], A c c e p t e d M a n u s c r i p t 8 von Storch and Navarra (1999), [START_REF] Björnsson | A Manual for EOF and SVD Analyses of Climate Data[END_REF], [START_REF] Venegas | Statistical Methods for Signal Detection in Climate[END_REF] and [START_REF] Hannachi | A Primer for the EOF Analysis of Climate Data[END_REF]).

The variability of geophysical fields is the result of complex nonlinear interactions between very many degrees of freedom, and a challenging task is to discriminate between signal and noise. The definition of signal and noise depends on the specific object of interest, but, in general, the signal can be a pattern in space and/or time, while the noise, physical or instrumental, consists of features that are not relevant for the signal.

Normally, the signal has larger temporal and spatial scales and fewer degrees of freedom than the noise.

The EOF analysis provides a compact description of the temporal and spatial variability of a data set of a single variable in terms of orthogonal components or, also called, statistical "modes". Conceptually, the decomposition in EOFs is analogous to, e.g., those in Fourier periodic functions or Legendre polynomials. The difference is that the EOF basis functions are not analytic but empirical; the members are not chosen on analytical considerations, but on maximization of the data projection on them. The first EOF is selected to be the pattern on which the data project most strongly. In other words, the leading EOF is the pattern most frequently realized. The second mode is the one most commonly realized under the constraint of orthogonality to the first one, the third is the most frequently realized pattern that is orthogonal to both higher modes, and so on.

Each EOF is associated to a percentage of the total variance (PVE, percent variance explained) of the original data set, which accounts for the relative importance of the corresponding mode of variability. This allows recognizing dominant modes, which are likely to compose the signal, leaving out the remainder, which is the noise.

A c c e p t e d M a n u s c r i p t Let f(x,t) be the variable of interest, defined at n locations (x 1 ,...,x n ) and at m times (t 1 ,...,t m ). The data are usually arranged in a m×n matrix (F), where each of the n columns represents the time series for a given location, and each of the m rows represents the spatial distribution at a given time. F can be written as:

F = UΓV T
Γ is a rectangular (m×n) matrix with positive or zero elements on the diagonal (γ k , k=1,….., min(m,n)) and zeroes elsewhere. U is a (m×m) quadratic matrix, whose columns (u k ) are orthogonal and represent the EOFs, (or modes or patterns, functions of space only). V T is a (n×n) quadratic matrix, whose rows (v Tk ) are also orthogonal and represent the principal components (functions of time only).

The original field can be reconstructed as a linear combination of all the modes:

F m (t) = Σ k u m k a k (t)
where a k (t) = γ k v Tk (t). Since, in the literature, the nomenclature of the quantities involved in the analysis is not univocal, we will call u k the "spatial patterns" and a k the "time components".

The SVD technique has the same purpose as the EOF, but each component describes a mode of coupled variability of two fields. Similarly to the EOF, each SVD is associated to a percentage of total covariance (SCP, squared covariance percentage).

The strength of the coupling for each SVD can be quantified by estimating the correlation coefficient between the relevant time components. In order to assess the It shall be pointed out that these analyses are purely statistical; therefore the physical interpretation of the individual components requires specific arguments.

Linear trends: GPS height, gravity and hydrology

All stations, except Loiano, which is on the uprising chain of the Apennines, are located in the Po Plain and northeastern Adriatic; this area is characterized by relevant natural subsidence. The largest rates are mainly of anthropogenic nature [START_REF] Zerbini | Combination of space and terrestrial geodetic techniques to monitor land subsidence: Case study, the southeastern Po Plain, Italy[END_REF]. Figure 3 presents the time series of the GPS heights in the ITRF2005 and the AG and SG gravity series of the stations in the network. Gravity has been multiplied by a factor (-1) in order to facilitate the visual comparison with the height data. The height time series show a decrease in the subsidence rates starting around 2004-2005; an exception is San Felice which, instead, exhibits an increase in subsidence. In Bologna (EUREF EPN station BOLG, Figure 3c), the linear trend turns out to be -10.57±0.25 mm/yr for the period 1999-2008 and -8.31±0.40 mm/yr in the time frame 2002-2008 corresponding to the period for which absolute gravity measurements are also available. The trend estimated from the gravity data is +4.17±1.75 µGal/yr, consistent with the GPS results, within the statistical error, if using the free-air conversion factor [START_REF] Vaníček | Geodesy: the concepts[END_REF]. The free-air relation mostly occurs locally [START_REF] Torge | Gravimetry[END_REF].

In the Bologna area, data acquired with the Interferometric Synthetic Aperture Radar (InSAR) technique are available. In a previous work [START_REF] Zerbini | Combination of space and terrestrial geodetic techniques to monitor land subsidence: Case study, the southeastern Po Plain, Italy[END_REF], a first comparison was made between the GPS and InSAR subsidence rates which showed the consistency of the two estimates. A recent analysis of InSAR observations [START_REF] Ferretti | Tele[END_REF], allows a further comparison supporting the reduction in the subsidence rate observed by GPS at the BOLG station. T.R.E. kindly provided us with the results of their Permanent Scatterers [START_REF] Ferretti | Permanent Scatterers in SAR interferometry[END_REF] analysis for two different periods and satellites images. The stations around the Venice area (Fig. 3 f, g, h, i and j) show reasonable agreement between the GPS and AG linear trends. However, it should be pointed out that the AG series are still limited in time and the associated errors, in particular for Marghera and Cavallino, located a few hundred meters from the sea shore, are rather large.

The main hydrological parameters used in this work are the integrated simplified hydrological balances and the precipitation series as described in section 2.2. In order to achieve an understanding of the long-period behavior of the regional hydrology, we estimated the linear trend of the hydrological balance for all stations over the period 1998-2007. It was found that, over this decade, only the northern stations at Treviso (+362±4 mm/yr), Voltabarozzo (+35±3 mm/yr) and Cavallino (+38±3 mm/yr) gain water, while in all other sites, the amount of water is decreasing. In particular, the stations located in the Plain, south of the Po River (Medicina, Marina di Ravenna, Bologna and Boretto), are characterized by large negative trends up to -500 mm/yr. The analysis of the precipitation series shows similar results.

Long-period oscillations

Many authors have shown that short-period (seasonal) height and gravity oscillations are related to seasonal mass variations (see, for example, [START_REF] Blewitt | A new global mode of Earth deformation: Seasonal cycle detected[END_REF]van [START_REF] Venegas | Statistical Methods for Signal Detection in Climate[END_REF][START_REF] Zerbini | Combination of space and terrestrial geodetic techniques to monitor land subsidence: Case study, the southeastern Po Plain, Italy[END_REF]. On the other hand, long-period (multi year) variations are not as yet well understood because of the limited length of most of the time series.

In order to identify long-period oscillations, both the GPS heights and the hydrological series were detrended to derive the residuals (Fig. 5). Possible correlations between these residuals were investigated. For the sake of clarity, we point out that we identify as hydrological balance and precipitation residuals the detrended series of the time integral values of the observed hydrological balance and precipitation respectively.

Linear regressions were calculated between the height and the corresponding hydrological balance and precipitation residual values.

By considering the Marina di Ravenna and Medicina residuals (Fig. 5b andd), over the long period (multiyear), we observe anticorrelation between the height and each of the two hydrological parameters. Anticorrelation can be interpreted as the loading effect on the Earth's crust caused by the hydrology. In particular, for the MSEL station, a regression between the hydrological balance and the height residuals over the period 1996-2008 provides a ratio of -24.7±0.7 (significant at 99% confidence level), which implies that for an increase of 24.7 mm of hydrological balance there is 1 mm decrease in height (loading effect). A similar result is obtained by using the precipitation data (-18.6±0.4, 99% confidence level). At Marina di Ravenna, on the Adriatic coast, the ratios turn out to be -7.9±0.4 and -8.1±0.2 for the hydrological balance and the precipitation residuals respectively. Since this site is on the coast, the magnitude of the loading effect turns out to be smaller than that of the inland site [START_REF] Van Dam | Atmospheric pressure loading effects on Global Positioning System coordinate determinations[END_REF], in

this case about a factor of 2-to-3 of that observed inland at Medicina. The BOLG station height residuals are characterized by a positive although weak correlation, +3.3±0.4 and A c c e p t e d M a n u s c r i p t 14 +4.1±0.2, with the hydrological balance and precipitation residuals respectively. All the northern stations exhibit positive correlation. Since correlation is not explained by loading, different physical mechanisms shall be invoked. A possible explanation could be provided by the buoyancy effect according to which an increase/decrease of the groundwater amount in the soil would produce an increase/decrease of height. Also clay swelling/shrinking on water uptake/loss may occur. The different observed behaviors can be related to the different geological settings and to the environmental conditions.

Regional analysis

The residual series of the GPS heights and gravity were analyzed both individually and in conjunction with hydrological parameters by means of the EOF and SVD approaches. We performed an EOF analysis on the individual variables, namely GPS height, precipitation and hydrological balance, in order to identify spatial and temporal variability features at local and regional scales. Additionally, at Medicina, the EOF analysis has been carried out between pairs of variables: GPS height and water table, gravity and hydrology, GPS height and gravity. The SVD analysis has been applied to detect common signals between GPS heights and precipitation as well as between GPS heights and hydrological balance. The length of the time series is an important factor because only three GPS stations are present for a decade or longer, while all the stations are simultaneously available only from 2004.

In order to reduce the high-frequency variability, the analysis is performed on weekly time series. The weekly values were obtained by averaging the original daily data. A

A c c e p t e d M a n u s c r i p t 15 weekly mean was computed when at least four daily values were available. Gaps shorter than three days were preliminarily interpolated using the Objective Analysis technique [START_REF] Gandin | Objective Analysis of Meteorological Fields[END_REF][START_REF] Bretherton | A technique for objective analysis and design of oceanographic experiments applied to MODE-73[END_REF].

The EOF and SVD analyses require that all series be defined at the same epochs; the time series were detrended over the common period. Since variables measured in different units may be involved, all the residual time series were standardized to zero mean and unit standard deviation.

In the following, several cases were considered. Since we are interested in understanding the long-period variability observed both in the height and gravity data, only the stations with observation periods longer than six years were analyzed. We point out that the gaps in the EOFs series are mostly due to missing GPS data at the Boretto station. We will discuss spatial patterns, i.e. the coefficients shown in Tables 2 through 5 for the EOF analysis and Tables 6 and7 for the SVD method, and time components, i.e.

the curves shown in Figures 6 and7 for the EOF and Figure 8 for the SVD. For each case, we will comment only the two largest EOFs (EOF1 and EOF2) and SVDs (SVD1 and SVD2) because, in general, they explain most of the observed variance. However, Tables 2 to 7 list the whole set of EOF and SVD spatial patterns.

GPS height series

Case A -Five GPS stations (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) A 6b). However, the EOF1 spatial pattern exhibits a different behaviour if we compare the northern and southern stations (Table 3): while the northern sites have winter height maxima, the southern ones, except BOLG, exhibit minima. Instead, the EOF2 spatial pattern is coherent across the whole region.

A c c e p t e d M a n u s c r i p t

Precipitation time series

The two cases described previously for the GPS heights are now considered for the precipitation time series. It shall be pointed out that, for the Medicina site, we deal with a unique precipitation series since the two GPS receivers (MEDI and MSEL) are separated by 30 m only. Moreover, for Marghera, the precipitation data are not available.

Therefore, we will analyze two cases with four and seven stations respectively. Both analyses show that EOF1 explains at least 75% of variance.

Four precipitation stations (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) The Medicina, Marina di Ravenna, Bologna and Boretto precipitation series are analyzed during the period 1999-2007. These precipitation stations are identified in the relevant Tables with the abbreviations ME, RA, BO and BR respectively. EOF1 is characterized by a coherent spatial pattern and by significant interannual fluctuations (Table 4 and Fig. 7a). Considerable precipitation minima are found in 1999, 2002, 2003 and 2007 and maxima during the winters 1999-2000, 2002-2003 to those of the previous case (Table 5 and Fig. 7b).

Hydrological Balance time series

We analyzed also the hydrological balance time series by following the same scheme used for the precipitation data described above. EOF1s, explaining about 90% variance, are characterized by spatial coherence over the entire region. The related time components exhibit a clear seasonality. The EOF2s are comparatively much smaller than those in the case of precipitation.

GPS height and precipitation

The coupled variability of GPS height and hydrological variables (precipitation and hydrological balance) is studied by using the SVD approach. Within the SVD analysis we also distinguish two cases.

Five GPS and four precipitation stations (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) We analyze the GPS heights of case A together with the Medicina, Marina di Ravenna, Bologna and Boretto precipitation series for the period 1999-2007. SVD1

A c c e p t e d M a n u s c r i p t explains 81% covariance; the precipitation spatial pattern is characterized by coherence (Table 6b). The GPS height spatial pattern is anticorrelated with precipitation at Medicina (both MSEL and MEDI) and Boretto, while RA and BOLG do not show significant correlation (Table 6a). Both time components exhibit a long-period fluctuation, with a minimum around 2003-2004 (Fig. 8a), superimposed to which are seasonal oscillations.

The anticorrelation can be interpreted as a loading/unloading effect on the Earth's crust caused by variations in the precipitation regime. The SVD2 spatial pattern of the GPS height mostly exhibits regional coherence. Both time components show a long-term fluctuation with a seasonal signal superimposed. GPS is characterized by decreasing heights until beginning of 2004 and by an increase afterwards.

The correlation coefficients between the relevant SVD time components turn out to be 0.45 for SVD1 and 0.65 for SVD2, SVD1 significant at 93% confidence level while SVD2 significance is greater than 99%.

Eight GPS and seven precipitation stations (2001-2007)

Here we analyze the GPS heights of case B together with the seven precipitation series mentioned above. The time span thus reduces to 2001-2007. The covariance explained by SVD1 amounts to about 92%. The precipitation spatial pattern is coherent (Table 7b), while the GPS height pattern is characterized by a dipole (Table 7a). As regards the northern stations, GPS heights are correlated with precipitation, while the southern ones are anti-correlated. SVD2 shows a low percentage of covariance (6%). A long-term fluctuation is clearly recognizable in both SVDs with characteristics similar to A c c e p t e d M a n u s c r i p t 20 those of the previous case (Fig. 8b). Here, the correlation coefficients are 0.45 for SVD1 and 0.64 for SVD2, the first is significant at 90% confidence level while the significance of the second is 98%.

Since SVD1 explains more than 90% of the signal and the GPS height pattern exhibits a well-defined south-north asymmetry, we further investigated the relationship between the GPS height and precipitation. After estimating the average standard deviations of the two parameters for the southern and northern parts of the region respectively, we have computed the relevant ratios indicative of the height response to precipitation variations. In the south, the ratio turns out to be 1 mm of height for 20 mm of precipitation. In the north the ratio is instead 1/35. Because, by definition, the standard deviation is a positive quantity, the signs are deduced from those of the relevant SVD1 spatial patterns (Table 7); therefore in the south the sign of the ratio is negative whereas in the north it is positive. These values are in agreement with those derived by means of linear regressions between height and precipitations residuals (see section 5.0). If we assume a 1 mm homogeneous increase of precipitation over the entire area, there would be about 0.05 mm height decrease in the southern part while a height increase of about 0.03 mm would occur in the north.

GPS height and hydrological balance

The coupled variance between GPS height and hydrological balance exhibits results quite similar to those obtained by using precipitation. In all these analyses, the SVD1s explain a slightly larger covariance compared to the same cases involving precipitation.

A c c e p t e d M a n u s c r i p t 21

The spatial patterns of the hydrological balance are always regionally coherent; however, as in section 6.4, the GPS heights of the northern stations are correlated with the hydrological balance, while the southern ones are anti-correlated except for BOLG.

Local scale analysis -Medicina

At

GPS height and water table

The two GPS height data sets, MSEL and MEDI, exhibit the same response to the water table variations (Fig. 9a andb, Table 8a andb

SG gravity and hydrology

We have compared the SG data series with those of precipitation and water table . The results provided by the two cases are similar (Fig. 10a andb; Table 9a andb), although the SG data are slightly more correlated with the water table variations (63.6% PVE) than to those of precipitation (52.7% PVE). Both EOF1s show a clear seasonal signal with summer minima and winter maxima likely resulting from the combination of two effects:

the variations in gravity due to the Newtonian mass attraction and to height fluctuations.

In fact, the EOF1s pattern coefficients have the same sign indicating on the one hand that to an increase/decrease of water mass corresponds an increase/decrease of gravity; on the other hand an increase of water/precipitation will load the crust inducing a height decrease, which turns into a gravity increase. In the case of the EOF2s, we can notice a long-period signal starting at the beginning of 2002-2003. The pattern coefficients of EOF2s have opposite sign; this mode could be explained by invoking a variation in gravity due to a height change resulting from a buoyancy effect and/or a soil consolidation process. In this latter case, a decrease of the water content in the soil will cause a height decrease which, in turn, will result into a gravity increase. The soil consolidation phenomenon is known to occur particularly in clayey soils such as that present at Medicina [START_REF] Romagnoli | Influence of soil consolidation and thermal expansion effects on height and gravity variations[END_REF]. 10a andb) suggesting that this mode could result from buoyancy and/or soil consolidation effects.

Comparisons

We compared, as an example, the original detrended GPS height series and the reconstructed signals by means of the main EOFs components for the two cases A and B (see section 6.1). According to our arbitrary choice, the reconstructed signals must explain at least 70% of the observed variance. Therefore, three EOFs were used in case A, thus accounting for 70.2% of the variance (Table 2 and Fig. 12). In case B, five EOFs were selected which represent 78.9% of the variance (Table 3 and Fig. 13). For the stations which are common to cases A and B, there are no major differences between the A c c e p t e d M a n u s c r i p t 24 two EOF reconstructions. However, it shall be noticed that in case A (Fig. 12) Marina di Ravenna shows, in the period 2004-2006 a minor discrepancy between the observed height residuals and the reconstructed signal, which does not appear in case B (Fig. 13).

This difference might be due to the fact that the reconstruction of the signal in case B is accomplished through a slightly larger variance (78.9%) with respect to A (70.2 %). In case B, the analysis is also covering a larger area with respect to A; this might imply the presence of an additional regional component.

Conclusions

Over northeastern Italy, starting around 2004-2005, a reduction of the subsidence rates has been observed in conjunction with a decrease of the hydrological load on the Earth's crust mostly due to a decrease in the amount of precipitation.

The EOF and SVD analyses of the residuals of the GPS heights, hydrological parameters and gravity made it possible to identify long and short period oscillations and also to recognize common features in different pairs of variables. These findings, obtained by means of purely mathematical approaches, are supported by sound physical interpretation that allows to point at the local/regional hydrology as one of the main contributors to the observed height oscillations. Since a percentage of variance or covariance is associated to each EOF or SVD component respectively, the mathematical decomposition of the time series, together with the physical interpretation, provides a means to attribute a relative weight to the different phenomena involved.

A c c e p t e d M a n u s c r i p t

25

In both the analyzed A and B cases (section 6.1), the time components of the main EOFs show a slope inversion starting around 2004-2005. In particular, a coherent height increase is found by EOF1 in case A (Fig. 6a, Table 2) and by EOF2 in case B (Fig. 6b, Table 3). Also the EOF1 spatial pattern in case B clearly indicates an opposite northsouth behavior of the station heights on the short as well as on the multi-year period (Table 3). For example, at seasonal scale, this means that summer height maxima of the stations in the southern Po Plain are mirrored by height minima in the northern part of the region. Case B, which involves a larger spatial domain with respect to case A, has made it possible to recognize a regional signal.

As regards precipitation, the EOF1s of both the four and seven station cases (section 6.2) show a marked decrease in the precipitation residuals starting at the end of 2005 (Fig. 7a and b and Tables 4 and5).

The SVD analysis of the GPS heights and precipitation series has demonstrated the high spatial covariance and the high time correlation between these two parameters.

The SVD1 in both analyses indicates that the GPS heights of the sites in the southern Po Plain are anticorrelated with precipitation, while the northern sites are positively correlated. We have interpreted anticorrelation as the loading/unloading effect on the Earth's crust exerted by the increase/decrease in precipitation, which occurs in the period 2002-2005 and 2006-2008 respectively (Fig. 8). The reduction in the precipitation regime observed, for example, from the end of 2005, unloads the crust and it turns into a height increase particularly at the stations located in the southern Po Plain which are characterized mostly by clayey soils.

A c c e p t e d M a n u s c r i p t

The opposite response to hydrological loading of the southern and northern parts of the region under investigation (see 6.1, case B and 6.4 "eight GPS and seven precipitation stations") results in a south-north height gradient, which turns out to be in the order of 0.1 mm over a distance of about 100 km, by assuming a 1 mm homogeneous increase of precipitation over the entire area. When observing and studying crustal deformation of a region, it is important to identify how the height of the stations changes in response to the hydrological load, in particular when weak tectonic signals are being sought.

At Medicina, a few comparisons were made, namely between GPS heights and water table, SG gravity data and hydrology and between GPS heights and SG gravity (Fig. 10). In all three cases, the EOF1s time components are characterized by a clear seasonal signal. This is due to loading/unloading of the crust in the case of the GPS height variations and to variations of Newtonian mass attraction and height in the case of gravity. GPS height and SG gravity show a marked common seasonal signal as well as a long-term fluctuation (Fig. 11) pointing to a height increase which is quite evident starting from the end of 2005. Table 1.

List of stations, North/South regional identifier, techniques available and linear trends of the GPS height, Superconducting Gravimeter (SG) and Absolute Gravity (AG) series.

The formal error associated to the estimated trends, both for GPS heights and gravity, have been multiplied by an arbitrary factor of 5 to provide conservative estimates of the errors. The superscript (a) indicates that only one AG measurement is available at LOIA.

The superscript (b) indicates that the MEDI data are EUREF weekly solutions.

Medicina

Water table (m)

Water table

Figure 2 A c c e p t e d M a n u s c r i p t

  such correlation, we need the number of independent data pairs (N ind ) involved (i.e. the number of degrees of freedom). In our case, a rigorous estimate is difficult, therefore we adopt the following definition: N ind = N tot T d / T. Here N tot is the number of available data pairs, T d the "decorrelation time", namely the first zero of the time autocovariance function of the involved variables, and T the length of the time interval covered by all the data. In the present analysis, the GPS data exhibit quite larger T d values (1.2-1.3 years) than the hydrological balance does (about 0.4 years) and precipitation (about 0.5 years). As a consequence, the number of independent data pairs is approximately reduced by a factor of 20, when the hydrological balance is involved, or by a factor of 30, in the case of precipitation.

  and gravity, the linear trends computed over the entire time span. The formal errors associated to the linear trend estimates both of GPS height and gravity have been multiplied by an arbitrary factor of 5 in order to provide a conservative estimate of the errors. Figure3dconcerns the Medicina station. It presents the MEDI and MSEL GPS height series as well as the SG and AG gravity observations multiplied by (-1). The two GPS series and those of gravity are comparable. They all show from about 2005 a different long-term behavior, a reduction of subsidence, with even a tendency to uplift during the last two years 2006-2007.

  Figure 4 (a) shows the InSAR results obtained using the ERS satelliteA c c e p t e d M a n u s c r i p t images for the period 1992-2000 and (b) the RADARSAT passes in the time frame 2003-2007. The comparison between the two images indicates a major reduction in the subsidence rate over the Bologna area in the latter period.

  , identified as case A, we analyse the GPS height time series of five stations: MSEL, MEDI, Marina di Ravenna, BOLG and Boretto. These sites are located in the southern part of the region and are characterized by the longest data records, spanning about eight years, from 1999 to the end of 2007. Concerning EOF1, the spatial pattern shows that all stations behave coherently (Table2). A height increase is observed starting 2004-2005. A clear seasonality is superimposed with summer maxima and winter minima (Fig.6a). The spatial pattern of EOF2 is characterized by opposite behaviours of Medicina (both MSEL and MEDI) and Boretto on one side, and Marina di Ravenna and BOLG on the other. The EOF2 time component also shows a long-term signal, the slope inversion starting around 2004.CaseB -Eight GPS stations (2001-2007) Here we consider the stations of the previous case plus Cavallino, San Felice and Voltabarozzo. In this case, called case B, eight stations are simultaneously available for a period of about six years, from 2001 through the end of 2007, thus allowing us to extend the analysis to the northern part of the region. Both EOF1 and EOF2 time components show a long-term change starting around[START_REF] Hannachi | A Primer for the EOF Analysis of Climate Data[END_REF] -2005 (Fig. (Fig.

  and during most of 2006. In Medicina and Bologna, EOF2 exhibits a decrease until 2003-2004, followed by an increase; the opposite is found in Marina di Ravenna and Boretto. Felice and Voltabarozzo precipitation time series are considered in addition to the four stations mentioned above; in the Tables they are identified by means of the following abbreviations: CA, SF and VO respectively. The results are quite similar

  components for MEDI and MSEL show a clear seasonal signal both in GPS heights and SG time series (Fig.11a and b) as well as a long-period feature starting at the beginning of 2005. The pattern coefficients are of opposite sign (Table10aand b), thus indicating anticorrelation between the two variables for the EOF1s. This means that the seasonal maxima and minima of GPS heights and SG gravity are opposite in phase as expected, and that, on a longer period of time, to the height increase starting in 2005 corresponds a gravity decrease. The EOF2s indicate the presence of a long-period oscillation with a minimum around 2002-2003. The pattern coefficients have the same sign (Table

  This work has been developed in the frame of the IGCP 565 project "Developing the Global Geodetic Observing System into a Monitoring System for the Global Water Cycle". The authors are grateful to Dr. Erika De Simone for having contributed to the analysis of the GPS data and to Dr. Hartmut Wziontek and Eng. Peter Wolf for the analysis of the SG gravity data. We would like to thank Prof. A. Ferretti and T.R.E. for having provided the InSAR results. We also would like to thank Drs. K. Szafranek. and A. Kenyeres for having provided the weekly EUREF GPS height solutions of the MEDI station. The authors are also grateful to two unknown reviewers the comments and suggestions of which helped us to improve the manuscript.

Figure 1 .

 1 Figure Captions

Figure 2 .

 2 Figure 2. Medicina: daily water table (cyan line), hydrological balance (red line) and

Figure 3 .

 3 Figure 3. Time series of GPS heights and gravity; they are ordered from south to north.

Figure 4 .

 4 Figure 4. InSAR maps of the terrain deformation in the Bologna area; (a) ERS

Figure 5 .

 5 Figure 5. GPS height, precipitation (black line) and hydrological balance (cyan line)

Figure 6 .

 6 Figure 6. EOF1 and EOF2 time components of GPS height residuals. Panel a) Five GPS

Figure 7 .

 7 Figure 7. EOF1 and EOF2 time components of precipitation residuals. Panel a) Four

Figure 8 .

 8 Figure 8. SVD1 and SVD2 time components of the variable pair GPS height (blue dots)

Figure 9 .

 9 Figure 9. EOF1 and EOF2 time components of the Medicina GPS height and water table

Figure 10 .

 10 Figure 10. EOF1 and EOF2 time components of the Medicina superconducting

Figure 11 .

 11 Figure 11. EOF1 and EOF2 time components of the Medicina GPS height and

Figure 12 .

 12 Figure 12. EOF reconstruction and observed height residuals. Five GPS stations: Panel a)

Figure 13 .

 13 Figure 13. EOF reconstruction and observed height residuals. Eight GPS stations: Panel

  Figure 3

  Figure 6

  Medicina, time series of various parameters are acquired continuously for more

than a decade now, GPS height

(MSEL and MEDI stations)

, SG gravity, water table data and other environmental information. This makes it possible to compare locally the different series. A few examples are provided, namely we correlate, by means of the EOF analysis, GPS height (both MEDI and MSEL) with water table, SG gravity with precipitation and water table respectively and GPS height and SG gravity. Tables 8 through 10 present the pattern coefficients while the time components are presented in Figures 9, 10 and 11.

  22height increase starting at the end of 2005. These long-period height variations are correlated with increase/decrease in the water table level as a consequence of changes in the precipitation regime.

	M a n u s c r i p t
	). Both EOF1s are clearly representative of the loading/unloading effect caused by the seasonal water table A c c e p t e d
	variations. We observe maxima in the summer when the water table decreases and
	minima in the winter when the water table rises close to the topographic surface. This
	comparison also allows identifying a height decrease starting in 2002-2003 followed by a

Table 2 .

 2 EOF analysis of five GPS height residuals series (MSEL, MEDI, RA, BOLG and BR). The second column lists the percentage of variance explained (PVE) by each EOF, the following columns describe the spatial patterns.

Table 3 .

 3 EOF analysis of eight GPS height residuals series (MSEL, MEDI, RA, BOLG, BR, CA, SF and VO). The second column lists the percentage of variance explained (PVE) by each EOF, the following columns describe the spatial patterns.

	Table 4. EOF analysis of four precipitation residuals series (ME, RA, BO and BR). The
	second column lists the percentage of variance explained (PVE) by each EOF, the
	following columns describe the spatial patterns.

Table 5 .

 5 EOF analysis of seven precipitation residuals series (ME, RA, BO, BR, CA, SF and VO). The second column lists the percentage of variance explained (PVE) by each EOF, the following columns describe the spatial patterns.

Table 6 .

 6 SVD analysis of the variable pair GPS height (a) and precipitation residuals (b).Five GPS stations (MSEL, MEDI, RA, BOLG and BR) and four precipitation stations (ME, RA, BO and BR). The two stations at Medicina, MSEL and MEDI, are coupled to the same ME precipitation data series. The second column lists the squared covariance percentage relevant to each SVD, the following columns describe the spatial patterns.

Table 7 .

 7 SVD analysis of the variable pair GPS height (a) and precipitation residuals (b).

Eight GPS stations (MSEL, MEDI, RA, BOLG, BR, CA, SF and VO) and seven precipitation stations (ME, RA, BO, BR, CA, SF and VO). The two stations at Medicina, MSEL and MEDI, are coupled to the same ME precipitation data series. The second column lists the squared covariance percentage relevant to each SVD, the following columns describe the spatial patterns.

Table 8 .

 8 EOF analysis of two time series, GPS height and water table residuals, at Medicina. The second column lists the percentage of variance explained (PVE) by each EOF, the following columns describe the spatial patterns. (a) MSEL and water table; (b) MEDI and water table.

Table 9 .

 9 EOF analysis of two time series, superconducting gravimeter (SG) and local hydrology residuals, at Medicina. The second column lists the percentage of variance

	explained (PVE) by each EOF, the following columns describe the spatial patterns. (a)
	SG and precipitation; (b) SG and water table.

Table 10 .

 10 EOF analysis of two time series, GPS height and superconducting gravimeter (SG) residuals, at Medicina. The second column lists the percentage of variance explained (PVE) by each EOF, the following columns describe the spatial patterns. (a) MEDI and SG; (b) MSEL and SG.

Table 1 .

 1

				Linear trend
	Station name	North/South Technique	(mm/yr) GPS
				(µGal/yr) SG, AG
	Loiano (a) (LOIA)	S	GPS	-1.41±0.25
			AG	---------
	Marina di Ravenna (RA)	S	GPS	-6.51±0.15
	Bologna (BOLG)	S	GPS AG	-10.57±0.25 +4.17±1.75
	Medicina (Gravity)	S	SG AG	+0.71±0.05 +0.52±0.65
	Medicina (b) (MEDI)	S	GPS	-2.44±0.25
	Medicina (MSEL)	S	GPS	-2.13±0.10
	Boretto (BR)	S	GPS	-5.00±0.10
	San Felice (SF)	N	GPS	-2.24±0.15
	Voltabarozzo (VO)	N	GPS AG	+0.36±0.15 -0.98±1.07
	Marghera (MA)	N	GPS AG	-0.88±0.40 -0.88±1.41
	Cavallino (CA)	N	GPS AG	-1.97±0.15 +0.29±1.46
	Treviso (TV)	N	GPS AG	0.00±0.25 -1.64±1.50

Acknowledgments

A c c e p t e d M a n u s c r i p t


C Prati
email: prati@elet.polimi.it

A Ferretti

D Perissin

Recent Advances on Surface Ground Deformation Measurement by means of Repeated Space-borne SAR Observations

Keywords: Synthetic Aperture Radar (SAR), INSAR, Permanent Scatterers, ground deformation

Space-borne Synthetic Aperture Radar Interferometry (INSAR) is a well known widely used remote sensing technique to get precise (sub centimetric) surface deformation measurements on large areas (thousands of km 2 ) and high spatial density of measurement points (hundreds per km 2 ). In this work the recent technological advances of this technique are presented. First, a short review of the INSAR basics is dedicated to readers who are not INSAR specialists. Then, an analysis of the improvement of ground motion measurement offered by multiple repeated space-borne SAR observations gathered by the new generation of high resolution SAR systems is given. An example obtained with the recent German TERRASAR-X system is shown and compared with the measurements obtained with the elder C-Band RADARSAT-1 system. Finally, a possible processing of multi-temporal analysis of SAR images that allow extracting ground motion information also from partially coherent targets is given. In this case the core idea is to relax the restrictive conditions imposed by the Permanent Scatterers technique. The results obtained in different test-sites show an increased spatial density of areal deformation

Introduction

Synthetic Aperture Radar (SAR) Interferometry [START_REF] Bamler | Synthetic Aperture Radar Interferometry[END_REF]) is based on the phase comparison of SAR images, gathered simultaneously or at different times with slightly different looking angles from space or airborne platforms. The phase of a SAR image contains the superposition of many terms including the distance of the radar from the illuminated targets on the ground. The phase difference of 2 SAR images gathered at different times contains a phase term proportional to the target motion occurring along the sensor-target Line-of-Sight (LOS) direction during that time interval. The measured phase difference shows an ambiguity cycle of 2π that corresponds to a 2-way travel path difference of λ (the used radar wavelength). Thus, in principle, SAR Interferometry (INSAR) has the potential to detect ground surface motion phenomena with the accuracy of a small fraction of the radar wavelength (usually from 3 to 24 cm) on large areas (thousands of km 2 ) with high spatial resolution (up to 1m with the space borne SARs of the last generation as the German TERRASAR-X and the Italian Cosmo-SKyMed -CSK). INSAR is thus a unique remote sensing tool for mapping volcano dynamics, co-A c c e p t e d M a n u s c r i p t seismic and post-seismic deformation along active faults, as well as slope instability and subsidence phenomena of different nature.

However, apart from cycle ambiguity problems, the routine use of INSAR is limited by the so called temporal and geometrical decorrelation, as well as to atmospheric artifacts.

Temporal decorrelation appears where the electromagnetic signature of the targets changes with time (e.g. in presence of dense vegetation or water basins). Geometric decorrelation is generated by the target reflectivity variation as a function of the incidence angle. Finally, atmospheric heterogeneity superimposes on each SAR acquisition an atmospheric phase screen (APS) that could have the same spatial pattern (low-pass with typically 1 km decorrelation length) and amplitude of the surface deformation phenomenon under investigation [START_REF] Beauducel | Volcano-wide fringes in ERS synthetic aperture radar interferograms of Etna (1992-1998): Deformation or tropospheric effect?[END_REF].

The Permanent Scatterers (PS) Technique [START_REF] Ferretti | Nonlinear subsidence rate estimation using Permanent Scatterers in differential SAR[END_REF][START_REF] Ferretti | Permanent Scatterers in SAR Interferometry[END_REF][START_REF] Kampes | Radar Interferometry: Persistent Scatterer Technique[END_REF][START_REF] Hooper | A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers[END_REF][START_REF] Berardino | A new algorithm for surface deformation monitoring based on small baseline differential interferograms[END_REF] developed in late nineties, is a powerful tool that overcomes the above said limitations of INSAR by simultaneously exploiting all the available SAR images gathered during repeated satellite passes (e.g. exploiting ESA, CSA, TERRASAR-X and, once available, the CSK archives).

After a short introduction to the PS Technique, recent examples of ground motion measurements showing the actual potential and accuracy of this remote sensing technique are presented in section 2. Also in this section the very first analysis of the improvement that can be achieved passing from the 20x5sqm resolution cell of C-Band SAR systems to the recently available 3x3sqm resolution X-band SARs is presented. All these satellite SARs observe the same area at different times from slightly different looking angles since repeated orbits do not follow exactly the same track (due to the atmospheric drag, solar wind, etc). The distance between two orbits is called interferometer baseline and its projection perpendicular to the line of sight (LOS) is called perpendicular baseline (see Figure 1).

Figure 1

A c c e p t e d M a n u s c r i p t

The phase difference between two SAR images is called interferometric phase. The following expression of the interferometric phase with respect to a selected ground reference point holds (Ferretti et al. 2007):

γ φ φ λ π θ λ π θ λ π φ ∆ + ∆ + + + = ∆ APS n n d s R B h R B 4 tan 4 sin 4 [Eq. 1]
In Eq.1, h is the relative ground elevation of targets referred to a horizontal reference plane, s is the relative slant range position of targets, d is the projection along the LOS of the relative displacement of targets, n B is the perpendicular baseline, R is the SARtarget distance, θ is the "off-nadir" angle of the LOS, The first two phase terms in Eq.1 can be computed and eliminated if a precise DEM is available together with precise orbital information (baseline, sensor-target distance and off-nadir angle). The residual phase component, the so called Differential SAR Interferometric (DINSAR) phase, is then proportional to the terrain motion component along the LOS plus atmospheric and decorrelation noise: . From Eq.2 it is also clear that the actual possibility of measuring the terrain motion with a precision of a small fraction of the radar wavelength is mainly limited by the atmospheric and decorrelation noise terms.

γ φ φ λ π φ ∆ + ∆ + = ∆ APS d d 4 [Eq. 2]
Since the incidence angle of the radar signal is usually less than 30 deg (for ERS the nominal figure is 23 deg), SAR data are usually very sensitive to vertical displacements.

In Figure 2 the INSAR phases showing the co-seismic deformation of the Bam (Iran) earthquake that occurred on Dec. 26 th , 2003 are shown. In this case the terrain motion is measured with a sub-centimetre precision all over the imaged area since the decorrelation phase noise terms are really small (due to lack of vegetation and small baseline) as well as the atmospheric components (due to the particular meteo conditions).

Figure 2

In many other situations, however, both decorrelation and atmospheric phase noise are so strong to severely affect terrain motion measurement. A typical example is illustrated in However, on the one year interferogram (a time span that would allow to identify a big landslide motion) also the in land area is strongly affected by phase decorrelation thus preventing ground motion measurements.

A c c e p t e d M a n u s c r i p t It is interesting to note, however, that many targets that are not affected by decorrelation noise could be present in areas showing an apparent complete decorrelation as on the right image of Figure 3, but they cannot be identified since surrounded by noisy phases.

The simple visual experiment illustrated in Figure 4 shows that up to 10% of randomly distributed perfectly phase coherent pixels cannot be identified when surrounded by totally random phase pixels.

Figure 4

In general, isolated phase coherent targets can be identified if more SAR scenes gathered over the area of interest are jointly analyzed. Moreover, creating many interferograms, it is usually possible to highlight problems related to spurious phase components due to DEM errors, noise and/or atmospheric signals and try to separate the different contributions. This is exactly the idea which underlies the "PS technique".

Moving from INSAR to PS

As long as a significant number and density of independent radar-phase stable points (i.e., Permanent Scatterers) exist within a radar scene and enough radar acquisitions have been collected, displacement time series and range-change rates can be calculated. Using the PS method, we can resolve surface motions at a level of ~0.5 mm/yr at about 1km from the reference point and can resolve very small-scale features, including motions of individual targets/structures (e.g. a bridge or a dam), not previously recognized in traditional SAR interferometry (Ferretti 2007[START_REF] Adam | Practical persistent scatterer processing validation in the course of the Terrafirma project[END_REF].

A c c e p t e d M a n u s c r i p t

The PS approach is based on a few basic observations. There are ground targets that maintain a coherent reflectivity to the radar in time even when observed from different looking angles (the Permanent Scatterers -PS). The interferometric phase in correspondence of these targets is not randomized by temporal and geometric decorrelation phenomena. The PSs can be identified by a statistical analysis of amplitudes and phases of a series of SAR images gathered at different times from slightly different looking angles (the larger the number of images the more reliable is the statistical analysis). Atmospheric phase artefacts show a strong spatial correlation within every single SAR acquisition, but they are uncorrelated in time. Conversely, the phase term due to the target motion is usually strongly correlated in time and can exhibit different degrees of spatial correlation depending on the phenomenon at hand (e.g. subsidence due to water pumping, fault displacements, localized sliding areas, collapsing buildings, etc.).

Finally the topographic phase term shows a linear behaviour with the normal baseline (as seen from Eq.1).

Thus, the phase terms coming from motion, atmosphere and topography can be estimated and separated from the others in correspondence of the PS locations by looking for their different behaviour in space and time.

Due to the high spatial correlation of the APS, even a low density of PS (3-4 PS/km 2 ) allows proper sampling of the atmospheric components. On each PS sub-metre elevation accuracy (due to the wide dispersion of the incidence angles available, usually 70 millidegrees with respect to the reference orbit) and millimetric terrain motion detection (due to the high phase coherence of these scatterers) can be achieved, once APS's are estimated and removed. In particular, relative target LOS velocity can be estimated with A c c e p t e d M a n u s c r i p t accuracy often better than 0.5 mm/yr, due to the long time span. The higher the accuracy of the measurements, the more reliable the differentiations between models of the deformation process under study, a key issue for landslide monitoring and risk assessment.

The final results of this multi-interferogram approach are the following:

• a map of the PS identified in the image and their coordinates;

• their average LOS velocity (in mm/yr);

• the estimated motion component of each PS as a function of time (in mm)

• the estimated Atmospheric Phase Delay along the LOS for every acquisition date

PS Data Precision

In the PS technique all the interferograms used in the processing are generated with respect to a single master scene, properly selected among all the available SAR acquisitions so that geometrical and temporal decorrelation phenomena are kept to a minimum. This strategy does not allow the exploitation of radar targets exhibiting PS behaviour only on a sub-set of images (as it will be discussed in Section 3), however, it allows:

i -the exploitation of all the available SAR images (thus maximizing the temporal sampling rate of the PS motion time series);

ii -the exploitation of large baselines (thus getting precise 3D PS positioning that is essential for estimating the atmospheric delay)

iii -the direct sampling of the actual PS motion time series.

A c c e p t e d M a n u s c r i p t

Moreover, it allows a first-order precision assessment of PS results using standard procedures, e.g. the computation of the a posteriori standard deviations of the PS velocity and elevation values. In Table 1 the typical precision of PS results obtained processing C-band data acquired by a multi-year ESA-ERS data-set is reported. It should be pointed out that new high resolution X band sensors can achieve actually better performance in both geocoding accuracy and sensitivity to target motion, nevertheless figures are noticeable and justify why InSAR data are gaining importance in space geodetic applications [START_REF] Hilley | Dynamics of Slow-Moving Landslides from Permanent Scatterer Analysis[END_REF][START_REF] Dixon | New Orleans Subsidence and Relation to Flooding after Hurricane Katrina, as measured by Space Geodesy[END_REF], Ferretti et al. 2004[START_REF] Berardino | A new algorithm for surface deformation monitoring based on small baseline differential interferograms[END_REF][START_REF] Kampes | Radar Interferometry: Persistent Scatterer Technique[END_REF][START_REF] Hooper | A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers[END_REF], Ferretti 2007) Table 1 Although a thorough statistical analysis of PS measurements is beyond the scope of this paper, it is worth pointing out that the characterization of the precision of a PS velocity field using just one statistical figure (e.g. the maximum standard deviation value within the area of interest) can be misleading. Since the precision strongly depends on the distance between the PS under study and the reference radar benchmark used in the analysis, it is actually the variogram (or the autocorrelation function) of the phase components that can provide a more comprehensive statistical characterization of the results.

As a rule of thumb, in any DInSAR analysis, low-frequency spatial components of the deformation pattern (such as glacial rebound effects) are always more difficult to estimate than local features (such as a landslide). In fact, both atmospheric effects and orbital fringes, i.e. the main noise sources in PSInSAR analysis apart from phase unwrapping If no prior information is available for the area under study, the use of many data-sets coming from different acquisition geometries (ascending, descending orbits, parallel tracks, etc.) or even multi-platform analyses can help clarifying the scenario and improving the quality of the results.

The effect of resolution improvement

The ESA-ERS archive, historically, was the first one to be exploited for PSInSAR A c c e p t e d M a n u s c r i p t

Since a Permanent Scatterer usually corresponds to a scattering center dominating the radar return of a resolution cell, the probability that a pixel of a SAR image exhibits PS behavior is certainly related to the spatial resolution of the SAR data used in the analysis.

In fact, the higher the resolution, the higher the probability that a certain scattering center Apart from atmospheric effects, severe limitations to DInSAR applications can arise from possible phase unwrapping errors. Indeed, phase aliasing can occur more likely than at Cband and X-band results should be carefully checked for data consistency, at least in areas exhibiting low coherence values.

In order to assess the theoretical analysis summarized in the previous paragraphs, two data-stacks acquired over the same area by both RADARSAT-1 and TERRASAR-X (TSX) have been used (Table 2). The area is in the Southern part of Italy. (the ratio between the two wavelengths) at X-band, the dispersion of the elevation values is higher than expected. This is almost certainly due to the fact that the PS pairs used in the comparison could belong to different parts of a building or even to two distinct structures. Moreover any geo-coding error would impact more strongly the dispersion of the elevation values rather than the average displacement rates, since the local velocity field turned out to be smooth.

It is worth mentioning that we haven't experience any particular problem in unwrapping the phase values of the X-band datastack. This is probably due to the fact that the area is not affected by severe displacement phenomena and the local topography is not rough.

Finally in Figure 7 a vertical motion time series (LOS measurements have been projected

on the vertical axis) obtained by combining RADARSAT-1 and TERRASAR-X data is

shown. Here a similar vertical velocity estimate dispersion can be noted together with the higher temporal sampling offered by TERRSAR-X with respect to RADARSAT-1.

Multi-master interferograms and targets life-time

In the classical PS analysis the interferometric phase is generated by referring all images to a common Master acquisition. In the normal baseline -temporal baseline space this configuration can be represented with a star graph as in Figure 8, where each point indicates one of the 84 ERS images of the area of Dossena in the Italian Alps and each connection an interferogram.

Figure 8

In this framework, the actual coherent scatterer motion is sampled at the acquisition time of all the available SAR images. Of course only targets that coherently interfere with the Master acquisition can be exploited to get meaningful measurements. These targets should be then almost point-wise dominant targets within the resolution cell (in order to maintain coherence for large normal baselines) and temporally stable along the whole acquisition time interval of the SAR images (years, in most cases). Alternative SAR images combinations have been proposed in literature trying to relax the conditions under which a target acts as a Permanent Scatterer. As an example in the Small Baseline approach (SBAS) [START_REF] Berardino | A new algorithm for surface deformation monitoring based on small baseline differential interferograms[END_REF], the entire set of available SAR images is this approach some of the available images may not be used, the resulting graph is often disconnected thus preventing the correct motion measurement without an available a priori model and the small baseline condition does not allow a precise target 3D location in space thus preventing a good separation of elevation, motion and atmospheric phase components.

What is proposed here is to search for the minimum best coherence graph connecting all the images of the data-set, without imposing any pre-defined decorrelation model (temporal or spatial), simply by maximizing the interferometric coherence. In second instance, the minimum graph can also be increased by adding more connections to make the estimate more robust. Figure 9 shows the obtained graphs for the same 84 SAR images data-set used for generating the "star graph" of Figure 8. As visible in the picture, many selected connections have small normal baselines, but not all interferograms follow this rule. Moreover, the graph is connected and can change from pixel to pixel.

Figure 9

During the generation of the time-baseline connection graph, the coherence of all the ( )

2 1 N N -
interferograms that can be generated with N SAR images should be analyzed.

From the analysis of these coherence maps it is possible to generate a map of the average coherent life-time of the imaged scatterers. The map is obtained by averaging the interferometric temporal baselines weighted by the corresponding interferometric coherence. By using data of the ESA-ERS archive that include also the so called "tandem" interferometric pairs gathered with just one day time interval, the measurable to red (more than 1000 days). The red (long life time) scatterers correspond to urban and exposed rocks areas (i.e. they are good PS candidates) whereas blue scatterers remain stable just for one day and cannot be exploited for ground motion estimation (noncoherent scatterers are not colored at all). In between there are very many scatterers with a medium lifetime shorter than that of PSs, but long enough (a few months) to allow ground motion measurement. Moreover, the "birth" and "death" dates of such medium life-time scatterers belonging to a small area affected by the same surface motion are different and by combining their individual life interval, the measurable time interval of the areal ground motion can be extended up to several years.

Figure 10

The core idea of the new technique is to jointly estimate elevation and deformation trend of a given point from the complex interferometric spatial coherence. In this way, the absolute value of the spatial coherence acts as a weight in the estimation process and only coherent interferograms play a role in the retrieval of the target parameters (uncoherent ones are weighed out). Note that the set of the actually used coherent interferograms can be different from point to point. Moreover, the phase used in the estimate is filtered by the coherence estimation process. This is an advantage for distributed targets but can lead to errors in case of high density point-wise targets areas whose elevation and motion are averaged. It has also been shown that the theoretically less precise areal measurements on a higher set of partially coherent targets can complement the more precise point wise scatterers time series achievable with the PS Technique.

The integration of these techniques with other independent measurements (e.g. GPS)

represents a powerful tool that can be exploited by geologists and geophysicists for investigating many ground deformation phenomena (Burgmann 2006).

Finally it should be pointed out that maps of the differential travel time delay are generated as a side product of these techniques. These maps could be usefully exploited to get information on the troposphere water vapor content with unprecedented spatial resolution.

A c c e p t e d M a n u s c r i p t

Finally

  , in the third section the results obtained by means of an INSAR data processing approach obtained by relaxing the strict conditions imposed by the PS technique are A c c e p t e d M a n u s c r i p t shown together with an analysis of what it is gained and what is lost with respect to the PS Technique. 2.0 INSAR and DINSAR (Differential SAR Interferometry) The European Space Agency (ESA) satellites ERS-1, ERS-2 and ENVISAT are collecting SAR data at C band (5.6 cm wavelength) since late 1991 along near polar orbits with a revisiting time of 35 days. Since 1995 the Canadian Space Agency (CSA) satellites RADARSAT-1 and RADARSAT-2 are collecting similar data with a revisiting time of 24 days. The Japanese Space Agency (JAXA) satellites J-ERS (till 1998) and ALOS have been collecting L band (23.6 cm wavelength) SAR images with a revisiting time of 44 and 46 days respectively. More recently, the German satellite TERRASAR-X and the Italian constellation of satellites Cosmo-SkyMed (CSK) started collecting SAR images at X band (3 cm wavelength) with a revisiting time of 11 and up to 4 days (with 4 satellites in orbit) respectively. An impressively huge archive of repeated SAR images is thus available for imaging and INSAR applications covering the entire Earth surface.

  on temporal and geometric decorrelation of targets.

  The sensitivity of the DINSAR phase to terrain motion depends on the SAR wavelength λ since

Figure 3

 3 Figure3where the INSAR phases generated with ERS images taken in one day and one

Figure 3

 3 Figure 3

  A c c e p t e d M a n u s c r i p t errors, are characterized by high spatial correlation, since their correlation length typically exceeds 1 km. Locally, spurious phase components are compensated for by the double difference computation inherent in any DInSAR analysis, but regional signals affecting hundreds or even thousands of square kilometres can be difficult to discriminate without prior information. That is the reason why PSInSAR and GPS are complementary data that can be used in synergy to map surface deformation phenomena. Indeed, geodetic observations from a few permanent GPS stations within the area of interest (a satellite frame is typically 100x100 km 2 wide) can be used successfully to calibrate PS data (all PS data a referred to a PS supposed motionless) and remove systematic errors in SAR observations (e.g. orbital components).

  analyses. The two ERS platforms followed by ENVISAT, operating at C-band and sharing a 35-day repeat cycle, created an historical archive of SAR acquisitions since 1991, when the ERS-1 platform was successfully launched. Since then, an ever growing archive of SAR data with a spatial resolution of about 20x5 m became available and made it possible the development of many algorithms and techniques to exploit phase information (ESA archives can be found at http://earth.esa.int/resources/catalogues/) A c c e p t e d M a n u s c r i p t RADARSAT-1 archives, launched in 1995 and still operational, offer another source of InSAR data at C-band with a 24-day repeat-cycle, but the availability of many acquisition modes for the radar sensor has reduced significantly the average number of homogeneous data-stacks for interferometric applications, although in some areas (e.g. Italy) it has been carefully planned to acquire interferometric data-sets since 2003. In January 2007, RADARSAT-2 was successfully launched: this should guarantee the continuation of all the applications triggered by RADARSAT-1 data availability and allows the acquisition of high resolution (<3 m) data at C-band.Common to optical sensors, the general trend of the aerospace industry is towards SAR sensors featuring an ever increasing spatial resolution and shorter temporal sampling.During the last few years, four new SAR satellites operating at X-band have been launched, all featuring high spatial resolution and short repeat cycles: three belong to the dual-use Cosmo-SKymed constellation operated by the Italian Space Agency, and one is the German TerraSAR-X. It is then reasonable to ask what will be the impact of highresolution data on PSInSAR analyses.Theoretically, the X-band wavelength (λ=3 cm) is equal to about half of the C-band (the band of the ERS, ENVISAT and RADARSAT satellites) resulting in a greater sensitivity (approximately double) to movement, that can lead to more accurate displacement measurements. Moreover, the revisiting time is much shorter (11 days for TerraSAR-X and up to 4 days for the CSK constellation, rather than 24 or 35 days) and data are characterized by higher spatial resolution (<3 m), making it capable of highlighting details never-before visible in satellite C-band images.

  can dominate the surrounding scatterers sharing the same resolution cell. Of course, Xband SAR interferograms over non-urban areas are more prone to temporal decorrelation phenomena than C-band or L-band data, due to the increased sensitivity of phase values to any change in scatterers distribution, but, in principle, over man-made targets or rocky areas the spatial density of PS should significantly increase with respect to lower resolution datasets. This should be expected also due to the improved orbital stability of the newer platforms, allowing lower baseline values compared to the satellite SAR sensors launched in the nineties and the shorter repeat-cycle, limiting temporal decorrelation effects, at least for interferogram generated using two successive acquisitions.

  Here 22 TSX images acquired in Strip-Map mode from April 2008 up to January 2009 (11 days revisiting time) have been processed and compared to a previous PSInSAR result A c c e p t e d M a n u s c r i p t obtained from 22 RADARSAT-1 S3 (Standard Beam) data gathered from December 2006 to September 2008 (24 days revisiting time). In both cases, all PS exhibiting atemporal coherence exceeding 0.85 have been selected. This threshold has been set to limit the number of possible statistical outliers generated by random noise to less than 10 - 5[START_REF] Colesanti | SAR Monitoring of progressive and seasonal ground deformation using the permanent scatterers technique[END_REF] and corresponds, in the hypothesis of phase gaussian noise, to a standard deviation of about 0.6 rad of the phase values.Two snapshots of the results are reported in Figures5 and 6. In general, we measured an increased PS density by a factor of 12 that can open brand new applications to PS data, such as those related to the monitoring of individual buildings and structures. In particular, looking into Figure6, it should be noted how, passing to high resolution SAR data, the highway crossing the area from North to South, can now be accurately monitored by means of PS data.Since no prior information was available neither on the local topography (apart from the SRTM DEM used for the generation of the differential interferograms), nor on the surface displacement field affecting the area of interest, it was not possible to carry out a quantitative precision assessment. Moreover, it should be noted that the two datasets span two different time windows, so -at least in principle -the area could be affected by different displacement fields in the two PS results. However, the two PS datasets have been compared by computing the difference of both average displacement rates and elevation values of a set of PS pairs (one from the RADARSAT dataset and one from the TSX dataset) less than 10 m apart. The standard deviation of the difference between the estimated displacement rate was 4.6 mm/yr while for elevation values we obtained a dispersion of 5.6 m. While the first figure is in A c c e p t e d M a n u s c r i p t agreement with what should be expected based on the time spans of the datasets and the impact of atmospheric components having a dispersion of 1 rad at C-band and 1.8 rad

  Figure 5

  of images showing a very small normal baseline thus relaxing one of the target conditions requested by the Permanent Scatterers technique. As drawbacks of

  time ranges between one day and some years (e.g. 10 years for the data set used in this example). In Figure10such a scatterer's life-time map is shown referred to the area of Dossena. The logarithmic colour scale of Figure 10 ranges from blue (1 day)

  As an example, the average ground areal velocity of the area of Dossena measured by exploiting the above said multi-master interferograms and partially coherent targets strategy is shown in Figure11. As a comparison the single target average ground velocity of the same area measured by means of the standard PS Technique in shown in figure12.By comparing Figures11 and 12, it is quite evident that the number of measurement points increases by more than one order of magnitude passing from PSs to partially coherent targets. A detailed comparison of the estimated ground velocities at the PSs locations has shown that they are often in good agreement (standard deviation around .5 mm/year), but for some cases where two or more neighbouring PSs with different velocities have been averaged during the coherence estimation step of the partially coherent targets processing.

  Figure11

Figure 1 -

 1 Figure 1 -INSAR acquisition geometry.

Figure 2 -

 2 Figure 2 -Example of co-seismic interferogram of the Bam (Iran) earthquake, Dec. 26 th ,

Figure 3 -

 3 Figure 3 -Example of interferograms of the town of Ancona (on the Adriatic Italian

Figure 4 -

 4 Figure 4 -Left: random spatial distribution of phase coherent samples (10% of the total)

Figure 5 -

 5 Figure 5 -Spatial distribution of the Permanent Scatterers identified in an urban area

Figure 6 -

 6 Figure 6 -Average Spatial distribution of the Permanent Scatterers identified in a semi-

Figure 7 -

 7 Figure 7 -Vertical velocity time series of a PS obtained by combining RADARSAT and

Figure 8 -

 8 Figure 8 -Normal baseline -temporal baseline space of 84 ERS SAR images of the area

Figure 9 -

 9 Figure 9 -Normal baseline -temporal baseline space of the same 84 ERS SAR images

Figure 10 -

 10 Figure 10 -Scatterers life-time map of an area of the Italian Alps. The logarithmic colour

Figure 11 -

 11 Figure 11average ground areal velocity of the area of Dossena measured by exploiting

Figure 12 -

 12 Figure 12single target average ground velocity of the area of Dossena measured by

Figure

  Figure

Figure

  Figure

Table 2 :

 2 Nominal parameters values of RADARSAT-1 and TERRASAR-X data used in the comparison analyzed in the text.