InSAR Time Series Analysis of Natural and Anthropogenic Coastal Plain Subsidence : The Case of Sibari ( Southern Italy )

We applied the Small Baseline Subset multi-temporal InSAR technique (SBAS) to two SAR datasets acquired from 2003 up to 2013 by Envisat (ESA, European Space Agency) and COSMO-SkyMed (ASI, Italian Space Agency) satellites to investigate spatial and temporal patterns of land subsidence in the Sibari Plain (Southern Italy). Subsidence processes (up to ~20 mm/yr) were investigated comparing geological, hydrogeological, and land use information with interferometric results. We suppose a correlation between subsidence and thickness of the Plio-Quaternary succession suggesting an active role of the isostatic compensation. Furthermore, the active back thrusting in the Corigliano Gulf could trigger a flexural subsidence mechanism even if fault activity and earthquakes do not seem play a role in the present subsidence. In this context, the compaction of Holocene deposits contributes to ground deformation. Despite the rapid urbanization of the area in the last 50 years, we do not consider the intensive groundwater pumping and related water table drop as the main triggering cause of subsidence phenomena, in disagreement with some previous publications. Our interpretation for the deformation fields related to natural and anthropogenic factors would be a comprehensive and exhaustive justification to the complexity of subsidence processes in the Sibari Plain.


Introduction
Many coastal and delta plains worldwide are affected by land subsidence phenomena [1,2], which often involve inhabited areas causing conspicuous economic costs.
Examples come from different geodynamic, geological, climatic and social contexts.In Italy, subsidence due to a combination of natural causes and anthropogenic activities is observed in the Venice coastland and the Po delta plain [3][4][5][6], while in the United States in correspondence of the coastal areas of Southern Louisiana and Mississippi [7][8][9][10][11].Other examples are observed in the Yellow River and Yangtze China Deltas [12,13], along the Taiwan coastland [14][15][16], in the Indonesian Semerang city [17], in Thessaloniki coastal plain and municipality region (Northern Greece) [18,19].
Over the last decades, the land subsidence monitoring has been significantly improved thanks to the earth observation techniques based on Interferometric Synthetic Aperture Radar (InSAR).Furthermore, the recent advances in radar satellite capabilities and techniques based on the interferometric analysis of large datasets [20][21][22][23][24][25][26][27][28][29] have allowed even better spatial and temporal resolutions.
We applied multitemporal differential SAR Interferometry technique and specifically the Small Baselines Subset (SBAS) approach [20] to investigate the ground deformations of the Sibari Plain (SP), located in the northeastern sector of the Calabria Region (Southern Italy) covering an area of ~500 km 2 .This Plain involves an important economic weight due to the agricultural production and a cultural appeal related to the presence of the ancient Sybaris, a powerful Greek colony in Magna Grecia founded in 720 BC.The area is also susceptible to flooding risk, reduced by means of a drainage channels network.
The used SAR datasets cover the temporal interval spanning from May 2003 and September 2013 for both ascending and descending orbits acquired by Envisat and COSMO-SkyMed satellites (for details see Table 1).Envisat is a European Space Agency (ESA) satellite operating in the C-Band (5.6 cm of wavelength) launched in the 2002 and operating up to 2012, with a revisiting period of 35 days, covering an area of about 100 ˆ100 km.The COSMO-SkyMed mission consists of a constellation from the Agenzia Spaziale Italiana (ASI) in the X-Band frequency (3.2 cm of wavelength).The first satellite of the constellation has been launched in June 2007 and the constellation (four satellites) has been fully operative by the end of 2010.The used acquisitions (Himage mode) are characterized by a swath of about 40 ˆ40 km 2 with a revisiting time period of 16 days.Moreover, the availability of ascending and descending datasets allowed us to discriminate the vertical and east-west displacement components.
Finally, the retrieved InSAR results were validated by means of GPS measurements then analyzed and interpreted considering the available geological and hydrogeological information as well as new data collected during our field surveys.

Geological Setting
The SP is located along the boundary between Calabrian Arc and Southern Apennines (Figure 1), with the Pollino massif to the north and Sila massif to the south.The Calabrian Arc represents a fault-bounded continental fragment within the western Mediterranean orogeny.The tectonic evolution is related to the subduction and rollback of Ionian oceanic lithosphere and the slow convergence between the Eurasian and African-Adriatic continental plates [30][31][32].The Southern Apennines are a NW-SE oriented segment of the Apennines thrust belt.It is characterized by a duplex structure, which consists of two thrusts belts overlapping the Apulian platform [33][34][35][36].
In the SP, the chain is prevalently composed by Apenninic terranes overlapped by thin Calabrian Arc units [37].The first ones, outcropping along the northwestern side of the SP, are made by Mesozoic-Tertiary carbonatic succession and "flyschiod" deposits overlapped by Plio-Pleistocene siliciclastic succession [38][39][40].Instead, the Calabrian Arc terranes, cropping out along the southern margin, consist of igneous and metamorphic rocks [41] (Figure 1).
In the central sector of the SP, the basement is overlapped by a thick Miocene-Holocene succession; the mio-pliocenic deposits represent the infilling of the Sibari-Corigliano Basin, which evolution is strictly connected with the Sangineto Line activity [42,43], and they are deformed by fault activity with formation of morphological/structural highs and lows [44].In [45], the authors suggest the existence of WNW-ESE trending shallow-crustal folds, developed within a recent and still active transpressional field.
oriented segment of the Apennines thrust belt.It is characterized by a duplex structure, which consists of two thrusts belts overlapping the Apulian platform [33][34][35][36].
In the SP, the chain is prevalently composed by Apenninic terranes overlapped by thin Calabrian Arc units [37].The first ones, outcropping along the northwestern side of the SP, are made by Mesozoic-Tertiary carbonatic succession and "flyschiod" deposits overlapped by Plio-Pleistocene siliciclastic succession [38][39][40].Instead, the Calabrian Arc terranes, cropping out along the southern margin, consist of igneous and metamorphic rocks [41] (Figure 1).

