Estimation of Subsidence in Po Delta Area ( Northern Italy ) by Integration of GPS Data , High-Precision Leveling and Archival Orthometric Elevations

Subsidence in a deformation area can be measured in various ways, examples being conventional high-precision leveling, differential InSAR and multi-temporal GPS surveys. Integration of methods can improve results, and is crucial to extract high-precision data. In particular, orthometric and ellipsoid elevations, surveyed at different moments in time, can be compared to yield information on vertical movements when geoid anomalies are known. However, a data checking procedure must be applied if archival orthometric elevations are used, because long-term measurements for many historical benchmarks may have been lost and/or replaced with other points, but at different elevations. This type of checking can be carried out over an area without gravimetric anomalies by modeling geoid undulations and vertical displacements in the time-span used for analysis, excluding points with anomalous values. This procedure was tested and applied in the Po Delta area (northern Italy), historically subject to high subsidence rates: the leveling benchmarks of 1983 were measured with the GPS technique in 2008. After checking of archival data and transformation from ellipsoid to orthometric elevations, comparisons of the same points and interpolations on the study area provided a subsidence map for the 1983-2008 period.


Introduction
Ground deformation in areas subject to vertical displacements can be studied by comparing the orthometric elevations of benchmark leveling measured at different times.This technique supplies high-precision data, but the measurements require long timings and are expensive to monitor vertical deformations.In the past, this method was the only one available, yielding high-precision data for assessment of land subsidence along the northern Adriatic coastland.
It was only in the 1990s that GPS (Global Positioning System) methods, both continuous (CGPS) and differrential (DGPS), began to be used for measuring sufficiently accurate ground movements, with the result that many permanent GPS stations, which acquire data continuously, are now available [1].Many GPS points have also been surveyed on ground which, when measured subsequently, yield the same differences: in this case, the precision with which GPS elevation values are determined must be estimated with care.Again since the 1990s, space-borne observation techniques based on Synthetic Aperture Radar (SAR) interferometry have also been used to detect subsidence and uplift along the Veneto coastland.First differential InSAR (DInSAR) and later Interferometric Point Target Analysis (IPTA) have been applied to measure and map displacements occurring since 1992, when SAR data first became available [2].Both GPS and SAR data can now provide analyses covering the last twenty years, but integration with conventional high-precision leveling data is not simple, due to differences in reference systems.Ellipsoid and orthometric elevations can be compared when the data refer to the same surface, when the geoid heights of the analysed area are known, and ellipsoid (GPS) elevations can now be transformed into orthometric ones [3] [4].GPS elevations applied to geoid heights, for direct comparisons between orthometric (conventional leveling) and transformed orthometric (GPS) elevations, performed at different times on the same benchmark, yield data on vertical displacements in the time-span in question.Thus, the congruence of elevation of the same benchmark can be checked over long periods of time.Archival leveling benchmarks, measured in the past but not repeated recently, must be verified, to ascertain that GPS measurements are performed on the same archival benchmark: only in this case can the two methods be integrated.Other techniques, such as aerial digital photogrammetry and laser scanning, are useful for measuring vertical movements by comparisons between DTMs (Digital Terrain Models) acquired on the deformation area at different moments in time.Photographic archives, which are available for the last 70 -80 years, are also an important source of historical data, enabling the ground surface to be modeled on a local scale over some decades.Historical images can be processed by archival photogrammetry, yielding multi-temporal metric information of the study area [5]- [10].However, these methods are only useful when the expected displacements are of the order of 50 cm or more, due to problems of data resolution and technical accuracy [11]- [13].
In this work, comparisons between orthometric and ellipsoid elevations are used to evaluate the subsidence of the Po Delta area in the Veneto Region (northern Italy; Figure 1) for the period 1983-2008, with expected deformations of a few tens of centimeters.In this area, subsidence was monitored by GPS in 2008 on archival benchmarks and compared with elevations measured in 1983.In addition, due to the long time-span (25 years), a procedure to verify the 1983 orthometric elevations was performed, because some benchmarks may have been lost and replaced with other points at non-homogeneous elevations.Subsequently, these data were integrated with the orthometric elevations of conventional leveling lines measured in 1984 and 2005.

