Different Ground Subsidence Contributions Revealed by Integrated Discussion of Sentinel-1 Datasets, Well Discharge, Stratigraphical and Geomorphological Data: The Case of the Gioia Tauro Coastal Plain (Southern Italy)

Abstract

Groundwater is the main water supply for agricultural and industrial needs in many coastal plains worldwide. Groundwater depletion often triggers land subsidence, which threatens manmade infrastructure and activities and aggravates other geohazards. We applied a multi-temporal interferometric synthetic aperture radar technique to Sentinel-1 datasets to detect ground motion in the Gioia Tauro plain (Calabria, Southern Italy) from 2018 to 2021. The InSAR data were analysed through the integrated use of groundwater head, stratigraphical and geomorphological data, and land use information to distinguish the potential subsidence divers. The results show that subsiding areas, with a mean rate of about 10 mm/yr, are in the middle of the plain, and their location is influenced by the spatial distribution of compressible sediments included in the shallow aquifer. Furthermore, the subsidence arrangement is spatially accordant with the main groundwater depression area, which can be ascribed to the ongoing and increasing water pumping for predominantly agricultural usage. We also observed that subsidence (up to 10 mm/yr) affects the western dock of the Gioia Tauro harbour, in front of which, in very shallow water, are two submarine canyon heads already affected by slides in the past.

1. Introduction

Coastal areas, representing the connection between terrestrial and marine environments, have huge social, economic, and biological value [1]. The economic value of goods and services provided by coastal ecosystems amounts to 77% of the total value given by all the rest [2]. Coastal zones are densely populated and host most of the world’s megacities [3]. Thanks to their features (fertility, water availability), coastal plains are sites of intensive agricultural development. Furthermore, coastal areas gather economic assets (e.g., tourism) and critical infrastructure, such as harbours (around 80% of the world’s trade is carried by maritime transport [4]).

Different hazards driven by natural processes and/or anthropogenic activities threaten coastal areas [5]. For this reason, the management of coastal zones and their sustainable development have been important topics for the last few decades [6,7,8,9,10,11]. The increasing anthropic pressure on coastal areas amplifies the negative impact of land use and hydrological changes, as well as groundwater overexploitation. The latter is often the driver of groundwater salinisation (e.g., [12,13]) and subsidence (e.g., [14,15]).

Over the last 20 years, land subsidence monitoring has been significantly improved thanks to the implementation of multi-temporal interferometric synthetic aperture radar (MT-InSAR), which permits the collection of time series of ground motion over wide areas, with millimetre accuracy [16,17,18,19,20]). In recent scientific literature, investigations into land subsidence using InSAR data is a popular topic, and it has been applied to areas worldwide that are characterised by different geological frameworks and various socio-economic backgrounds (e.g., [21,22,23,24]). The aim of this study was to analyse the ground motion of the Gioia Tauro coastal plain in the southeastern part of the Calabria region (Southern Italy) based on Sentinel-1 images and ancillary data (geological, hydrogeological, and land use information).

The Gioia Tauro plain is located along the western coast (facing the Tyrrhenian Sea) of the Calabria region (Figure 1) and is 500 km² wide. The plain is inhabited by about 160,000 people divided by towns, villages, and many rural locations. The area is characterised by a highly developed agricultural sector and two thermoelectric industries (waste-to-energy and thermoelectric power plants). Furthermore, the plain hosts two critical components of infrastructure: the Gioia Tauro harbour (the largest transhipment terminal in Italy and one of the most important container traffic hubs in the Mediterranean basin) and the A2 “Mediterranea” motorway (main road line of the region).

The ground motion of the Gioia Tauro plain has been examined in previous studies [25,26,27], which, using SAR data (ERS1/2 and Envisat), have identified subsidence with rates of up to 23 mm/yr from 1992 to 2006, without analysing the drivers. In our study, we monitored the plain subsidence until 2021, allowing us to observe the evolution of the phenomenon by comparison with the previous works. Moreover, in our research, underground and field data, as well as land use observation, were used in an integrated approach to investigate the potential subsidence drivers.

Figure 1. Simplified geological and structural map of the Gioia Tauro Plain; faults from [28] (CF—Cittanova fault; SLF—Sant’Eufemia–Laureana fault; PGF—Palmi–Gioia Tauro fault; NGFZ—Nicotera–Gioiosa fault zone; PLFZ—Palmi–Locri fault zone). Magnitude and epicentral location of earthquakes (from ISIDe) [29]; landslides (from PsAI-Rf) [30].