Figure 1.
Simplified lithological map of the study area with projection of the faults of ITaly HAzard from CApable faults database [46] and of the Sybaris fault zone (f.z.) [47].[46] and of the Sybaris fault zone (f.z.) [47].
The boundaries between the SP and Pollino and Sila massifs are marked by dislocations produced along plio-holocenic high angle faults.The latter is represented in the NW sector by the NW-SE normal fault of the Sangineto Line [48,49], while in correspondence of the Sila Massif it consists of the WNW-ESE Rossano-Corigliano fault (Figure 1).The recent activity of these faults is still debated [49][50][51][52][53].
The existence of NE-SW trending normal fault in the area of the Crati Delta is suggested by [48]; recent coseismic evidences related to a NE-SW structural lineament, maybe active during the period between the II and the VII-IX century A.D., are observed in Parco del Cavallo and Casa Bianca archaeological sites [47].
Five to eleven different order of marine terraces are recognized along the outer limit of the plain and uplift rates greater than 3.5 mm/yr are estimated [50,54,55].
Holocene deposits are the result of the post Last Glacial Maximum (LGM) transgression and the following (started about 6 ka B.P.) normal regression related to the Crati Delta progradation [56,57].The Holocene evolution occurred in a tectonic-controlled setting, which drives the thickness and the facies association lateral variability of the deposits [57].The post LGM succession consists of marine clayey-silty deposits, overlapping a late Pleistocene coarse-grained unit, passing upward to costal, deltaic and continental ones [57][58][59].
In the Northern sector, the main geomorphological elements are the alluvial fans of Raganello River, Satanasso Fiumara and Saraceno Fiumara.