State-of-the-Art of Subsidence in Northern Italy and Po Delta
Long-term natural subsidence is mainly the result of deep tectonics, glacial cycles and geodynamic movements [14]- [19].Carminati et al. [14] indicated the downflexure of the Adria plate as one cause of movement, at a maximum rate of 2.5 mm/year.The possible effects of tectonic activity along the northern Apennine thrust front buried under quaternary sediments was thoroughly discussed in [20].More sophisticated measurements now show that the recent deformation pattern of the Apennine belt and Po Plain is due to the northward movement of the Adria plate [21].
Subsidence in the study area has been periodically monitored by leveling measurements since 1897 [22], and comparisons between the leveling lines for the period 1897-1957 indicate subsidence rates of a few millimeters per year.Total subsidence dramatically increased in the second half of the 20 th century, when economic development in the region almost quadrupled the amount of groundwater extracted for industry, agriculture and domestic purposes, compared with previous years [23] [24].The most recent dataset, collected by the regional authorities, indicates that subsidence in the last years of the 20 th century was in some areas at least one order of magnitude higher than that produced by natural geodynamic and geological processes.
Along the Veneto coastland, human activities-pumping water and gas from underground reservoirs-still creates many problems in inhabited areas only a few meters above sea level, since they are continuously subject to frequent flooding [25] [26].One clear-cut example is the ancient city of Venice, which lies in a barrier island lagoon system just north of the Po Delta and is regularly subject to flooding by high tides.Both natural and anthropogenic land subsidence and eustasy in the northern Adriatic have caused 23 cm of relative land subsidence, with respect to the mean sea level over the last 100 years [15]- [18] [27]- [29].Spirit leveling and CGPS and DInSAR methods, applied since 2000, have shown that subsidence no longer occurs in the center of Venice and in its industrial zone, but is still ongoing in the northern and southern lagoons and bordering land [30] [31].Baldi et al. [32], using CGPS method, provides data on maximum vertical displacements of some mm/year at Modena, Reggio Emilia, Bologna, Ravenna and the Po Delta.
Deltas and their estuaries have great ecological and economic value and support both high populations and extensive agriculture [33].The influence of the sea (waves, storm surges, flooding) and the reduction of river sediment supplies, due to natural processes and human interventions, make the most important Mediterranean deltas (Nile, Rhone, Po and Ebro) highly vulnerable to rises in sea level [34].
The Po Delta covers an area of about 400 km 2 , extends seaward for about 25 km (Figure 1) [35], and is characterized by a complex multi-aquifer freshwater system [36].Subsidence is an important factor in the altimetry of the Po Delta.It has both natural and artificial origins: natural, caused by soil lowering due to the compaction of lithological layers of "young" soil and oxidation of peat, and artificial, due to draining of wetlands, land reclamation and, in particular, the extraction of methane water, which was carried out with the greatest intensity from the 1940s until 1961, when the Italian government suspended extraction in the area [37].Natural subsi-dence has been estimated at 1 -2 mm/year in the Delta.Reclamation has also caused local subsidence, but cannot completely explain the falling rates over the last 60 years [38]- [44].
Due to the significant economic impact of subsidence, vertical movements have been monitored since 1950 [45] by repeated high-precision leveling, undertaken by the IGMI (Istituto Geografico Militare Italiano) and later by other local authorities, companies and institutions [22] [39] [44] [46]- [49].Caputo et al. [39] and Borgia et al. [47] reported a maximum subsidence rate of 250 mm/year in the central part of the Po Delta for the period 1951-1957 and 180 mm/year between 1958 and 1962.Later, from 1962 to 1967, these rates fell to 33 mm/year, matching the progressive reduction of extraction [39] and to 37.5 mm/year from 1967 to 1974 [40].These last data clearly show the benefits gained by halting extraction, begun in the late 1970s.In subsequent years, the subsidence rate decreased, although the actual observed subsidence in some areas is at least one order of magnitude higher than that due to natural processes.The GPS technique has more recently been used to integrate conventional leveling with episodic and continuous measurements on spirit leveling lines [50].Recent studies have shown that subsidence, although reduced, is still ongoing [44] [51]- [54].Most of the Po Delta area now lies below sea level, sometimes by as much as 6 -8 meters, and subject to flooding of river branches and other water courses.Combined subsidence, due to natural and anthropogenic causes [55] [56] may lead to serious environmental problems with consequent economic and social implications, especially as regards the mean sea level increase caused by global climate variations [57].The irregular lowering of the terrain also definitely modifies the drainage of the secondary hydraulic network, and increases the rise in the salt wedge for many kilometers: for these reasons, constant monitoring of plano-altimetric changes throughout the Po Delta is essential in implementing territorial defense systems against flooding [14]