Figure 1. Simplified geological and structural map of the Gioia Tauro Plain; faults from [28] (CF—Cittanova fault; SLF—Sant’Eufemia–Laureana fault; PGF—Palmi–Gioia Tauro fault; NGFZ—Nicotera–Gioiosa fault zone; PLFZ—Palmi–Locri fault zone). Magnitude and epicentral location of earthquakes (from ISIDe) [29]; landslides (from PsAI-Rf) [30].

2. Study Area

The Gioia Tauro plain is a NNE–SSW-oriented half-graben of the Upper Pliocene–Pleistocene age [31] whose border consists of three main fault systems. These include (i) the NNE–SSW Cittanova fault along the inner and eastern edge, a high-angle normal and active fault; (ii) along the northern edge is the Nicotera–Gioiosa fault zone, striking NW–SE with right lateral kinematics; (iii) the NW–SE Palmi–Locri fault zone along the southern edge, characterised by mainly strike-slip kinematics [29]. In addition, other major faults also cross the plain—(i) the Palmi–Gioia Tauro fault near the coastline, a high-angle, NE–SW-striking, west-dipping fault, and another parallel east-dipping fault; (ii) the Sant’Eufemia-Laureana fault, a NE–SW-oriented and west-dipping normal fault (Figure 1) [28].

The eastern edge of the plain is marked by the juxtaposition between the crystalline basement (mainly made up of granitoid rocks [32]) and the Neogene–Quaternary sedimentary infill along the Cittanova fault [28,33]. The latter is an active fault as verified by paleoseismological studies which identified the fault ruptures in the 1783 catastrophic earthquake (Mw = 7.0) and in the 4th century AD [34,35]. Corresponding to its north-western edge, the plain is confined by the Capo Vaticano promontory consisting of granitoids and migmatites [32]. Along the southwestern edge, the plain is bounded by the Palmi high, made up of granitoids and gneiss [32], which has been uplifted and eastward tilted [33]. A structural high, with an approximate N–S orientation, is present along the western sector of the plain close to the coastline [28,36].

The Neogene–Quaternary sedimentary infill started with an upper Miocene siliciclastic and carbonate succession (cropping out mainly along the Palmi high), which covers the crystalline basement unconformably. Toward the top, the succession passes to the Trubi formation [37], consisting of Pliocene clayey and silty deposits a hundred metres thick [33], with a pumice-rich horizon in the upper part [38]. The succession continues with the Calcareniti di Vinco [37], late Pliocene–Pleistocene cross-bedded sands and calcarenites, the infill the Siderno paleo-strait [39]. Upward, along the eastern edge of the plain, gravelly and sandy alluvial fan deposits (Sintema di Taurianova of the Late Pleistocene age [37]) cover the Calcareniti di Vinco along an erosive contact. The rest of the Quaternary succession includes Late Pleistocene–Holocene terraced marine deposits made up of sands and gravels, Holocene aeolian sediments of ancient inland dunes, and recent alluvial and coastal deposits [37,40].

The Gioia Tauro groundwater system consists of three main hydrogeological units. From top to bottom, these include (i) the shallow and unconfined aquifer that includes all the Quaternary sediments; (ii) the aquitard, made up of the Pliocene clayey and silty deposits of the Trubi Formation; (iii) the deep aquifer, consisting of Late Miocene clastic deposits and limestone [41]. The plain is characterised by a mean precipitation rate of about 1000 mm [42].

Previous studies [25,26,27] have investigated the subsidence of the Gioia Tauro plain using data from the ERS 1/2 (1992–2001) and ENVISAT (2002–2006) satellites processed by the permanent scatterer interferometric synthetic aperture radar (PS-InSAR) technique [16]. They observed subsidence rates of up to 10–15 mm/yr between 1992 and 2006, recording, in the middle sector of the plain, peaks reaching 23 mm/yr in the period of 1992–2000. The two main subsidence bowls are recognised in the areas of the Rizziconi and Rosarno villages. The authors generically related the subsidence to the aquifer compaction due to the groundwater overexploitation (based only on the presence of numerous wells located in the plain). Computing an acceleration index between the two considered periods [27], they observed a general deceleration (about 65%) in the deformation rate, related to a supposed decrease in the groundwater pumping except for some areas in the southern part, which are, however, characterised by a low persistent scatterer density.