Historical Subsidence
Subsidence is a well-known phenomenon in the Sybaris archaeological area (Figure 2), testified by the presence of three overlapping ancient towns, the Greek Sybaris (720-510 BC), the Hellenistic Thurii (444-203 BC) and the Roman Copiae (193 BC), presently buried between 7 and 3.5 m [60].
In [61], the author calculates subsidence rates ě0.The subsidence in the archaeological area starts in the late Pleistocene up to the Holocene and is due by a combination of causes as neotectonic activity, glacio-eustatic variations and sediments compaction [62][63][64].In the archaeological area the subsidence variability is correlated to the lateral variation of facies [65] and the most of geotechnical subsidence can be ascribed to a very compressible clay layer, laterally discontinuous, between 35 and 40 m of depth [62,65].Since the 1950s, ~20 cm of subsidence have been recorded and ascribed to groundwater exploitation and related primary consolidation [65,66].The subsidence in the Casa Bianca and Parco del Cavallo sites is confirmed by [67], which instead calculates a mean uplift rate of ~0.5 mm/yr in the last 11.2 ka for the Stombi site and hypothesizes that subsidence started 4ka B.P. and is ascribed to the deposition of fine compressible sediments, so excluding the tectonic contribution.The author of [68] highlights a temporal variability for the mean subsidence rates showing a peak of ~5-6 mm/yr in the Early Holocene whilst stability or small uplift in historical age (~3-1.

InSAR Adopted Tecnique
The classical differential SAR Interferometry (InSAR) is a technique that allows us to estimate the ground surface movements occurred between two different passes of the satellite over the same area using Synthetic Aperture Radar (SAR) data.The phase difference at each pixel is calculated after proper image co-registration resulting in a new image, called interferogram, an interference pattern made up of fringes [69].Each fringe represents a ground movement along the sensor line of sight (LOS) equal to λ/2 (where λ is the adopted radar wavelength).
There are many physical phenomena contributing to phase measurements: phase variations within a pixel, the contribution of orbital variations, topography, atmosphere, and displacement.These

InSAR Adopted Tecnique
The classical differential SAR Interferometry (InSAR) is a technique that allows us to estimate the ground surface movements occurred between two different passes of the satellite over the same area using Synthetic Aperture Radar (SAR) data.The phase difference at each pixel is calculated after proper image co-registration resulting in a new image, called interferogram, an interference pattern made up of fringes [69].Each fringe represents a ground movement along the sensor line of sight (LOS) equal to λ/2 (where λ is the adopted radar wavelength).
There are many physical phenomena contributing to phase measurements: phase variations within a pixel, the contribution of orbital variations, topography, atmosphere, and displacement.These contributions are estimated and then removed so to obtain only the effective ground displacement.
Our goal is to investigate the temporal evolution of the detected deformations.To this aim several different approaches have been proposed in the last decade to generate time series of ground displacement, capable of measuring the ground velocities rates with accuracies of ~1 mm/yr [70].In this work, we applied the Small Baseline Subset InSAR technique (SBAS; [20]).The used SAR images were coupled based on constraints of small temporal and spatial baselines.The SBAS algorithm combines these acquisitions, also included in different interferometric subsets, using a minimum norm criteria combination of the velocity deformation, based on the Singular Value Decomposition (SVD) method.During the SBAS processing the estimation and removal of temporally decorrelated atmospheric artefacts are performed using double-pass filtering in the temporal and spatial domains, as explained in [20].Moreover, we used the SRTM DEM (90-m resolution) for the subtraction of topography phase (http://www2.jpl.nasa.gov/srtm).Once the displacement time series is retrieved, the mean ground velocity in the time period covered by the data was calculated.
In detail, we applied the SBAS technique to four image datasets, two acquired from the European Space Agency (ESA) Envisat sensor and two from the Italian Space Agency (ASI) COSMO-SkyMed satellites.The first group contains 32 Envisat images acquired on the descending orbit (track 222, frame 2803), in the period 2004-2010 imposing 700 days as maximum temporal baseline and 350 meters for the maximum perpendicular one and considering 100 pairs.Instead, the ascending dataset (track 86, frame 789) is formed by 38 ASAR images acquired in the period 2003-2010 setting 900 days for the maximum temporal baseline and 350 meters for the maximum spatial one resulting in 135 interferograms.The second group (COSMO-SkyMed data) was composed of 28 images from the ascending orbit considering 180 days for the maximum temporal baseline and 700 meters as maximum spatial baseline resulting in 137 interferograms, spanning the interval 2011-2013 and 12 images from the descending orbit, acquired in the interval 2011-2012 setting 100 days for the temporal and 1000 meters for the spatial maximum baselines respectively forming 47 interferograms.We used the Sarscape software (Sarmap, CH) applying a looking factors equal to 20 and 4 for the azimuth and range direction, respectively for the ascending and descending Envisat processing obtaining a ground resolution of 80 meters and 11, 12 for the COSMO-SkyMed case with a final ground posting of 25 meters; later downgraded to 100 meters for the post-processing analysis.The initial performed multi-looking operation reduced the image radar noise (speckle) and increased the signal to noise ratio.
During the processing the precise orbital files (Envisat case) were used to estimate and remove orbital errors and Ground Control Points (GCPs) were selected, especially along the frame borders (in slant range geometry) and possibly over stable areas, to remove residual ramps.A stable point (under the geological point of view) was chosen as reference point (red star) for the retrieved mean ground velocity maps (Figure 3).
Initially, we validated our results comparing the ascending and descending velocity map after shifting the former with respect to the reference point of the latter preventing the subsidence areas.We also validated the Envisat results comparing with the available GPS in the whole SAR frames (Figure 3).Firstly, we projected the GPS velocities onto the ascending and descending SAR LOS respectively then we averaged the InSAR velocities in correspondence of the GPS benchmarks considering a buffer area of 200 meters.The comparison between GPS and InSAR velocities is shown in Figure 3a,b for the ascending and descending tracks.Except for very few GPS benchmarks, the velocities differences are within 1.5 mm/yr, which is the associated error to the multi-temporal SAR outputs.Unfortunately, it has not been possible to perform the same validation method for the CSK results because not enough CGPS are available.However, considering the common similar patterns and trends for both Envisat and CSK ascending and descending retrieved maps, we convinced the good quality of the obtained CSK results.Note that for the CSK processing we focused our study on the Sibari Plain.
Remote Sens. 2015, 7 where u2d, e2d, n2d, u2a, e2a and n2a are the descending and ascending cosine directors respectively, and Asc and Dsc are the ascending and descending mean velocity maps.Concerning the third condition, we imposed the North equal to the interpolated CGPS one.The final Vertical and East component (Figure 4) result quite sparse, this mainly depends on the difficulty to obtain a large coverage for the ascending and descending velocity maps due to the land cover and agricultural use of the investigated area leading to a fast coherence loss.Moreover, about the East component some Concerning the third condition, we imposed the North equal to the interpolated CGPS one.The final Vertical and East component (Figure 4) result quite sparse, this mainly depends on the difficulty to obtain a large coverage for the ascending and descending velocity maps due to the land cover and agricultural use of the investigated area leading to a fast coherence loss.Moreover, about the East component some overestimated values are present in the North Eastern part probably due to unwrapping errors due to the presence of the mountains.
We performed the previous post-processing steps for the full SAR frames and only at a later stage; we focused our analysis and interpretations referring to the East and Vertical component about the Sibari Plain area.

Remote Sens. 2015, 7 8
We performed the previous post-processing steps for the full SAR frames and only at a later stage; we focused our analysis and interpretations referring to the East and Vertical component about the Sibari Plain area.

InSAR Results
In the following, we report the deformation pattern of SP at a large scale retrieved by InSAR processing and discuss the results in detail using available tectonic, stratigraphic, geomorphological and hydrogeological information.We investigate natural and/or human causes of the vertical and horizontal surface displacements.
The distribution of Vertical component shows a coastal area (from Villapiana Lido to Marina di Schiavonea) characterized by subsidence with velocity up to ~´20 mm/yr (Envisat and COSMO-SkyMed sensors); on the contrary, boundary areas of SP show ~0 mm/yr on the vertical velocity component (Figure 4).
The horizontal Eat-West component shows low displacement values that are inside the errors range and are not significant.
In the archaeological area, the subsidence rate is ~2 mm/yr (Envisat) and ~3 mm/yr (COSMO-SkyMed) for the Parco del Cavallo site and ~2 mm/yr (Envisat) for the Casa Bianca site.These values are comparable to the velocities calculated by previous authors [61,67].

Discussion
We analyzed different causes contributing to subsidence phenomena.To this aim, we compared the InSAR results in time and space with the geodynamic and structural settings, the seismicity distribution, the spatial variation of the Plio-Quaternary and Holocene deposits successions thickness and compaction, depth variation of the water table and urban sprawl.

Geodynamic Setting
The SP is located at the junction of the Calabrian Arc and the Southern Apennine.Moving few kilometers towards SE from the SP, the geodynamic setting is complicated by the presence of the junction (Apulia Escarpment) (Figure 1) between the thin Ionian basin crust (southward) and the thick crust of the Apulian platform (northward).In the Ionian Sea, in front of the SP an active oblique-contractional belt (the Amendolara Ridge) is recognized.The latter is due to the combined effects of subduction retreat of the Ionian slab underneath Calabria and stalling of Adriatic slab retreat underneath the Apennines [71].The belt consists of the alignment of three anticlines (Amendolara, Rossano and Cariati) bounded by a main SW-verging back-thrust (Figure 5a).
A flexural subsidence, probably triggered by a back-thrusting [72], is considered one of the mechanisms to explain the ground subsidence detected by InSAR results along coastal sector of the SP (Figure 5b).

Structural Setting and Earthquakes
Possible relationship between SP subsidence and recent tectonic activity relies on the analysis of structural lineaments reported into the Italy Hazard from Capable faults database [46].According to the main faults bounding the SP, we subdivided it in three sectors (Figure 5a) and with foreseen different deformation patterns for the three sectors from InSAR ascending and descending time series (Figure 6).The analysis of the averaged Envisat displacement time series show very similar trends for the three investigated areas, although with different absolute values.On the other hand, COSMO-SkyMed time series show a different deformation pattern for each sector.The bathymetric map (modified from [73]) shows the presence of three morphological highs (Amendolara, Rossano and Cariati highs).The main structural elements (thrusts and anticlines) are extracted from [71].In the map is showed the hypocenters location of the earthquakes occurred between 2003 and 2010 (ISIDe [74]) with the relative magnitude and hypocentral depth.Onshore capable faults from [46] and focal mechanisms from [45].Yellow dashed lines are the boundaries of three selected sectors (1.Corigliano; 2. Crati Delta; 3. Sibari-Villapiana) for Envisat and COSMO-SkyMed time series analysis (Figure 6).(b) Sketch (based on the interpretation of the seismic profile F75-89 of [45]) representing the trigger of the flexural subsidence due to the overlap of the Amendolara Ridge back thrusts on Calabrid and Southern Apennines units.
In detail, the ascending time series reveal a general trend characterized by an increase of the negative displacements but with greater absolute values in the southern sector.The descending time series highlight an increase of negative displacements in the Corigliano sector with values between 10 mm and −5 mm in the Villapiana-Sibari sector, and an expected increasing positive trend in the Crati Delta sector.This trend can be related to the presence of a horst bounded by Crati and Timparelle faults [48].[73]) shows the presence of three morphological highs (Amendolara, Rossano and Cariati highs).The main structural elements (thrusts and anticlines) are extracted from [71].In the map is showed the hypocenters location of the earthquakes occurred between 2003 and 2010 (ISIDe [74]) with the relative magnitude and hypocentral depth.Onshore capable faults from [46] and focal mechanisms from [45].Yellow dashed lines are the boundaries of three selected sectors (1.Corigliano; 2. Crati Delta; 3. Sibari-Villapiana) for Envisat and COSMO-SkyMed time series analysis (Figure 6).(b) Sketch (based on the interpretation of the seismic profile F75-89 of [45]) representing the trigger of the flexural subsidence due to the overlap of the Amendolara Ridge back thrusts on Calabrid and Southern Apennines units.
In detail, the ascending time series reveal a general trend characterized by an increase of the negative displacements but with greater absolute values in the southern sector.The descending time series highlight an increase of negative displacements in the Corigliano sector with values between 10 mm and ´5 mm in the Villapiana-Sibari sector, and an expected increasing positive trend in the Crati Delta sector.This trend can be related to the presence of a horst bounded by Crati and Timparelle faults [48].
A comparison between ground deformation pattern and seismicity was performed taking into account the earthquakes occurrences and magnitude [74] during the time covered by SAR data (Figure 5a).No significant earthquakes struck the study area excepting for one shallow offshore M 4.5 event on 27 June 2006 [74] with none ground deformation effects (Figure 6).
Therefore, we can exclude the presence of seismically-induced displacement and we hypothesize that the structural settings could indirectly influence the ground subsidence phenomena.
the earthquakes occurrences and magnitude [74] during the time covered by SAR data (Figure 5a).No significant earthquakes struck the study area excepting for one shallow offshore M 4.5 event on 27 June 2006 [74] with none ground deformation effects (Figure 6).
Therefore, we can exclude the presence of seismically-induced displacement and we hypothesize that the structural settings could indirectly influence the ground subsidence phenomena.Figure 6.Envisat and COSMO-SkyMed both ascending and descending time series (error = ±1 mm) for the three selected sectors delimited by ITHACA [46] faults (location in Figure 5).Black dotted line represents the earthquakes (M = 4.5) occurred on 27 June 2006.

Role of the Plio-Quaternary Succession Load
Possible relationships between subsidence and spatial variation of the Plio-Quaternary succession thickness in the Crati Delta area (Figure 7a) as reconstructed from multichannel seismic profiles and well data [75] was investigated.The thickness observed in the seismic profile has been converted from Two Way Travel time to meter using an average velocity of 2000 m/s for the Pliocene and Pleistocene deposits based on the sonic logs, available for Thurio and Ogliastrello wells, analyzed by [43].
The Envisat Vertical component superimposed to the Plio-Quaternary succession isopach map (Figure 7b) shows a prevalent distribution of vertical displacement rates toward the coastline, where the deposits are thicker.
We compared ground velocity profiles, for the Envisat Vertical component, with the corresponding profiles of the Plio-Quaternary deposits thickness.The comparison (Figure 7c) shows a good correlation between the general subsidence trend and the Plio-Quaternary succession thickness.We observed higher subsidence values in the E and NE sectors of Crati Delta, where the Plio-Quaternary succession is

Role of the Plio-Quaternary Succession Load
Possible relationships between subsidence and spatial variation of the Plio-Quaternary succession thickness in the Crati Delta area (Figure 7a) as reconstructed from multichannel seismic profiles and well data [75] was investigated.The thickness observed in the seismic profile has been converted from Two Way Travel time to meter using an average velocity of 2000 m/s for the Pliocene and Pleistocene deposits based on the sonic logs, available for Thurio and Ogliastrello wells, analyzed by [44].
The Envisat Vertical component superimposed to the Plio-Quaternary succession isopach map (Figure 7b) shows a prevalent distribution of vertical displacement rates toward the coastline, where the deposits are thicker.
We compared ground velocity profiles, for the Envisat Vertical component, with the corresponding profiles of the Plio-Quaternary deposits thickness.The comparison (Figure 7c) shows a good correlation between the general subsidence trend and the Plio-Quaternary succession thickness.We observed higher subsidence values in the E and NE sectors of Crati Delta, where the Plio-Quaternary succession is thicker.In the SW sector, the reduced thickness corresponds with lower subsidence values.Lower values are also detected in the S-SE sector (Corigliano harbor area), where the presence of shallow igneous-metamorphic bedrock drives the Plio-Quaternary succession thickness.
Furthermore, we considered the ascending and descending ground velocity values in correspondence of the seismic profiles and wells locations, correlating the ground velocity with the thickness of Plio-Quaternary deposits (Figure 8).Both ascending and descending velocities increase with the growth of the succession thickness.
Based on the above analysis results, we found a direct correlation between subsidence spatial trends and distribution of the thickness of Plio-Quaternary deposits, which amplify the subsidence phenomenon as well documented during the evolution of the sedimentary basins [76][77][78][79][80][81].

Holocene Deposits
A possible correlation between subsidence and Holocene sediments compaction was investigated focusing on the lithological features, thickness and stratigraphic architecture of the Holocene deposits.
We analyzed 200 boreholes (Demanio Idrico Prov.Cosenza, [82,83]) and literature data [57,61,65] to define lithology, vertical-lateral variations and thickness of the post-LGM deposits.An isopach map of the Holocene succession superimposed with the Envisat Vertical component was generated (Figure 9a).The greater vertical displacement values are located in the area of the Holocene deposits accumulation but the larger subsidence rates do not correspond to the major depozones.
To well investigate the second topic, we compared each Vertical displacement value of with the corresponding thickness of the Holocene deposits (Figure 9b).We also performed a comparison between thickness and vertical component in correspondence of each borehole, if SAR data are available (Figure 9c).
The analysis reveals that the Holocene sediments compaction does not play a dominant role in the subsidence of the SP area, on the contrary as described for other delta plains [5,7].
Subsidence rates of 3-5 mm/yr induced by to the Holocene deposits consolidation due to phreatic and confined water tables drop are calculated by [66] and are not sufficient to justify the vertical displacements inferred by SAR data.
Although the minor role of Holocene deposits compaction, the boreholes analysis highlights a clear influence of the lateral variations of the sediments lithology for the subsidence phenomena.In detail, we took into account the ground deformation patterns from ascending and descending Envisat data, and the boreholes into the Saraceno and Satanasso alluvial fans areas (Figure 10a), considering the thickness of fine-grained (silt and clay) deposits (Figure 10b,c).Greater ground velocities are usually retrieved close to the boreholes recording a major thickness of fine-grained sediments while the lower values are in correspondence of the coarse-grained (gravel and sand) deposits which made up the alluvial fan bodies.Thus, the lateral variation of the fine-grained materials thickness creates differential land subsidence.

Land Use and Historical Evolution
SP represents a geomorphological system with a rapid evolution in historical and recent times controlled by geological processes and anthropogenic activities.The present morphological setting results from the works of a 1960s-1990s land reclamation project aimed to convert a marshy area to an agricultural zone.Several widespread agricultural areas (e.g., orchards, arable soils) with localized urban settlements [84] are present in SP.
In order to define land use variations, we analyzed the 1954 and 1998 aerial photos; for a detailed analysis of urban area growth effects.We considered three areas (Marina di Sibari; Laghi di Sibari and C.da Ricota; see Figure 4), characterized by the same subsoil stratigraphy, and compared the Envisat ascending and descending displacement time series concerning urban and agricultural zones respectively (Figure 11).
In Laghi di Sibari and C.da Ricota sites, we observe a progressive temporal differentiation in the subsidence values that increase faster in urban areas reaching a maximum difference compared with agricultural ones of ~35 mm in August 2010 (for the descending track) and of ~25 mm in September 2010 (for the ascending one).The subsidence values are very similar for the two considered sectors in the Marina di Sibari zone, where the descending time series show a little increase (~5 mm) of ground deformation values into the urban area.
We observed a diffuse subsidence irrespective of the land use, with some increments of deformation velocity in correspondence of Laghi di Sibari and C.da Ricota settlements.
Therefore, we suggest that the urban-induced loading represents an incremental factor of subsidence.
Marina di Sibari zone, where the descending time series show a little increase (~5 mm) of ground deformation values into the urban area.
We observed a diffuse subsidence irrespective of the land use, with some increments of deformation velocity in correspondence of Laghi di Sibari and C.da Ricota settlements.
Therefore, we suggest that the urban-induced loading represents an incremental factor of subsidence.

Water Table Variations
SP is characterized by two well-defined aquifers: the shallower one (from soil surface to −20/−30 m above sea level) separated by clayey and silty-clayey layer and the deeper one (from −50/−60 m a.s.l.) [85].Both the aquifers are characterized by intense water exploitation.The piezometric level variations of the shallow one shows 5 meters of drop between 1930s and 2002 [86,87].
We investigated to find possible correlation between groundwater exploitation and ground deformation comparing the spatial distribution of phreatic water table variations in the time interval

Water Table Variations
SP is characterized by two well-defined aquifers: the shallower one (from soil surface to ´20/´30 m above sea level) separated by clayey and silty-clayey layer and the deeper one (from ´50/´60 m a.s.l.) [85].Both the aquifers are characterized by intense water exploitation.The piezometric level variations of the shallow one shows 5 meters of drop between 1930s and 2002 [86,87].
We investigated to find possible correlation between groundwater exploitation and ground deformation comparing the spatial distribution of phreatic water table variations in the time interval 1930s-2002 [86] and the Vertical velocity distribution from Envisat and COSMO-SkyMed data in the following 2003-2013 time interval (Figure 12).In the most recent period, we have observed that subsidence is present in areas showing rise, drop and stability of the water table.
Considering the effect of water table variation between 2002 (data from [86]) and 2013 (this work), we analyzed the distribution of the Vertical velocity (from Envisat and COSMO-SkyMed data, Figure 12).The subsidence location does not show a strict correlation with the water table drop.

Remote Sens. 2015, 7 16
1930s-2002 [86] and the Vertical velocity distribution from Envisat and COSMO-SkyMed data in the following 2003-2013 time interval (Figure 12).In the most recent period, we have observed that subsidence is present in areas showing rise, drop and stability of the water table.
Considering the effect of water table variation between 2002 (data from [86]) and 2013 (this work), we analyzed the distribution of the Vertical velocity (from Envisat and COSMO-SkyMed data, Figure 12).The subsidence location does not show a strict correlation with the water table drop.We selected 6 wells (Figure 1), with piezometric records in 2002 and 2013, and we analyzed the Envisat time series for both ascending and descending orbit (Figure 13) considering a circular buffer of 500 m radius.For the two wells No. 2 and 3, showing a groundwater drawdown, we observe a displacement trend similar to the No. 1, 4, 5 and 6 ones characterized by water table rise.Moreover, the time series show smaller displacements rates in correspondence of the well No. 2. There is no evidence of correlation between water table variation and displacement time series, probably because the progressive groundwater level stabilization from 2002 to 2013.In contrast, the groundwater level drop We selected 6 wells (Figure 1), with piezometric records in 2002 and 2013, and we analyzed the Envisat time series for both ascending and descending orbit (Figure 13) considering a circular buffer of 500 m radius.For the two wells No. 2 and 3, showing a groundwater drawdown, we observe a displacement trend similar to the No. 1, 4, 5 and 6 ones characterized by water table rise.Moreover, the time series show smaller displacements rates in correspondence of the well No. 2. There is no evidence of correlation between water table variation and displacement time series, probably because the progressive groundwater level stabilization from 2002 to 2013.In contrast, the groundwater level drop between 1930s and 2002 [86,87] could have been caused by a more recent consolidation process of soft materials as observed by [88] in the Alto Guadalentin Basin.
We selected 6 wells (Figure 1), with piezometric records in 2002 and 2013, and we analyzed the Envisat time series for both ascending and descending orbit (Figure 13) considering a circular buffer of 500 m radius.For the two wells No. 2 and 3, showing a groundwater drawdown, we observe a displacement trend similar to the No. 1, 4, 5 and 6 ones characterized by water table rise.Moreover, the time series show smaller displacements rates in correspondence of the well No. 2. There is no evidence of correlation between water table variation and displacement time series, probably because the progressive groundwater level stabilization from 2002 to 2013.In contrast, the groundwater level drop between 1930s and 2002 [86,87] could have been caused by a more recent consolidation process of soft materials as observed by [88] in the Alto Guadalentin Basin.For the confined aquifer, the lack of historical data did not allow us to analyze the effects of its variations on subsidence.
Therefore, we conclude that water table variations are only an incremental factor of subsidence in some areas as the Sybaris archaeological area (see Figure 2 to localize Sybaris site), where a well-points system (80 l/s) works every time to maintain the excavations above the water table.Moreover, in these areas, a water depletion in compressible soils could trigger an amplified and long-term subsidence as shown in other Quaternary basins [88,89].

Conclusions
InSAR results from data acquired by Envisat and COSMO-SkyMed satellites, and spanning from 2003 to 2013, are able to detect ground displacements in the SP area (Calabria, Southern Italy).The ground deformation field is dominated by a widespread-subsidence with rate up to ~20 mm/yr.
We can summarize that the location of the SP in an area marked by geodynamic and geological complexity precluding the attribution of the ground deformation triggering factor to a unique subsidence mechanism.
We suppose that the active oblique-contractional belt in the Ionian Sea in front of SP [45] can trigger a flexural subsidence mechanism due to the back thrusting load above the underlying Calabrid and Southern Apennines Units.Furthermore, the observed correlation between Plio-Quaternary sediments thickness and ground deformation suggests that the isostatic compensation plays a major role in the SP present subsidence.We highlight the absence of correlation between the occurring earthquakes close to the main capable faults, so excluding the contribution of faulting activity and earthquakes for the present SP subsidence phenomena, according to [67] (based on geomorphological data); we hypothesize that the structural settings could indirectly influence the ground subsidence phenomena.
Another minor triggering mechanism of subsidence is represented by the compaction of the compressible fine-grained sediments related to the Holocene post-LGM transgression and following Crati Delta progradation.We show the effect of the lateral transition between fine-and coarse-grained Holocene sediments on the ground deformation spatial pattern.
Anthropogenic activities can be considered as an incremental factor of the SP subsidence.In fact, we observe higher subsidence velocities in the more anthropized areas, where soil compaction is amplified by urban-induced loading.In some area showing excessive water withdrawals (e.g., Sybaris archaeological site), water table variations may increase the deformation velocity because can be the origin of delayed consolidation process in soft materials.
The presence of subsidence phenomenon in the whole SP and the complicated identification of a main triggering mechanism, due to the peculiar geological setting, is also in agreement with [90].
The results of this study show as the deformation field of complex basins may be due to coexisting regional and local natural processes and anthropic activities.To better assess and monitor the hazard(s) related to the detected ground deformation phenomena, a multiparametric monitoring of SP area would be recommended.

Figure 1 .
Figure 1.Simplified lithological map of the study area with projection of the faults of ITaly HAzard from CApable faults database[46] and of the Sybaris fault zone (f.z.)[47].
3 ka B.P.) for the Casa Bianca and Parco del Cavallo sites, no subsidence in the Stombi site and uplift (~1.5 mm/yr) 6 km toward SE from archaeological area; this differential subsidence is mainly controlled by local tectonic structures.Remote Sens. 2015, 7 5 for the mean subsidence rates showing a peak of ~5-6 mm/yr in the Early Holocene whilst stability or small uplift in historical age (~3-1.3ka B.P.) for the Casa Bianca and Parco del Cavallo sites, no subsidence in the Stombi site and uplift (~1.5 mm/yr) 6 km toward SE from archaeological area; this differential subsidence is mainly controlled by local tectonic structures.