Available Data and GPS Measurements
Vertical displacements in the Po Delta have been analysed with geodetic data acquired in the Po Delta in the last few decades.In particular, two conventional leveling lines were measured by the IGMI at different times: line 19 (with 10 benchmarks) at the northern boundary, and line 174_D2 (22 benchmarks) which follows the main stream of the "Po di Pila" in a west-east direction, in the center of the Delta.Line 19 was measured in 1950, providing long-term data over the same benchmarks, although some points have been lost in the meantime.Line 19 also marginally involved the study area.Unfortunately, spirit leveling was used to measure line 174_D2 only in 1978, 1984 and 2005.As regards the accuracy of data, the leveling measurements of the above periods were double run in the flat Delta plain, to exclude the possibility of serious systematic errors affecting leveling [59] [60].They indicate an accuracy of a few mm per km, corresponding to an accuracy of a few cm for elevations of benchmarks separated by distances of up to 20 km, as in this study (Figure 2).This estimate may be considered rather pessimistic of their true uncertainty level, because measurements were double run and the possibility of errors is therefore low.
Comparison of the orthometric elevations of the same benchmarks measured at different time yielded information on vertical movements (Figure 2).Throughout the dataset, point elevations decreased, mainly in the period 1950-1970, when the slope of the straight lines was greater.This effect was simultaneous with or slightly delayed with respect to the most intense phase of methane extraction [37] [39] [47] [52].In the last period, 1970-2005, after extraction ceased, the phenomenon was still not exhausted [61].The differences measured in the various periods are also certainly greater than the accuracy of this survey.
Other archival orthometric elevations of benchmarks from the Veneto Region regard data used to produce the regional CTR (Carta Tecnica Regionale), scales 1:5000 and 1:10,000, from the 1983 photogrammetric flight over the region.The maximum distance between benchmarks is about 40 km, providing an estimate of accuracy double that of previous measurements (a few cm).The 1983 regional benchmarks were well distributed over the area, but site descriptions were missing or not available due to lack of data.The GPS campaign of February 2008 measured the position of regional and/or IGMI benchmarks: double-frequency Leica GPS Systems 1200 were used and the survey was performed in static mode at a sampling rate of 15 s and acquisition time of 60 min (elevation mask 15˚).Benchmarks well distributed in the analysed area were chosen, and their points were selected to permit nearly ideal conditions (e.g., unobstructed horizon view, avoidance of multipath effect).As some of the 1983 benchmarks could no longer be directly measured, points very close to the old ones were recorded (usually no more than a few meters away).Elevation differences between the old points and the new GPS ones were then measured twice with a spirit level and a staff; the latter is accurate to within 2 mm.Some points on the Veneto GPS network planned and measured in 2006 were used as reference stations (masters), and baselines with distances of less than 6 -7 km were measured.At the end of the campaign, GPS data of 43 leveling benchmarks of the 1983 orthometric elevations were obtained.The data were processed by Leica Geo Office software, yielding the coordinates of the UTM-WGS84 reference system.The ellipsoid heights of the 1983 benchmarks were subsequently obtained from the elevation differences between GPS points and benchmarks as deduced from spirit leveling.This last correction did not introduce any significant errors, as the difference in the geoid height of points up to a few meters and could be ignored for the purposes of this study.The technique used to compute GPS ellipsoid heights has a nominal accuracy (standard error) of 2 GPS σ described by the following formula [62]: 0.6 cm 0.5 ppm where S is the distance in kilometers between master and rover receivers, up to 7 km in this case, giving a nominal value for 1 cm

