Received 6 March 2014; revised 2 April 2014; accepted 28 April 2014

ABSTRACT

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 me- thods 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 infor- mation on vertical movements when geoid anomalies are known. However, a data checking pro- cedure must be applied if archival orthometric elevations are used, because long-term measure- ments 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 trans- formation from ellipsoid to orthometric elevations, comparisons of the same points and interpola- tions on the study area provided a subsidence map for the 1983-2008 period.

Keywords:

Subsidence, Po Delta, Spirit Leveling, GPS, Local Geoid Modeling

1. Introduction

Ground deformation in areas subject to vertical displacements can be studied by comparing the orthometric ele- vations 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 differ- rential (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 cov- ering 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 compari- sons 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 measure- ments 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 move- ments 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. Histori- cal 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 de- formations 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.

2. 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 devel- opment in the region almost quadrupled the amount of groundwater extracted for industry, agriculture and do- mestic purposes, compared with previous years [23] [24] . The most recent dataset, collected by the regional au- thorities, 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 subsi- dence, 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², extends seaward for about 25 km (Figure 1) [35] , and is cha- racterized 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 rec- lamation 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 can- not 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 magni- tude higher than that due to natural processes. The GPS technique has more recently been used to integrate con- ventional 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 be- low 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 environ- mental problems with consequent economic and social implications, especially as regards the mean sea level in- crease 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 imple- menting territorial defense systems against flooding [14] -[17] [19] [27] [29] [33] [58] .

3. 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 infor- mation on vertical movements (Figure 2). Throughout the dataset, point elevations decreased, mainly in the pe- riod 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 stat- ic mode at a sampling rate of 15 s and acquisition time of 60 min (elevation mask 15˚). Benchmarks well distri- buted 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 de- scribed by the following formula [62] :

where S is the distance in kilometers between master and rover receivers, up to 7 km in this case, giving a no- minal value for. However, as real errors may be several times higher [63] , a more realistic value was adopted, i.e.,.

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

as the leveling error at very short distances is minimum (up to 1 - 2 mm) relative to the error in height determi- nation according to GPS. This procedure extracted double elevations (orthometric 1983 and ellipsoid 2008 val- ues) for 43 benchmarks with σ of about 5 - 6 cm.

4. The Methodology

When geoid heights are known, the values of vertical displacements can be obtained with the following equation (Figure 3):

where:

: vertical displacement between 1983 and 2008;

: ellipsoid elevation measured with GPS in 2008;

: orthometric elevation measured in 1983;

: 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 subsi- dence 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, GPS surveys can yield the parameter which, if the of archival benchmarks are known, yields values. 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 simpli- fication: even in the Fennoscandinavian uplift, which expresses a large-scale effect in both space (mantle dom- ing 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 pre- viously 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-homo- geneous 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 (i.e.) parameter to be modeled. Due to territorial morphology and the absence of gravimetric anomalies, the shape of the geoid sur- face (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). Testing was carried out near Chioggia, a few kilometers north-west of the Po Delta, for which coeval orthome- tric 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², and 30 others were used as check-points to evaluate models pre- cision. 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.

The following interpolation surfaces were analysed:

(1)

(2)

(3)

(4)

(5)

where:

: parameters;

: planimetric coordinates;

: geoid height.

The Kriging geostatistical method, IDW (Inverse Distance Weight) and RBF (Radial Basis Function) deter- ministic interpolation methods were also applied [66] -[72] .

Twenty double-elevation control points were used to calculate parameters of models (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 pro- vided 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 stan- dard deviation of 1.6 cm. These data show that Equation (5) clearly represents the progress of N and suggests using the surface for local modeling of values in the Po Delta area.

The orthometric elevation of the 1983 benchmarks was analysed by modeling the 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 pro- cessing) can be removed.

Unfortunately, points with substantial relative differences in values are also removed during statis- tical processing. However, the gradients could later be recovered in the differences (with correct) by GPS measurements. In this case, Equation (5) becomes:

where:

: nine parameters;

: planimetric coordinates;

(i.e.) 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 uni- formly distributed over the Delta. The 11 previously excluded points were then assigned the orthometric eleva- tions of 1983, according to the surfaces calculated on the check-points. At the end of the procedure, 43 bench- marks with 1983 and 2008 corrected orthometric elevations were available.

5. 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⁶ points; only the area enclosed by benchmarks was analysed, as the characteristics of the interpolation method provides unre- liable data outside the study area (Table 1) (Figure 5(b)).

The map shows most vertical displacements (30 - 35 cm in 25 years) in the south-eastern part of the Po Delta,

bordering the Emilia Romagna Region. These results were confirmed by comparisons with two leveling bench- marks 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 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 availa- ble) 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 in 1984-1987, 1987-1999 and 1999-2005, respectively.

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

work, these values indicate the increased risk of flooding throughout the area.

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 low- impact ([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 high- risk 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.

6. Conclusions

Vertical displacements in Po Delta area of the Veneto Region were estimated for the 1983-2008 period accord- ing to a procedure based on GPS measurements of archival leveling benchmarks and integration with recent spi- rit leveling values. Due to the uncertain elevations of “historical” benchmarks, a method of geoid height check- ing 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 use- ful 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, inter- cepting 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, border- ing 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, espe- cially considering the medium/long-term risk of flooding due to sea level rise.

References