Figure 2 .
Figure 2. The archaeological sites of the ancient Sybaris (see inset in Figure 1 for its location).

Figure 2 .
Figure 2. The archaeological sites of the ancient Sybaris (see inset in Figure 1 for its location).

Figure 3 .
Figure 3. (a) Ascending and (b) descending line of sight (LOS) ground velocity maps.Red star represents the stable point (stable means that the mean velocity is 0 mm/yr) used as reference point.In the maps are also plotted the GPS benchmarks used during the validation of the Envisat results.The diagrams show the comparison between the LOS Synthetic Aperture Radar (SAR) velocities and GPS.The correlation coefficient is respectively 0.51 and 0.65 for the ascending and descending cases, while the mean squared error is 1.1 mm and 1.2 mm.

Figure 3 .
Figure 3. (a) Ascending and (b) descending line of sight (LOS) ground velocity maps.Red star represents the stable point (stable means that the mean velocity is 0 mm/yr) used as reference point.In the maps are also plotted the GPS benchmarks used during the validation of the Envisat results.The diagrams show the comparison between the LOS Synthetic Aperture Radar (SAR) velocities and GPS.The correlation coefficient is respectively 0.51 and 0.65 for the ascending and descending cases, while the mean squared error is 1.1 mm and 1.2 mm.