GPS σ ≅
. However, as real errors may be several times higher [63], a more realistic value was adopted, i.e., cm 5

GPS σ ≅
. Since leveling and GPS data are uncorrelated, the error in the elevation of each individual benchmark, following the law of error propagation, is given by the following equation:

The Methodology
When geoid heights are known, the values of vertical displacements can be obtained with the following equation (Figure 3):  N : difference between WGS84 ellipsoid and geoid (geoid height).In this case, the conventional equation of GPS leveling [4] [62] must be corrected to take into account subsidence occurring in the period between the two surveys.N values can be obtained from data on double elevations (surveyed at the same time) of many points uniformly located in the study area and interpolating the data or from gravimetric measurements.Subsequently With this reasoning, we assume that ground subsidence affecting only some near-surface strata has no effect on geoid height N.This is a reasonable simplification: even in the Fennoscandinavian uplift, which expresses a large-scale effect in both space (mantle doming over a region of 1000 -2000 km) and in time (thousands of years), the rate of change in geoid height is one order of magnitude smaller than the rate of change in topography [64].
However, the equation is only correct in theory: in fact, over a period of 25 years, many human interactions have occurred in the territory, sometimes altering some benchmarks and replacing them with other points very close to the old ones, but at different orthometric elevations.In these cases, the 2008 ellipsoid elevations of points arising from unauthorised re-installations replacing the real 1983 benchmarks, now destroyed but previously often at different heights, and then at different orthometric elevations, were re-measured.
As described by Psimoulis et al. [62], this is the major source of errors in benchmarks, because non-homogeneous data are obtained.A procedure to evaluate the congruence of the 1983 orthometric elevations thus had to be applied.Due to the small area and no gravimetric anomalies are documented [65], the N values must be homogeneous.
Some equations were tested to find the best surface allowing the ) parameter to be modeled.Due to territorial morphology and the absence of gravimetric anomalies, the shape of the geoid surface (geoid height) is close to an inclined plane in the analysed area, so that initially simple equations were tested and the data were later improved with increasingly complex forms (until the second order of 1983-2008

N S −
). Testing was carried out near Chioggia, a few kilometers north-west of the Po Delta, for which coeval orthometric and ellipsoid elevations were available (spirit leveling and GPS data).The area is also characterized by the same morphology and geoid surface forms as the Delta.Twenty benchmarks with double elevations were used as control points over an area of about 80 km 2 , and 30 others were used as check-points to evaluate models precision.The benchmarks were located both inside (15) and outside (15) the interpolation area up to a distance of 10 km, so that the extrapolation limits of the models could also be evaluated (Figure 4).These data had the same accuracy as those discussed above.
Twenty double-elevation control points were used to calculate parameters , , , , 1), ( 2), (3), ( 4) and ( 5) with the least-squares method.The N values, obtained with the models used for the check-points, were then compared with the inner and outer measured check-points (Table 1).
Table 1 shows the most accurate data for the model obtained with Equation ( 5): on benchmarks about 10 km from the study area, the maximum difference between N calculated and N measured was 5.5 cm.The data provided by quadratic surface (5), inside the study area, were also compared with Italian gravimetric geoid ITALGEO 2005, characterized by declared precision of about 4 cm, and provided a mean of −0.6 cm and standard deviation of 1.6 cm.These data show that Equation ( 5

N S
− values: for interpolation, the best model obtained from the tests was used, corresponding to a second-order surface with nine parameters (5).Points with anomalous values (those exceeding a threshold of 3 cm in least-squares processing) can be removed.
Unfortunately, points with substantial relative differences in 1983-2008 S values are also removed during statistical processing.However, the gradients could later be recovered in the 2008 1983 h H − differences (with 1983 H correct) by GPS measurements.In this case, Equation (5) becomes: where: ) value.The analysis required excluding 11 out of 43 points (as there were benchmarks with anomalous geoid heights and/or high differential vertical displacements in the period 1983-2008).However, the other 32 points were uniformly distributed over the Delta.The 11 previously excluded points were then assigned the orthometric elevations of 1983, according to the surfaces calculated on the check-points.At the end of the procedure, 43 benchmarks with 1983 and 2008 corrected orthometric elevations were available.