3. Data and Methods

Differential interferometric synthetic aperture radar (DInSAR) has been exploited with great success over the last twenty years to monitor ground displacement associated with, for example, earthquake deformation, volcanic movement, landslides, subsidence, and mining activities [19,43,44,45]. It shows limitations due to the spatial and temporal decorrelation and atmospheric phase disturbance. To address the disadvantages of DInSAR due to low temporal sampling, undesired atmospheric contribution, and temporal/spatial decorrelation, MT-InSAR techniques have been applied to overcome such limitations. MT-InSAR technique can be divided into two main groups, as follows:

(1) Persistent scatterer interferometry (PSI) technique—only one single master image is used [16,19,46]. The PSI technique is more suitable for monitoring local deformation phenomena over highly coherent structures, such as buildings, rail tracks, and bridges [16,46]. It uses only one master image to generate stacks of interferograms. The PSI method aims to identify coherent radar targets exhibiting high phase stability over the whole duration of the observations. These targets are only slightly affected by temporal and geometrical decorrelation, and often correspond to manufactured structures or bare soil.

(2) Small baseline subset (SBAS) approach [18,47] involves the usage of multi-master interferograms, and only the so-called distributed scatterer (DS) pixels are selected for further analysis. Interferometric pairs are chosen to minimise temporal and geometric decorrelation, allowing deformation time series to be retrieved for distributed scatterers, not being dominated by a single scatterer, and showing similar backscattering properties. Compared to single stable scatterers, an averaging operation is typically required to retain a good level of coherence, and thus, to achieve high measurement accuracies.

In the best case, both PSI and SBAS multi-temporal InSAR methodologies can reach accuracies as high as 1–2 mm/yr [48,49].

We adopted the multitemporal InSAR methodology named small baseline subset (SBAS, [18]), which is, nowadays, one of the most efficient and consolidated techniques to investigate many possible ground motions with millimetric accuracy. This technique was applied here to investigate local subsidence phenomena affecting the investigated area and mainly induced by water pumping activities. In detail, to map the ground displacement over the Gioia Tauro plain, 93 C-band Sentinel-1 (S1) images acquired between 11 September 2018 and 01 October 2021 along the ascending orbit (track 44), and 90 Sentinel-1 images acquired between 19 June 2018 and 22 May 2021 for the descending orbit (track 124), were analysed. The IW scenes were accessed through the GEP as Level 1 single-look complex (SLC) data (Figure 2), consisting of focused SAR products in zero-Doppler slant-range geometry, and georeferenced using satellite orbit and altitude data. Both datasets consisted of Sentinel-1 images acquired in the TOPSAR mode. Temporal and normal baselines thresholds were set to couple the SLC acquisitions to mitigate the spatial and temporal decorrelation phenomena (accordingly with the SBAS theory) and formed 262 and 254 interferometric pairs for the ascending and descending processing, respectively. Each SLC (for which the VV co-polarisation mode was selected) was oversampled by a factor of 20 in the azimuth and by 5 in the range direction, resulting in the application of a ground pixel spacing of 90 m to mitigate the decorrelation noise in the interferograms and generate the coherence maps. The topographic phase was derived using the 1 arc-second (30 m) DEM acquired from the Shuttle Radar Topography Mission (SRTM, [50]) and then removed from the formed differential interferograms. The interferograms were then noise-filtered using the Goldstein filter with a weight of 0.5. Temporal-spatial filtering was carried out to estimate the atmospheric phase components. It is generally accepted that both tropospheric turbulence and tropospheric stratification phases are correlated in space and not in time. Hence, a double filtering operation was carried out in both time and space, and 1200 m and 365 days were adopted as the values for the high and low pass filters. The SRTM DEM was also exploited to geocode the MT-InSAR results from the radar coordinate system into the Universal Transverse Mercator (UTM) coordinate system (Figure 2).

To assess the quality of the obtained outcomes, we estimated the precision velocity (V_prec) relative to the ascending and descending maps using the following formula [51]:

$$V_{prec}=\sqrt{\left(1-{\gamma }^2\right)/2{\gamma }^2}\times \left(\lambda /4\pi \right)$$

where γ and λ are the interferometric coherence and the wavelength, respectively. The higher this value is, the lower the measurement precision is.

Hence, we computed the statistics, obtaining mean values of 2.1 mm/yr and 1.7 mm/yr, and standard deviations of 0.5 and 0.8 for the ascending and descending cases, respectively, thus confirming the good reliability of the InSAR results.