Figure 4 .
Figure 4. Mean deformation velocity components computed from (a,b) Envisat and (b,c) COSMO-SkyMed dataset.(a,c) Vertical velocity component; positive values indicate uplift and negative values subsidence.(b,d) East component; positive values indicate eastward displacement and negative values indicate westward displacement.The area covered by COSMO-SkyMed data is represented by inset (1) in the panel a.

Figure 4 .
Figure 4. Mean deformation velocity components computed from (a,b) Envisat and (b,c) COSMO-SkyMed dataset.(a,c) Vertical velocity component; positive values indicate uplift and negative values subsidence.(b,d) East component; positive values indicate eastward displacement and negative values indicate westward displacement.The area covered by COSMO-SkyMed data is represented by inset (1) in the panel a.
Remote Sens. 2015, 7, 16004-16023 Remote Sens. 2015, 7 10 investigated areas, although with different absolute values.On the other hand, COSMO-SkyMed time series show a different deformation pattern for each sector.

Figure 5 .
Figure 5. (a) The Amendolara Ridge located in the Ionian Sea in front of the Sibari Plain.The bathymetric map (modified from[73]) shows the presence of three morphological highs (Amendolara, Rossano and Cariati highs).The main structural elements (thrusts and anticlines) are extracted from[71].In the map is showed the hypocenters location of the earthquakes occurred between 2003 and 2010 (ISIDe[74]) with the relative magnitude and hypocentral depth.Onshore capable faults from[46] and focal mechanisms from[45].Yellow dashed lines are the boundaries of three selected sectors (1.Corigliano; 2. Crati Delta; 3. Sibari-Villapiana) for Envisat and COSMO-SkyMed time series analysis (Figure6).(b) Sketch (based on the interpretation of the seismic profile F75-89 of[45]) representing the trigger of the flexural subsidence due to the overlap of the Amendolara Ridge back thrusts on Calabrid and Southern Apennines units.