Discussion
The vertical displacements obtained with the procedure described above were integrated with data from spirit leveling measurements (32 benchmarks), interpolated and extrapolated linearly to yield a dataset for the same 1983-2008 period.The 75 differences in the points of the study area (Figure 5(a)) were interpolated according to the IDW deterministic interpolation method, which is able to intercept subsidence peaks which Equation ( 5) cannot record adequately.In this way, a 10-m grid size was extracted by processing about 5 × 10 6 points; only the area enclosed by benchmarks was analysed, as the characteristics of the interpolation method provides unreliable data outside the study area (Table 1) (Figure 5(b)).
The map shows most vertical displacements (30 -35   bordering the Emilia Romagna Region.These results were confirmed by comparisons with two leveling benchmarks at Goro and the forest of Mesola (Figure 1), measured in 1984, 1987, 1999 and 2005 with spirit leveling by the Emilia Romagna regional authorities [73].Direct comparisons for the same period between measured subsidence on the two benchmarks and data extrapolated from the map of Tosi et al. [2], merging various displacement measurements obtained from spirit leveling, DGPS, CGPS and SAR-based interferometry [74], also produced a map of land displacement rates in the period 1992-2002 which, in the Po Delta, provides vertical movements of about 10 mm/year.Baldi et al. [32], using the CGPS method, recently reported subsidence in the eastern part of the province of Ferrara near the Delta of about 6 mm/year in 2006-2011, although their data only partially cover our study area.Thus, it seems realistic to adopt the same subsidence trend for our study area, with considerable vertical displacements in the 1980s, whereas they fell to almost half in the last period.
Figure 6 plots subsidence rates, together with tidal records from the Venice and Ravenna stations (from 1986 to 2013), as these are the stations nearest to the Po Delta with long-term sea level datasets.They show that the sea level has been rising for the last 27 years: with the isostatic and tectonic subsidence rates obtained in this The problem of sea level rise, both past and future, especially in coastal areas with subsidence which enhance the risk of flooding still further, have been discussed by many authors [16]- [19] [27] [75]- [76].Lambeck et al. [29] provide a vertical displacement scenario for the Delta in the year 2100, considering the contribution of isostasy, tectonics and relative sea level rise calculated in 2010, ranging from 315 mm to 1535 mm for the lowimpact ( [78]; +180 mm sea level rise) and high-impact scenarios ( [79]; +1400 mm sea level rise), respectively.However, these values do not include the contribution of recent soil compaction and fluid withdrawal (gas and water).The latter data, extrapolating the trend of Figure 6 until 2100 for artificial subsidence, give new values of 370 and 1590 mm, respectively.In these conditions, Syvitski et al. [58], who place the Po Delta in the highrisk group of "virtually no aggradation and/or very high accelerated compaction", estimate an increase of 50% in the risk of flooding during the 21 st century if global sea level continues to rise rapidly.The results of this work do seem to confirm this risk value and indicate that, due to the ongoing increase in storm events, reinforcement of all kinds of flood defenses will have to be examined carefully.

Conclusions
Vertical displacements in Po Delta area of the Veneto Region were estimated for the 1983-2008 period according to a procedure based on GPS measurements of archival leveling benchmarks and integration with recent spirit leveling values.Due to the uncertain elevations of "historical" benchmarks, a method of geoid height checking was applied.Eight models were compared in a study area near the Po Delta.Equation ( 5) provides the best results here; in areas with different morphology, documented gravimetric anomalies, etc., the equation of local geoid modeling may be different.Thus, after ellipsoid elevations had been transformed into orthometric ones and inhomogeneous benchmarks had been removed or corrected, vertical displacement values could be obtained by comparing homogeneous elevations surveyed over various periods of time.
The above method is particularly efficient as regards time and costs (GPS measurements are faster) than spirit leveling.Full integration with archival orthometric elevations of benchmarks can also be made, recovering useful historical data.Monitoring of the area can then be performed simply with periodical GPS surveys on the same benchmark, which provide vertical displacements for comparison with homogeneous ellipsoid elevations.
The application of this method to the Po Delta yielded subsidence map for the period 1983-2008, interpolating 75 points.In the future, increasing benchmarks will better describe vertical movements in more detail, intercepting any other peaks of subsidence.The major vertical displacements measured in the southern part of the study area were confirmed by comparing the two leveling benchmarks of the Emilia Romagna Region, bordering the Delta.The integration of more recent data [2] [32] also indicates mean vertical subsidence of about 18 mm/year from 1983 to the early 1990s, 10 mm/year from then until the early years of the new century, and 6 mm/year the early 2000s to the present.These data confirm that artificial subsidence is still active, although the rate has decreased (with no linear trend) in the last decade.For this reason, monitoring is still necessary, especially considering the medium/long-term risk of flooding due to sea level rise.

Figure 1 .
Figure 1.(a) Location of Po Delta area; (b) Effects of subsidence, also linked to changes in coastline, causing submergence of large areas; (c) Regional boundaries and indication of study area, of about 600 km 2 (dashed cyan line).

Figure 2 .
Figure 2. Location of two leveling lines (nos.19 and 174_D2) measured by IGMI: analysis of elevations from surveys in different periods on same benchmarks show decreasing values for whole dataset (mainly in period 1950-1970).Colors: elevation trends for each benchmark of above leveling lines.
as the leveling error at very short distances is minimum (up to 1 -2 mm) relative to the error in height determination according to GPS.This procedure extracted double elevations (orthometric 1983 and ellipsoid 2008 values) for 43 benchmarks with σ of about 5 -6 cm.

Figure 3 .
Figure 3. GPS leveling with subsidence: conventional equation of GPS leveling is integrated, to take into account changes occurring between surveys.H: 1983 orthometric elevation; h: 2008 ellipsoid elevation; N: geoid height; S: subsidence in period 1983-2008.
Y e N e XY e XN e YN e X e Y N e

Figure 4 .
Figure 4. Tests area located NW of Po Delta (Veneto Region); 50 points of coeval double elevations: 20 benchmarks and 30 points (15 inside and 15 outside study area) used as checkpoints to evaluate precision of local geoid modeling and limits of surface extrapolation.
) clearly represents the progress of N and suggests using the surface for local modeling of 1983-2008 N S − values in the Po Delta area.The orthometric elevation of the 1983 benchmarks was analysed by modeling the 1983-2008

Figure 5 .
Figure 5. (a) Location of 32 leveling benchmarks and 43 GPS points measured in study area; (b) map of subsidence, 1983-2008, from 75 points calculated with IDW deterministic method (10 m grid size): only area within benchmarks was analysed.Maximum differences (30 -35 cm) located in south-eastern part of study area.

Figure 5 (
b) show differences of about 2 mm, i.e., less than the precision of the method.The three leveling benchmarks in Emilia Romagna (Goro, Mesola, and the mouth of the Po di Goro-for the last, only vertical displacement value for 1999-2005 is available) close to that part of the Po Delta near the Veneto Region, showed no linear displacements between 1983 and 2008: the benchmarks provided values of about 18, 17 and 10 mm/year, measured in1984-1987, 1987-1999  and 1999-2005, respectively.

Figure 6 .
Figure 6.(a) Location of Ravenna and Venice tidal stations, with long-term data of sea level nearest to Po Delta area; (b) Plot of subsidence rates, with tidal records of Venice and Ravenna stations (1986-2013, data from www.mareografico.it).

Table 1 .
cm in 25 years) in the south-eastern part of the Po Delta, Comparisons between N values calculated with various models and data measured at check-points inside and outside study area.