Finally, thanks to the availability of both the ascending and descending ground velocities (Figure 2), the east–west and upward displacement components were retrieved (Figure 3), taking into account the relative incidence angles (θ_a, θ_d) of the ascending and descending processed S1 frames and the data according to [52].

From the two displacements (D_d and D_a are the descending and ascending velocity maps, respectively) measured along the LOS, we can recover the vertical and horizontal components of the displacement d_z and d_e, respectively. A straightforward geometric calculus leads to the following expressions:

$$d_e=\frac{D_d\cos {\theta }_a-D_a\cos {\theta }_d}{\sin ({\theta }_a+{\theta }_d)}$$

$$d_z=\frac{D_d\sin {\theta }_a+D_a\sin {\theta }_d}{\sin ({\theta }_a+{\theta }_d)}$$

4. Results

The distribution of the upward component (Figure 3a) shows the presence of subsiding areas mainly localised in the middle sector of the plain.

Southward, the largest subsidence bowl can be observed in the village of Rizziconi, with recorded upward velocities of up to −9 mm/yr. Northward of Rizziconi, in the village of Drosi, the velocity amounted to −11 mm/yr. Southward of Rizziconi, close to the Porto degli Ulivi shopping centre and the Cirello settlement, the maximum recorded velocity is −8 and −5 mm/yr.

Lower mean subsidence rates affect some rural locations, such as Russo (−3.3 mm/yr), Spina (−3 mm/yr), Marotta (−2.8 mm/yr), and Spartimento (−2.7 mm/yr), located in the middle sector of the plain.

Northward, upward velocities of up to −9 mm/yr characterise the San Ferdinando industrial area, while the velocities amounted to −9 and −5 mm/yr in the Fabiana and Ponte Annegato quarters (Rosarno village surroundings), respectively. Along the A2 motorway, a segment of 1.5 km over the viaduct (Mesima IV) on the Mesima River showed a velocity greater than −10 mm/yr.

Along the coastline, a subsidence area with a rate of up to −10 mm/yr was observed in correspondence to the western dock of the Gioia Tauro harbour.

An upward velocity of up to −5 mm/yr was observed in the town of Polistena (in the inner sector of the plain) close to the stadium area.

Scattered pixels and pixels located in riverbeds with high and negative upward velocity values were not considered.

Concerning the east component (Figure 3b), no significant area affected by horizontal ground motion is present. Most of the pixels characterised by high (positive and negative) east velocity were scattered and located in correspondence to anthropic activities (e.g., velocities of up to −6 mm/yr recorded northeast of the Polistena stadium fall on an aggregate quarry) and in the riverbeds. Along the A2 motorway, mean east velocities of 5 and −5.5 mm/yr were recorded in correspondence to the Mesima IV and Barletta viaducts, respectively. Between the villages of Rizziconi and Drosi, few pixels with velocities up to −6 and 8 mm/yr were present along slopes affected by localised and shallow landslides [30]. Close to the southern sector of the town of Cittanova, a negative east velocity of up to −5 mm/yr characterises the right bank of the Serra stream affected by a shallow landslide [30].

The main towns of the plain (Gioia Tauro, Rosarno, Nicotera, Taurianova, and Cittanova) are characterised by stability. For the town of Gioia Tauro, the absence of vertical ground motion was also confirmed by the time series of the GNSS station “Gioia”, which operated in the period of 2010–2019 and recorded a mean vertical velocity of −0.83 ± 0.67 mm/yr [53,54]. In correspondence to the GNSS station, the S1 upward velocity was equal to 1.4 mm/yr.

The use of ERS 1/2 (1992–2001) and Envisat (2002–2006) data, also used by previous studies [25,26,27], permits further remarks. The S1 data (2018–2021) confirm the presence of the main subsidence bowls in the middle sector of the plain. However, the mean subsidence rates (5–10 mm/yr) recorded in the S1 data (2018–2021) are lower than the values observed between 1992 and 2006. Furthermore, the surface affected by subsidence also seems to be reduced. The Sentinel data reveal a vertical land drop (up to −10 mm/yr) only for the Gioia Tauro harbour’s western dock, which was not recorded by the ERS and Envisat datasets.