Figure 5 .
Figure 5. (a)The Amendolara Ridge located in the Ionian Sea in front of the Sibari Plain.The bathymetric map (modified from[73]) shows the presence of three morphological highs (Amendolara, Rossano and Cariati highs).The main structural elements (thrusts and anticlines) are extracted from[71].In the map is showed the hypocenters location of the earthquakes occurred between 2003 and 2010 (ISIDe[74]) with the relative magnitude and hypocentral depth.Onshore capable faults from[46] and focal mechanisms from[45].Yellow dashed lines are the boundaries of three selected sectors (1.Corigliano; 2. Crati Delta; 3. Sibari-Villapiana) for Envisat and COSMO-SkyMed time series analysis (Figure6).(b) Sketch (based on the interpretation of the seismic profile F75-89 of[45]) representing the trigger of the flexural subsidence due to the overlap of the Amendolara Ridge back thrusts on Calabrid and Southern Apennines units.

Figure 6 .
Figure 6.Envisat and COSMO-SkyMed both ascending and descending time series (error = ˘1 mm) for the three selected sectors delimited by ITHACA [46] faults (location in Figure 5).Black dotted line represents the earthquakes (M = 4.5) occurred on 27 June 2006.

Remote Sens. 2015, 7 13 Figure 7 .
Figure 7. (a) Location of the used wells and seismic profiles.(b) Isopach map of the Plio-Quaternary succession with the distribution of the Envisat Vertical component.(c) Comparison between the thickness of the Plio-Quaternary succession and ground velocity, based on the Envisat Vertical component, along one of the analyzed profiles (its trace in Figure7a).The discontinuous lines for the InSAR profile are due to missing data in the ascending and/or descending original velocity maps due to coherence lack.

Figure 7 .Figure 8 .
Figure 7. (a) Location of the used wells and seismic profiles.(b) Isopach map of the Plio-Quaternary succession with the distribution of the Envisat Vertical component.(c) Comparison between the thickness of the Plio-Quaternary succession and ground velocity, based on the Envisat Vertical component, along one of the analyzed profiles (its trace in Figure7a).The discontinuous lines for the InSAR profile are due to missing data in the ascending and/or descending original velocity maps due to coherence lack.

Figure 9 .
Figure 9. (a) Envisat Vertical component overlapped to the Holocene deposits isopach map.(b) Scatter plot Envisat Vertical component vs. Holocene deposits thickness inferred from isopach map.(c) Comparison between Holocene sediments thickness and Envisat Vertical component in correspondence of each borehole.

Figure 9 .
Figure 9. (a) Envisat Vertical component overlapped to the Holocene deposits isopach map.(b) Scatter plot Envisat Vertical component vs. Holocene deposits thickness inferred from isopach map.(c) Comparison between Holocene sediments thickness and Envisat Vertical component in correspondence of each borehole.

Figure 10 .
Figure 10.(a) Analyzed boreholes between Saraceno and Satanasso alluvial fans (see Figure 1 for their location).(b) Correlation among Envisat ascending and (c) descending velocity and the thickness of fine-grained deposits (silt and clay).