The general decrease in the subsidence rates agrees with the deceleration between 1992–2001 and 2002–2006 in the deformation rate observed by [27], in which it was related to a decrease in the water requirement. Still, the reason for this deceleration cannot be ascribed to a reduction in groundwater pumping considering the increase in the water demand during the last decades due to drought periods [55,56] and land use change (see later). Therefore, in our opinion, the reason for the deceleration can be found in geotechnical causes such as the progress of the primary consolidation.

5. Discussion

In the following, we present our analysis of different potential contributions to the subsidence. Then, a detailed comparison of the InSAR results is presented with the available tectonic, stratigraphic, geomorphological, and hydrogeological information.

The plain is bounded and crossed by active (last activity during the Holocene) and mainly normal faults [35,57], but the low magnitude (maximum local magnitude of 2.8) of the earthquakes (Figure 1) that took place since 2017 [29] allowed us to exclude the presence of seismically-induced displacements.

We analysed the spatial pattern of the vertical velocity depending on the shallow aquifer thickness. The plain is characterised by a shallow and unconfined aquifer consisting of Late Pleistocene and Holocene marine and alluvial sediments, an aquitard made up of Pliocene clayey and silty rhythmites, and a deep aquifer made up of late Miocene deposits [41].

In detail, we plotted the S1 vertical velocity on an isopach map of the shallow aquifer (Figure 4). The thickness was estimated by means of the data of 131 boreholes from the literature [38,58,59,60,61], integrated by 133 vertical electrical soundings (VES), and inferred by SP26 [62] and field surveys.

The isopach map (Figure 4) shows that the shallow aquifer reaches its main thickness in the central plain, in the middle of the so-called Gioia Tauro extensional fault system [28]. We assessed the fine-graded index “If” as the relative thickness of the fine-grained sediments for each borehole or point as a percentage ratio of the total thickness of the silt and clay-type layers divided by the total shallow aquifer thickness.

Except for the A2 motorway and the Gioia Tauro harbour, the recognised subsidence bowls are localised in this central sector. However, close to the harbour and located westward of the PGF, the shallow aquifer thickness amounts to 80 m. The absence of subsidence can also be observed close to the town of Taurianova, where the shallow aquifer has its maximum thickness. Still, it does not include compressible sediments, showing a null or negligible If value, as inferred by the available boreholes (Figure 4a). Instead, the presence of compressible deposits characterises the shallow aquifer in the subsiding areas between Rizziconi and Cirello, and the industrial area of San Ferdinando. Compressible sediments are also included in the shallow aquifer in areas not affected by subsidence, such as the area around the village of Melicucco. We assessed the If for each borehole and the S1 upward velocity of the overlapping (n.33) or nearest (n.91) pixel analysing the correlation between the subsidence pattern and the thickness of compressible deposits. Considering the median value of the upward velocity, we observed a good correlation between the overlapping and nearest pixels, with R-squared values equal to 0.83 and 0.63, respectively (Figure 4b).

Close to the town of Rosarno, it is possible to observe different ground deformation patterns between the stable urban centre, built on Pleistocene sands and Holocene conglomerates, and the subsiding surroundings (Fabiana and Ponte Annegato quarters) located on the Mesima River (carrying mostly fine-grained sediments) flood plain (Figure 4c), which is a marshy area reclaimed during the first decades of the 1800s [63]. Another subsiding site located in correspondence to the Mesima flood plain is the viaduct of the A2 motorway.

Despite the main subsidence bowls being localised in the middle of the plain (located in the so-called Gioia Tauro extensional fault system), there is no linear correlation between the spatial subsidence pattern and the shallow aquifer thickness. This mismatch suggests that the isostatic sediment loading does not play an active role differently from other basins (e.g., [64,65,66]). Instead, a controlling factor of ground displacement seems to be represented by the thickness of fine-grained sediments (included in the shallow aquifer), as observed in other many Quaternary basins (e.g., [67,68,69,70]).

We investigated the correlation between groundwater exploitation and subsidence rates, comparing the vertical velocity distribution from the S1 data and the water table variation of the shallow aquifer (bearing all the groundwater abstractions for agricultural, industrial, and urban uses) in the time interval of the 1970s (data from [62]) to 2021 (new data collected from 88 wells) (Figure 5a). The main subsidence bowls are in areas characterised by a water level drop (−30/−35 m close to Rizziconi and Drosi, −5 m northward of Rosarno and up to −15 and −35 m in the industrial area of San Ferdinando and along the A2 motorway, respectively). The water table of the area between the Porto degli Ulivi shopping centre and Cirello oscillates between increase and decrease, while stability is recorded close to the Gioia Tauro harbour.