Figure 10 .
Figure 10.(a) Analyzed boreholes between Saraceno and Satanasso alluvial fans (see Figure 1 for their location).(b) Correlation among Envisat ascending and (c) descending velocity and the thickness of fine-grained deposits (silt and clay).

Figure 11 .
Figure 11.Envisat ascending and descending time series (error = ±1 mm) comparison between urban and agricultural areas for three selected sectors (see Figure 4 for their location).

Figure 11 .
Figure 11.Envisat ascending and descending time series (error = ˘1 mm) comparison between urban and agricultural areas for three selected sectors (see Figure 4 for their location).

Figure 13 .
Figure 13.Envisat (both ascending and descending orbit) time series of the wells No. 2 and No. 4 (see Figure 1 for their location) recording respectively a watertable drop and rise.In table are reported the water table variation between 2002 and 2013 for each well.

Figure 13 .
Figure 13.Envisat (both ascending and descending orbit) time series of the wells No. 2 and No. 4 (see Figure 1 for their location) recording respectively a watertable drop and rise.In table are reported the water table variation between 2002 and 2013 for each well.

Table 1 .
Used image datasets.Ground resolution values are obtained after a multi-looked operation.
table drop and rise.In table are reported the water table variation between 2002 and 2013 for each well.
table drop and rise.In table are reported the water table variation between 2002 and 2013 for each well.