The widespread water depletion in the plain (Figure 5a) can be ascribed to the groundwater overexploitation testified by the presence of numerous water wells, which also feature 24-h pumping stations (Figure 6). The growing water requirement in the plain can be linked to the increase in industrial settlements and to the land use change that happened during the last two decades. In detail, the main change is due to the intensive and ongoing kiwifruit farming that requires a high-water supply (6000–10,000 m³/ha [71]) compared to the previous citrus and olive crops (4000–6000 m³/ha for citrus crops [72] and almost zero for olive [73], considering the mean annual rainfall of about 1000 mm in the plain [42]). For example, only in the area surrounding Rizziconi are 5.6 km² of kiwifruit orchards planted during the last decade (Figure 6).

We analysed multitemporal aerial orthophotos (2000 and 2008 available on [75]) and satellite images (2003, 2011, 2015, 2016, 2017, 2018, 2020; available on Google Earth platform) of the subsiding area to investigate the potential urban area growth effects on subsidence. The only significant variation was observed in correspondence to the Porto degli Ulivi shopping centre, which was inaugurated in 2007.

We analysed the LOS S1 time series (reporting the ground deformation trend in time) for both the ascending and descending orbits of the subsiding areas (surrounding Rizziconi and Rosarno and the industrial area of San Ferdinando) located in the middle sector of the plain. The time series from both orbits (Figure 7) show a general trend characterised by an increase in negative displacement, shown as absolute values. Furthermore, some of them, both ascending and descending, seem to display a seasonal trend with positive and negative displacements with the wet (October–April) and dry (May–September) seasons. In detail, the top five time series (Figure 7) show oscillations of displacement values characterised by a positive trend (values above the regression line) during the wet season and by a negative trend (values below the regression line) during the dry months.

This alternation also coincides with the aquifer recharge and dry periods (corresponding to the time of maximum groundwater exploitation for agricultural use), characterised by water table variations (referring to our monitoring in 2021) between −0.5 and −3 m.

In correspondence to the subsiding sector of the A2 motorway, the time series of both orbits (Figure 6) are characterised by a clear linear trend of ground motion, with values of cumulative displacements of up to −7.5 and −5.5 cm for the ascending and descending tracks, respectively.

Our results show an important relationship between the piezometric level changes and subsidence (Figure 5a). For this reason, we observed that the decrease in the groundwater level is the main trigger and control factor of temporal and spatial evolution of subsidence over the plain, as observed in other case studies worldwide, such as the Beijing plain (China) [76], the Konya plain (Turkey) [77], the Alto Guadalentín aquifer (Spain) [78], and the tectonic valleys of Central Mexico [79].

The seasonal trend of ground displacements (shown by the time series) that are related to piezometric level fluctuations is in agreement with observations made by other authors for test sites located, for example, in the United States [80,81], Mexico [82], and Italy [70]. Furthermore, considering that the shallow aquifer mainly includes coarse-grained sediments with high permeability [41], the seasonal trend of ground deformations reflects the compaction and dilation of the aquifer system in response to the seasonal fluctuations of the groundwater level linked to natural recharge and exploitation.

Concerning the subsiding sector of the A2 motorway in correspondence to the Mesima IV viaduct, the linear trend of ground motion could be ascribed to the permanent load due to the road infrastructure and traffic, along with the piezometric level drop and the presence of soft sediments in the Mesima flood plain, and the construction characteristics of the consolidation of the road embankment. This motorway sector is also affected by eastward movement (contrary to the flow direction of the Mesima River) of up to 5 mm/yr; in situ monitoring and information on the construction would be needed to understand if the horizontal motion concerns the deck and/or substructure.

A particular analysis was reserved for the subsidence (up to −10 mm/yr) affecting the western dock of the Gioia Tauro harbour, where the Sentinel time series for both the ascending and descending orbits (Figure 7) show non-linear trends (accelerations and decelerations of the movements) of ground motion. In this area, in addition to the usual potential subsidence triggers (e.g., groundwater exploitation), it is necessary to consider submarine morphological evolution. In front of the harbour mouth and western dock, the two headwalls of the Gioia Canyon are present, arriving at shallow water depths < 10 m (Figure 8) and flowing into the Stromboli Canyon 20 km away from the coast [83]. The hazard due to the slope failures at the submarine canyon headwall is high considering the nearness to the coast and the retrogressive erosion characterising the canyon heads [83]. A slope failure involving 5 million m³ of material occurred at the head of the Gioia Canyon on 12 July 1977 during the construction of the harbour; this event was not triggered by an earthquake or flood [83,84,85].

We suggest that the retrogressive activity of the submarine canyon heads, potentially driving or aggravating deep-seated deformations of the coastal areas, contributes to the local subsidence.

Subsidence triggers related to seaward deformations close to the submarine canyon heads were also proposed for the Vulcano [86] and Lipari [83] islands of the Aeolian archipelago (southern Italy) and the coastal sector between Cirò Marina and Punta Alice in the Calabria region [87].

Figure 8. Morphobathymetry (MAGIC project [88]) in front of the Gioia Tauro harbour showing the two heads of the Gioia Canyon in very shallow water. The map also shows the position of the 1977 landslide as reported by [84].

Figure 8. Morphobathymetry (MAGIC project [88]) in front of the Gioia Tauro harbour showing the two heads of the Gioia Canyon in very shallow water. The map also shows the position of the 1977 landslide as reported by [84].

6. Concluding Remarks

InSAR results from data acquired by the S1 mission spanning from 2017 to 2021 allowed us to detect subsiding areas in the Gioia Tauro Plain. The primary subsidence bowls, with a mean rate of ~10 mm/yr, are in the middle sector of the plain, while along the coastline, the only subsidence (up to 10 mm/yr) area is the western dock of the Gioia Tauro harbour.

Except for the harbour, the spatial pattern of subsidence depicted by the Sentinel-1 data is in agreement with that of the ERS 1/2 (1992–2001) and Envisat (2002–2006) data previously analysed by other authors [25,26,27]. Furthermore, the general deceleration of ground motion observed between 1992–2001 and 2002–2006 [27] is confirmed by the S1 data.

After locating the subsiding area, we compared the InSAR results with ancillary data (geological, hydrogeological, and land use information) to identify the driving forces of the subsidence. Our analysis shows that subsidence mainly occurs where compressible sediments are included in the shallow phreatic aquifer. However, the groundwater level drop (up to 40 m) plays the main role, which characterises the middle sector of the plain, as documented by the comparison between the historical (1970s) and present 2021 piezometric levels. This water level drop is due to the ongoing overpumping of groundwater for agricultural use (the most prevalent), household needs, and thermoelectricity production. In addition, during the last decades, the increased water requirement was due to a widespread land use change consisting of extended kiwifruit planting (the water demand of kiwifruit is higher compared to previously farmedcitrus and olive crops). The deceleration clearly shows that not only groundwater exploitation is the driving force (as reported in previous works), as groundwater exploitation has definitely not been reduced, while subsidence is decelerating. For this reason, the deceleration can be found in geotechnical causes.

The ongoing aquifer compaction started at least three decades ago, as testified by the subsidence recorded by the InSAR data of this and previous works, producing an irreversible loss of aquifer porosity and storativity. In the long term, this effect can reduce management margins and aquifer performance with consequent damage to the communities depending on the water resource and to environmental equilibrium. Furthermore, the enduring subsidence can produce severe damages to civil structures and infrastructure, as observed in other subsiding areas worldwide (e.g., [79,89,90]). In our opinion, the managing authority must plan the license of new water wells, taking into account this issue and evaluating alternative water supplies to satisfy the growing water requirement. By reducing the groundwater pumping by wells, the mitigation of subsidence will occur, as observed, for example, in Queretaro City (Mexico) after a long pipeline was built, bringing surface water allowing for half the number of water wells [82].

We also observed subsidence (up to −10 mm/yr) in correspondence to the Gioia Tauro harbour’s western dock. We suggest that the retrogressive activity of the Gioia Canyon headwalls (extending to water depths < 10 m) are the driving force of the ground motion. Local and multi-sensors monitoring the dock displacements and integrated by observing the evolution of the submarine canyon headwalls through multi-temporal bathymetric surveys can provide the necessary information to understand the phenomenon better and to anticipate the occurrence of hazards.

Our future aims consist of more detailed analyses of the identified subsidence areas using long-term time series acquired by different sensors, geotechnical data, and permanent groundwater level monitoring. The final challenge is the disentanglement of different subsidence triggers and driving forces and the estimation of their relative contributions.