Sentinel-1 BigData Processing with P-SBAS InSAR in the Geohazards Exploitation Platform: An Experiment on Coastal Land Subsidence and Landslides in Italy

: The growing volume of synthetic aperture radar (SAR) imagery acquired by satellite constellations creates novel opportunities and opens new challenges for interferometric SAR (InSAR) applications to observe Earth’s surface processes and geohazards. In this paper, the Parallel Small BAseline Subset (P-SBAS) advanced InSAR processing chain running on the Geohazards Exploitation Platform (GEP) is trialed to process two unprecedentedly big stacks of Copernicus Sentinel-1 C-band SAR images acquired in 2014–2020 over a coastal study area in southern Italy, including 296 and 283 scenes in ascending and descending mode, respectively. Each stack was processed in the GEP in less than 3 days, from input SAR data retrieval via repositories, up to generation of the output P-SBAS datasets of coherent targets and their displacement histories. Use-cases of long-term monitoring of land subsidence at the Capo Colonna promontory (up − 2.3 cm/year vertical and − 1.0 cm/year east–west rate), slow-moving landslides and erosion landforms, and deformation at modern coastal protection infrastructure in the city of Crotone are used to: (i) showcase the type and precision of deformation products outputting from P-SBAS processing of big data, and the derivable key information to support value-adding and geological interpretation; and (ii) discuss potential and challenges of big data processing using cloud/grid infrastructure.


Introduction
The last three decades have witnessed the increasing and widespread use of satellite synthetic aperture radar (SAR) imagery to retrieve key information on surface deformation processes affecting the Earth's surface using two-and three-pass differential interferometric SAR (InSAR) (e.g., [1][2][3]), and multi-pass/multi-temporal InSAR methods, such as persistent scatterers interferometry (PSI) (e.g., [4,5]) and Small BAseline Subset (SBAS) (e.g., [6,7]). Many SAR constellations have been developed since the 1990s and, as part of their background observation scenarios, have been tasked to acquire images regularly over major cities, targets with economic and/or strategic relevance (e.g., key infrastructure, polar regions, cultural heritage sites), regions affected by specific Earth hazards (e.g., tectonic zones and volcanoes), and entire nations or continents (e.g., [8,9]). This means that the volume of SAR data collected to date, archived, and available for InSAR applications by far falls within the "big data" remit for such type of imagery (typically, 15-20 SAR scenes are needed as minimum input of an advanced InSAR analysis, e.g., [5], and usually a few to several tens of scenes are used in most InSAR applications), thus creating novel opportunities and opening new challenges for Earth sciences and observation (e.g., [10,11]). For instance, since the beginning of its operations in late 2014, the Copernicus Sentinel-1 SAR constellation acquired over each area of Europe around 300 passes for each geometry, i.e., ascending and descending, with an initial repeat pass of 12 days with Sentinel-1A only (~180 scenes in 6 years), which was improved to 6 days in late 2016 thanks to the addition of Sentinel-1B to the constellation (another~120 scenes in 4 years). This can be verified by consulting the ESA Copernicus Open Access Hub [12].
Thus, researchers willing to undertake multi-temporal InSAR studies of long-term (e.g., more than 3-5 years) dynamics of deformation processes face two main challenges. First, they need to find efficient combinations of computing infrastructure, algorithms, and software capable to run and optimize the big SAR data processing workflow. Second, they need to use tools to handle and analyze the outputs of such processing, which (proportionally to the processed area extent and its land cover) may consist of hundreds of thousands or millions of persistent or coherent targets, each equipped with information of geolocation, quality, annual displacement velocity, and a full deformation record, the length of which depends on the number of input SAR scenes processed.
Parallelized InSAR processing chains running on cloud and grid computing infrastructure help to address the first challenge. They relieve researchers of demanding tasks for big data download and processing in local PCs, and drastically reduce processing times (e.g., [13,14]). On the other side, InSAR dataset screening approaches and statistical methods to detect outliers, trend changes, accelerations/decelerations, and other deviations from the expected deformation behavior contribute to addressing the second challenge [15][16][17][18].
In this respect, the present study has three main goals: 1.
First, to demonstrate that through parallelized InSAR processing chains running on data exploitation platforms, it is possible to process big SAR time series made of an unprecedented number of Sentinel-1 images, i.e., in the order of hundreds per single observation geometry.

2.
Second, to showcase that big InSAR datasets can be generated across geographic areas of either urban or rural land cover (or a mixture of them), and later be post-processed to derive value-added products. These InSAR datasets are as "big" in the time-domain as the amount of input SAR images, and therefore as the number of observations composing the deformation time series of each persistent or coherent target obtained at the end of the multi-interferogram processing. Furthermore, they are "big" in the space-domain on the basis of the total number of persistent or coherent targets themselves that the interferogram network created during InSAR processing has allowed to identify over the long time span investigated (despite the relatively small region of interest of this work). 3.
Third, to prove that such big InSAR datasets can be handled through post-processing and data integration methodologies that have been established by the geohazard community in the last decade, and be transformed into higher-level geospatial information to characterize geological hazards of potential concern of local stakeholders. Without claiming to fully develop specific user-oriented or user-driven applications, this paper does not restrict to the presentation of InSAR results only (which is the dominant feature of most of the published literature on InSAR processing using cloud/grid computing), but also encompasses selected use-cases of geological interpretation.
To address the first two aims, in this study, six full years of Sentinel-1 SAR ascending and descending mode imagery compose the big data volume used as input. Multi-temporal InSAR processing was carried out using the parallelized implementation of the SBAS differential InSAR approach, also referred to as Parallel-SBAS (P-SBAS) [19], which was developed at the Institute for Electromagnetic Sensing of the Environment (IREA) of the National Research Council (CNR) of Italy and integrated in the European Space Agency (ESA)'s Geohazards Exploitation Platform (GEP). The latter is a cloud-based environment sourced with data (both optical and SAR), tools, and processing services, aimed to support satellite Earth observation (EO) exploitation for geohazards [20]. Developed as an ESA Thematic Exploitation Platform (TEP) in 2015, and conducted by a consortium led by Terradue Srl, the GEP follows the SuperSites Exploitation Platform (SSEP), which was originally initiated in the context of the Geohazard Supersites and Natural Laboratories initiative. Over the years, the platform has expanded to address broader objectives of the geohazards community, including several activities and projects of the Committee on Earth Observation Satellites (CEOS) Working Group on Disasters [21]. The GEP is now used by a large community of researchers and stakeholders working in the field of EO, with a focus on mapping hazard-prone areas and monitoring terrain deformation. Information on how to access the GEP and exploit its on-demand and systematic services is available on the dedicated website (https://geohazards-tep.eu), along with a description of all thematic apps, services, datasets, and other resources.
To the best of our knowledge, the present big data processing experiment with P-SBAS is the one exploiting the longest-ever Sentinel-1 SAR interferometric stacks used so far for a single study site, not only among the studies employing the P-SBAS method, but also globally, as the longest Sentinel-1 stacks published in the literature consisted of no more than 140-190 acquisitions each (e.g., [22][23][24]) or up to 240-260 (e.g., [25]). Therefore, in terms of the size of the input SAR data stack processed, this study aspires to act as a pioneering proof-of-concept of big SAR data analysis with multi-interferogram algorithms, and in particular the P-SBAS implementation.
Compared to the more paroxystic and catastrophic hazards that are commonly observed with the GEP (e.g., floods and tectonic deformation), this study focuses on more "silent" and subtle hazards, i.e., slow-moving landslides and coastal subsidence, whose regular monitoring and early-stage characterization could provide evidence-based input into land management, land-use planning, and environmental assessments. For this purpose, a coastal study area in southern Italy ( Figure 1) was selected to include the city of Crotone and its province, the Capo Colonna promontory and the namesake archaeological site in Calabria Region, where ground instability, including both natural geological processes and human-induced hazards, have been recorded (e.g., [26][27][28]). Instances of surface deformation matching with known landslides, long-term subsidence, and modern coastal protection infrastructure are used to showcase the type and precision of deformation products that can be generated with P-SBAS processing of Sentinel-1 big data using the GEP, and also the key information retrievable to support the geological interpretation and analysis of the process under investigation. This includes value-adding through the combination of multiple P-SBAS datasets for the retrieval of the vertical and east-west deformation fields (Section 5.1), and discussion of the geo-information and knowledge enhancement through ground deformation datasets outputted from a SAR big data analysis (Section 5.2). topography from NASA's SRTM digital surface model, with indication of hydrocarbon-production infrastructure and concession areas (data source: [29]); and (c) simplified geology (modified after [30]). D.C 1.AG, D.C 2.AG, D.C 4.AG, and F.C 1.AG are the unique identifiers of the concession areas.

Study Area
The study area in Calabria Region, southern Italy, extends 300 km 2 across the Crotone Basin (Figure 1), which lies on the Ionian side of the Calabrian Arc, i.e., one of the highly active Mediterranean tectonic regions. From a geological point of view, the basin is generally considered a fore-arc basin [31]. It started opening between Serravallian and Tortonian, on top of the SE verging Calabrian accretionary wedge, the latter being grown along the Neogene subduction zone between the sinking Ionian lithosphere and the overriding European plate [32]. Two NW-trending sinistral shear zones laterally bound the basin [31]. According to the model proposed by Van Dijk [33], such basement wrench faults caused the partitioning of the basin into separate sub-basins subject to differential subsidence and mutual displacements [34].
Along the coasts of the Crotone province, several marine terraces nowadays provide the evidence of the strong uplift that affected the Calabrian Arc since the Middle Pleistocene. At these terraces, a transgressive erosional surface can be recognized, above which shallowmarine sediments unconformably overlie a Plio-Pleistocene slope succession known as the Cutro Clay [31]. In absence of a general consensus on their number and identification [31], and following the identification made by Zecchin et al. [35], these marine terraces are referred to as (Figure 1c): the older Cutro terrace (1st order, Middle Pleistocene) and the Late Pleistocene Campolongo (2nd order), Capo Cimiti (3rd order), Capo Rizzuto (4th order) and Le Castella (5th order) terraces [30,36]. Being placed at progressively lower elevations (around 180-200 m a.s.l.) and dislocated by Quaternary faults [31], these terraces are a distinctive feature of the local topography ( Figure 1b).
Of these, it is worth focusing on the terrace a few kilometers SE from Crotone that occupies the 3 km EW-elongated promontory of Capo Colonna (Figure 1b), i.e., the easternmost piece of land in the study area, projecting into the Ionian Sea with an eastward dip of 1 • [31]. A topographic gradient characterizes the promontory, from~50 m a.s.l. elevation at the N-trending landward border representing a paleo-coastline, up to 10-15 m a.s.l. at the coastal cliffs [31]. The internal margin of the terrace is covered by a colluvial layer, while seaward up to 10 m thick carbonate and siliciclastic sediments are found. This complex sedimentary body exhibits a non-homogeneous behavior, from both hydrogeological and mechanical points of view [37].
While in other terraces the steep cliffs are carved in terrace deposits, at Capo Colonna they are carved directly into older units [38] and are subject to severe marine erosion due to strong wind-driven waves undercutting the newly exposed base of the cliffs. Large slabs of Pleistocene fossiliferous calcarenite consequentially fall down [27], thus contributing to severe scarp retreat that threatens the structural stability of properties and the preservation of renown archaeological monuments located close to the cliff edge [27,37,39].
Coastal erosion has been affecting this area for centuries, together with subsidence and tectonic-induced submergence of sites that were once sub-aerially exposed along the coast. Zecchin [40] has demonstrated the relationship between such fault-controlled subsidence and preservation of shallow-marine small-scale cycles visible in the northern Crotone Basin. More south in our study area, Stanley and Bernasconi [27] and the studies that they reviewed have demonstrated that the positions of modern Calabrian shorelines considerably differ from the former in Greek and early Roman times, with a landward retreat in the order of 200 m. Several sites submerged due to a number of tectonic pulses over time. Excluding world (eustatic) sea-level rise, Stanley and Bernasconi [27] estimate that during the past~3000 years in the late Holocene, subsidence has ranged from~0.5 mm/year to the remarkably high value of~4 mm/year. At present, other factors that may cause land subsidence and surface ground deformation across the region are not only hotspots of groundwater pumping for local agricultural and industrial activities, but also hydrocarbon exploitation [41].
The Crotone Basin and its offshore expansion have been extensively investigated for oil exploration [32]. Natural gas seepages were discovered in Cutro and Isola Capo Rizzuto, while a group of fields with thermogenic gas is located offshore of Capo Colonna ( Figure 1b) [42]. The Luna-Hera Lacinia-Linda field (1971, gas 370-400 million barrels of oil equivalent) is indeed among the several offshore fields discovered in Italy during the decade 1965-1974, and has been attributed to the Tertiary Supercycle-Oligocene-Miocene thermogenic gas petroleum system family [43]. Figure 1b shows the location of hydrocarbon production infrastructure and concession areas as per the Italian Ministry of Economic Development's database [29]. For the purposes of interpretation of the ground motions found in the Sentinel-1 data InSAR processing (see Section 4), it is worth focusing on the concession area called "Capo Colonne," which extends nearly 65 km 2 over the southern part of the Crotone province and encompasses the namesake promontory. Granted on September 1997 for 30 years, and once renewed in 2007 for additional 10 years, "Capo Colonne" has recently been operated by Ionica Gas, and then by the national hydrocarbons authority ENI (Ente Nazionale Idrocarburi), and includes the Hera Lacinia collection and treatment center (Figure 1b), as well as the production well Hera Lacinia 001 [29].
The potential impact of such exploitation activities on ground stability and local anthropogenic assets has been assessed based on an ad hoc geodetic-geophysical monitoring network installed in 2005 in the Archaeological Park of Capo Colonna, in the framework of a collaboration between ENI Spa E&P Division and the Superintendence for Archaeological Heritage of Calabria [37]. This integrated system included satellite SAR surveys among the shallow monitoring techniques employed by ENI in their typical subsidence monitoring networks [44]. SAR was also used in research papers for high-resolution [45] and multi-temporal [46,47] assessment of land subsidence in the last 25 years.
InSAR was helpful to investigate landslides that are distributed across the Crotone province and update their inventory [48]. Through integration with seismic profiles, morpho-bathymetric and GPS data, InSAR was brought among the evidence to support the discussion on the existence of a large-scale gravitational megalandslide off the Crotone Basin [28].
Limiting to inland susceptibility to landslides, which is one of the geological focuses of the present study, it is worth recalling that in the central and southern parts of the Crotone province hill-slopes carved into clayey and marly-clayey deposits are subject to slope-wash processes (sheet, rill, and gully flow), and lead to the development of erosion landforms called "calanchi" (i.e., a type of badland formed in clayey bedrock) [49]. These landforms are found in the coastal area south of Crotone; for instance, between the settlements of Vrica and Semaforo, and more inland toward Farina IV and Trafinello ( Figure 2).

Input SAR Data
A total of 579 Sentinel-1 SAR scenes were used for the analysis. These were acquired by the radar sensors operating in C-band (5.405 GHz central frequency f, and 5.547 cm wavelength λ) onboard the two satellites Sentinel-1A (launched on 3 April 2014) and Sentinel-1B (launched on 25 April 2016), which are phased at 180 • to each other along the same near polar, sun-synchronous orbit at~693 km altitude [50].
The data were acquired in right-looking mode over a 6 year-long period between October 2014 and October 2020, along ascending track T146 (296 scenes: 15 October 2014-13 October 2020, acquired at 18:48 local time) and descending track T51 (283 scenes: 21 October 2014-13 October 2020, at 06:55 local time). The spatial footprints of the two data stacks are shown in Figure 1a and, with zoomed view, in Figure 3. As per the mission observation scenario over Europe [9], until September 2016, the site revisit at the study area was of 12 days (for each geometry), while an improvement to 6 days was later achieved thanks to the addition of Sentinel-1B. This means that each data stack consists of an image every 12 days in the period October 2014-September 2016, and every 6 days in September 2016-October 2020. The frequency distribution of the available scenes along both tracks is shown in the time-bar in Figure 3. All scenes acquired in the period October 2014-October 2020 were used for the analysis. The employed acquisition mode was the Interferometric Wide (IW) swath, which is based on the Terrain Observation with Progressive Scans (TOPS) ScanSAR mode [51]. The IW swath width was 250 km, and the incidence angle θ of the line-of-sight (LOS) varied between 31 • and 46 • from the near to the far range (i.e.,~39 • at scene center). The pixel spacing was 2.3 m in slant range and 14.1 m azimuth, resulting in ground range and azimuth resolution of 5 m and 20 m (single-look), respectively.
The IW scenes were accessed through the GEP as Level-1 Single Look Complex (SLC) data ( Figure 3), consisting of focused SAR products in zero-Doppler slant-range geometry, georeferenced using satellite orbit and attitude data.
While Sentinel-1 IW data were acquired with dual polarization (co-polarized VV and cross-polarized VH), for this study the VV imagery was selected to run the P-SBAS processing, with the aim to guarantee continuity with previous co-polarization InSAR studies over the same region (see Section 5.1.1). Future research on this area will also attempt the exploitation of the VH data, both for comparison and in combination with VV. While the use of the VH data alone provides lower performance in terms of candidates and final persistent targets identification due to the generally lower amplitude of cross-polarized imagery (e.g., [52]), their combined exploitation with VV data through optimization methods proved capable to enhance the performance (e.g., [52,53]). This could be a way forward for future research specifically aiming to compare the InSAR processing performance and accuracy at different polarizations.
In this framework, the published literature shows a yet relatively limited number of applications that aimed to exploit this service for geological investigations (e.g., [23,24,61,63,64]), since the release of the P-SBAS Sentinel-1 on-demand service within the GEP in early 2019. Therefore, this study also contributes to this research direction.
The processing in GEP began with the retrieval of the input SLC data from the image archives, which is mainly based on the repositories of the ESA Open Access Hub, ONDA DIAS, and CREODIAS. Each SLC data stack was co-registered at the single-burst level ensuring very high coregistration accuracy (in the order of 1/1000 of the azimuth pixel size), as required for TOPS data due to the large along-track Doppler centroid variations [73]. Small temporal and perpendicular baseline pairs were selected to form the interferograms, with the aim to reduce interferometric phase decorrelation [74] by limiting the maximum number of consecutive acquisitions to be used to form the pairs [54], and by exploiting the narrow orbital tube (i.e., 100 m diameter during nominal mode operations; [50]).
A complex multi-looking operation in azimuth and range (5 by 20) was applied to mitigate the decorrelation noise in the interferograms and generate the coherence maps, thus reaching ground pixel spacing of~90 m, as per the default parameter setting of the current P-SBAS service version in the GEP. Topographic phase components were derived and subtracted based on NASA's 1 arcsec (i.e., 30 m) resolution Shuttle Radar Topography Mission (SRTM) Digital Surface Model (DSM). Multi-looked burst interferograms were mosaicked to sub-swaths, and then the two sub-swaths covering areas of land were mosaicked together to reach two-thirds of the swath. i.e.,~170 km (see stack coverage in Figure 3). The interferograms were then noise-filtered using the Goldstein filter with 0.5 weight [75], and optimized triangular networks of interferograms were identified [76], including 861 interferograms for the ascending stack, and 828 for the descending stack. Temporal and spatial phase-unwrapping was implemented with the extended minimum cost flow (EMCF) algorithm [77]. Across the large expanses of rural and vegetated land in the northwestern sector of the processed region, the phase-unwrapping step was more challenging, as typically happens across these land cover types, where errors are more likely to occur due to widespread phase decorrelation.
Least square solution retrieval was used for coherent targets using the singular value decomposition (SVD) method, and residual orbital ramps were estimated (as the linear components of phase artifacts; [78]) and then removed from the wrapped interferograms [19]. A second iteration of temporal and spatial phase-unwrapping was implemented for the orbit error-free interferograms, followed by SVD, as per the P-SBAS Sentinel-1 workflow [54]. Temporal coherence was computed, and the threshold on the minimum allowed temporal coherence to select coherent targets was set at 0.85. Atmospheric phase components were identified and removed. The output datasets including geolocation (i.e., latitude and longitude in the WGS84 system, and elevation above the reference ellipsoid), annual LOS velocity, and date-by-date displacement histories (time series) for each coherent target were exported in comma-separated values (.csv) ASCII format, according to the European Plate Observing System-Implementation Phase (EPOS-IP) project specifications. The corresponding metadata, key properties, and low-resolution quick-looks of the retrieved annual The stable point for both ascending and descending P-SBAS processing was set at the same location within the city of Crotone, at 39.078 • N, 17.120 • E, and the annual LOS velocity values and time series were referenced accordingly. The use of a common reference point allowed internal calibration of the two output datasets, and enabled their combination as described in Section 3.3.

Vertical and East-West Deformation Field Estimation
As increasingly done in the Sentinel-1 InSAR literature (e.g., [24,61,79,80]), the ascending and descending datasets were combined to retrieve the vertical V U and east-west V E velocity field. A fishnet of 90 m square elements was used to resample the point-wise datasets to a regular grid, and to link the output datasets into a unique digital layer. When both P-SBAS outputs (ascending and descending) were available at the same location i, the combination was achieved under the assumption of negligible north-south velocity, V N = 0. This assumption is typically employed in InSAR studies to account for the relatively poor visibility of north-south horizontal movements that the sensor LOS can achieve (e.g., [81,82]).
Given the known values of the LOS deformation velocity in the ascending (V A ) and descending (V D ) mode at each location i, the V U and V E were estimated as follows (e.g., [65]): where E, N, and U are the east-west, north-south and vertical directional cosines or unit vectors that define the 3D orientation of the Sentinel-1 IW LOS, and depend on the incidence angle θ, and the ground track or heading angle α of the satellite flight direction (i.e.,~345-350 • and~190-195 • for the ascending and descending modes, respectively). Their values can be estimated at each location as: E i = −cosα i *sinθ i , N i = sinα i *sinθ i , and U i = cosθ i , and accounting for θ variations across the IW swath. When only the ascending (or descending) mode dataset was available at location i, the additional assumption of null east-west velocity was required, V E = 0, and the vertical velocity was then estimated by simply dividing the LOS velocity value by U i .

Deformation along the Steepest Slope Direction
In hilly and mountains sectors of the study area that could be affected by mass movements, a simplified geomorphologic scheme was also considered (e.g., [83]), based on which deformation is more likely to occur along the steepest slope direction. This is true for translational slides, while it may not be suitable for the whole body of rotational landslides, topples and other types of mass movements, including erosion and calanchi landforms such as those present in some sectors of the study area.
Under the assumption that displacements occur along the steepest slope direction, the LOS estimates were projected onto it by considering the local values of terrain slope (β) and aspect (γ) as estimated from the SRTM DSM. These were used to identify the orientation of the steepest slope and its directional cosines at each location i, i.e., E Si = sinγ i *cosβ i , N Si = cosγ i *cosβ i , and U Si = −sinβ i .
The velocity along such direction (V S ) was then estimated as follows (e.g., [84]):

Results
The P-SBAS jobs for the whole processing area lasted 2.5 and 2.9 days for the Sentinel-1 IW ascending and descending datasets displayed in Figure 3, respectively. Therefore, these were the total times to process the stacks at full frame scale. In this particular experiment, two sub-swaths were processed for each geometry, given that the third covered only open waters (see Figure 3). The duration also includes the time required for sourcing the Sentinel-1 images to input into the processor from ESA's repositories. Therefore, based on the above numbers, it can be well appreciated how advantageous and well-performing the P-SBAS parallel processing and GEP infrastructure were to ingest such high number of images, co-register, multi-look, build highly interconnected interferogram networks, and analyze them to extract thousands of coherent targets and their associated deformation time series.
Within the whole processed area (Figure 4), the retrieved number of targets using equal temporal coherence threshold of 0.85 was of 47,400 across 3800 km 2 land area for the ascending mode (~12.5 targets/km 2 ), and 12,600 across 2000 km 2 for the descending mode (~6.3 targets/km 2 ). The distribution and density of coherent targets reflected the land cover of the region ( Figure 5), with vast expanses of forested areas of the Sila Massif, rural, and agricultural lands in the Crotone basin (where the step of phase-unwrapping was also more challenging), determining smaller densities of targets (e.g., 30-40 targets/km 2 ; Figure 5b) with respect to densely covered urbanized sectors, coastal infrastructure and shipyards (e.g., 180-200 targets/km 2 ), e.g., the city of Crotone ( Figure 6). Coherent targets were also found at four offshore platforms (see Section 2; Figure 1b) thanks to their good coherence across the differential interferograms.  LOS velocity in 2014-2020 spanned between −3.7 and +1.6 cm/year in the ascending dataset, and between −2.0 and +1.9 cm/year in the descending dataset. Negative velocity values indicated movement away from the satellite sensor (yellow to dark red points in Figure 6a,b), while positive values indicated movement toward the sensor (light to dark blue points). Standard deviations with respect to a linear deformation model were 0.01 cm/year on average across the dataset of coherent targets, for both the ascending and descending mode. Maximum values of 0.09 and 0.05 cm/year were recorded for the ascending and descending datasets, respectively. The deformation histories were characterized by 0.3 cm standard errors on average across both datasets. These values provided an estimate of the precision of the P-SBAS datasets, and guided through the setting of the velocity range to identify coherent targets that could be considered stable (green points in Figure 6a,b). Standard errors for the time series were influenced by the presence of nonlinearities in the long-term deformation trend, as discussed in Section 5.
The vertical velocity field after combination of the two datasets showed widespread presence of negative velocities indicating subsidence, with observed rates as high as −4.7 cm/year. The largest subsiding zone was located at Capo Colonna promontory, where the rates reached −2.3 cm/year (Figure 7a). Along the Ionian coast and involving the harbor infrastructure another subsidence pattern can be recognized (Figure 8a). Within an active landslide area south of Crotone, the vertical rates also reached nearly −3.0 cm/year across a large sector (Figure 9).    (Figure 8b), where eastward motion of +0.2 to +0.6 cm/year was found, with peaks reaching +1.0 cm/year.

Geological Interpretation
The P-SBAS results highlight the presence of surface deformation processes that are related to known natural and anthropogenic hazards characterizing the study area (see Section 2): land subsidence in areas of gas extraction; landslides and erosion landforms; and subsidence due to settlement of modern infrastructure, intermingled with wider-scale motion patterns. It is worth noting that any long-wavelength deformation signals potentially affecting the study area are likely to have been subtracted from the datasets during the removal of residual orbital ramps (see Section 3.2), given that all long-wavelength signals (either deformation-or orbit-related) are typically treated as orbital errors, due to their similar spatial characteristics (e.g., [87]). This does not represent an issue for investigations of local scale processes, while it could impact studies of processes occurring across large scales (e.g., regional scale movements linked with tectonics), as their signals were suppressed during the above processing step.

Coastal Land Subsidence at the Capo Colonna Promontory
The vertical and horizontal displacement field at Capo Colonna in 2014-2020 unequivocally highlights that the promontory is affected by land subsidence (Figure 7). The spatial configuration (albeit due to different causative factors) is, indeed, not dissimilar to other subsiding areas in Calabria region (e.g., Gioia Tauro [88]). V U increases from inland and from the promontory coastline toward the central agricultural area, which is bounded by the spatial distribution of the gas wells (Figure 7a). The full displacement history for a coherent target at the center of the bowl is shown in Figure 10a through its LOS time series at the nearly 300 dates composing the dataset. Despite the lower density of coherent targets across the rural landscape, the east-west horizontal field makes evident that the motions converge toward this agricultural area, and this is the location around which the subsidence bowl appears to be centered (Figure 7b), i.e., west of the Hera Lacinia 001 production well. The evidence of such an east-west deformation field at the promontory is similar to earlier observations based on COSMO-SkyMed 2011-2014 data [28].
These results highlight a spatial association between land subsidence and gas extraction infrastructure that has been observed in previous studies. In this regard, in chronological order of publication, Table 1 summarizes journal papers, proceedings, and technical reports that include InSAR analyses of land subsidence at Capo Colonna. Of these, the majority focused on InSAR datasets retrieved from one to two SAR sensor datasets. Whereas [46] provided not only the most recent assessment using COSMO-SkyMed ascending and descending (2011)(2012)(2013)(2014)(2015) and Sentinel-1A ascending (2014-2016) data, but also a 25-year-long review of local land subsidence through the combination with other C-and X-band InSAR datasets from 1992 to 2010. Comparison of the Sentinel-1 P-SBAS outputs with past observations available in the literature allowed corroboration of the spatial patterns and the order of magnitude of the deformation rates affecting the study area. In particular, the present results corroborate the spatial patterns of land subsidence that were found with ERS-1/2 ascending and descending data processed with SBAS [38] and Permanent Scatterers InSAR (PS-InSAR) [39,46], ENVISAT ascending data with SBAS [38] and Persistent Scatterers Pairs-DIFferential In-SAR (PSP-DIFSAR) [39,46], TerraSAR-X ascending and descending data with SBAS [45,46], and COSMO-SkyMed ascending and descending data with PSP-DIFSAR and the Coherent Pixel Technique-Temporal Sublook Coherence (CPT-TSC) [46]. It is worth noting that a quantitative comparison of 2014-2020 deformation rates with those observed in these past studies is not viable due to their different time coverage.  Such corroboration is achieved through a significant improvement in terms of both InSAR processing and post-processing of coherent targets and their time series. In particular, most of the previous studies on Capo Colonna relied on LOS displacement values only [38,39,46] and, when the authors of those papers attempted the extraction of deformation components, they limited it to either the vertical [45,46] or east-west [28] deformation field. Therefore, the present study is the first that, at least for the last six years, is able to fully constrain the geometry of the land subsidence process at the promontory not only along the LOS, but also along both the vertical and east-west directions.
Furthermore, a closer look at the temporally dense InSAR time series highlights that the displacement trend of land subsidence has been almost constant over the study period (Figure 10a), with a magnitude of total displacement that makes a clear separation between the subsiding points and the more stable ones (Figure 10b). Pending an in-depth time series analysis (also through statistical methods, which is beyond the scope of this research), the nonlinearities visible throughout the time series do not seem to coincide with any significant deceleration or inversion of the land subsidence trend. Future research will also need to focus on the temporally confined trend variations occurred in July-October 2016 and September-October 2018, which will have to be further investigated through comparison with independent in situ monitoring data (not available for this research).
In this regard, the present analysis allows a revisit of the assessment made based on an earlier Sentinel-1 dataset, in which a slight decrease in the annual LOS velocities of subsidence around the gas extraction sites was suggested [46]. That dataset consisted of 46 ascending mode images only, collected by the Sentinel-1A satellite alone, every 12 days, from October 2014 to May 2016, and was processed with an intermittent SBAS [89] processing method capable of extending the coverage of motion results across the full range of land-cover types by relaxing the approach to selecting image pixels to process [90]. Therefore, despite the high density of coherent pixels over the agricultural area provided in [46], the limited length of the image datasets somehow might have affected the estimation of the deformation velocity by tying it to a shorter time period. On the other hand, with nearly 300 displacement estimates composing the InSAR time series in each Sentinel-1 observation geometry (Figure 11), the current P-SBAS analysis provides a more robust evidence of the uninterrupted land subsidence long-term behavior and its rate of occurrence.

Landslides and Erosion Landforms
South of the city of Crotone and northwest of the Capo Colonna promontory, the V U map highlights two elongated patterns and a more circular one of vertical motions (Figure 9a) that cannot be attributed to land subsidence, given their shape and underlying topography, but instead to mass movement processes. As both ascending and descending mode data were exploited in the analysis, the directional variety of mass movements in this zone was taken into account. In turn, both east-facing and west-facing slopes were visible to the Sentinel-1 LOS and covered with targets from at least one P-SBAS output dataset.
The first pattern of vertical motions, generally up to −0.75 cm/year with a few peaks of −1.0 to −1.4 cm/year, covers the narrow coastline section where, according to the Hydrogeological Setting Plan (Piano Stralcio per l Assetto Idrogeologico, PAI) updated in 2016 by the River Basin Authority of Calabria Region [86], a dormant landslide area with a medium hazard level (i.e., hazard index IP2) is mapped (i.e., polygon A in Figure 9a). Within the boundary of this polygon, the velocity along the steepest slope direction is generally lower than 5 cm/year, with some peaks at 30-50 cm/year. As per the simplified geology reported in Figure 1c [30], the bedrock lithology is made of Pliocene-Early Pleistocene clays, which are known for their susceptibility to local slope-wash processes, particularly in non-vegetated areas (Figure 2; [49]; see Section 2). According to the Copernicus Land Monitoring Service Coastal Zones dataset (Figure 5a; [85]), the land cover at the slope foot is mainly "industrial, commercial, public and military units (other)", "continuous urban fabric", and "construction sites". These classes match with the Google Earth satellite observation of residential buildings, leisure estates, and local commercial activities distributed along the SP49 road connecting Crotone with the Capo Colonna promontory and the southern part of the province.
The second elongated pattern of vertical motions, up to nearly −3.0 cm/year, is located on top of a superficial landslide zone (i.e., polygon B in Figure 9a) classified in PAI 2016 [86] as active with a high hazard level (i.e., hazard index IP3), and where calanchi landforms (Figure 9b) have been recognized [49]. The velocity along the steepest slope direction reaches 10-50 cm/year at several locations within the mapped landslide zone, though in this case this may not be the predominant deformation direction for the whole landslide body, given the existence of some accumulation zones, as well as some scarps. Displacement time series for the coherent targets within this zone show a clear trend reaching up to −15 cm of total LOS displacement in the observation period (Figure 10c). More to the west, a third cluster of coherent targets moving vertically up to −2.4 cm/year defines a nearly circular pattern indicating deformation at another calanchi landform (i.e., polygon C in Figure 9a) [49] that is classified in the PAI 2016 [86] as a dormant superficial landslide zone with a medium hazard level (IP2). Velocity values along the steepest slope direction are lower than 5 cm/year across most of the feature. Both the landslide areas exhibit a rural land cover ranging from sparse vegetation to scrubland and grassland (Figure 5a).
Since the P-SBAS results provide more recent information than the data used for the PAI 2016 classification, the above velocity values suggest the opportunity to update the state of activity of these landslide and calanchi landforms, by applying the PSI-based matrix approach [84,91] that was originally tested in northern and central Calabria. By referring to the displacement rates observed in 2014-2020 at each mass movement area, it turns out that, while the active state of the IP3 landslide body (polygon B in Figure 9a) is confirmed and the intensity of the process belongs to the category of very slow processes (i.e., velocity of 16 mm/year to 1.6 m/year; according to [92]), the hazard classification of the IP2 landslide bodies (polygons A and C in Figure 9a) should be revised from dormant to active with extremely slow (i.e., velocity less than 16 mm/year) and very slow intensity, respectively, if considering their vertical rates (while estimates along the steepest slope direction suggest that both processes occur with very slow intensity). On the other side, given that most of these landslide bodies are calanchi, the proposed reclassification based on the found motions has to be interpreted in the context of the typical behavior of these landforms, which are characterized by washing and weathering processes, and conditioned by local rainfall regime.
Future in-depth analysis of these processes could also account for local precipitation records to investigate any potential correlation with the deformation observed in the satellite results and highlight impacts on landslide triggering and/or reactivation.

Ground Deformation in the City of Crotone
Across the whole P-SBAS datasets, the other zone with negative LOS velocity values that can be clearly detected is located over the harbors of the city of Crotone ( Figure 6). In particular, this pattern comes out very distinctively even with the naked eye in the ascending geometry dataset (Figure 6a), where LOS velocity is estimated at up to −1 cm/year. The rest of the city instead appears overall stable in both geometries, as expected for an historic urban settlement (Figure 6a,b). However, this evidence links with the reference point that was set within the urban area (see green star in Figure 6), and so any deformation signal is relative to such point.
By analyzing the P-SBAS data in more detail, and specifically the vertical displacement field (Figure 8a), coherent targets with V U values reaching −0.75 to −1.3 cm/year are distributed along the 1.7 km-long breakwater of the new harbor, as well as along the smaller breakwater and over the commercial and residential buildings overlooking the old harbor. A reasonable explanation of such subsidence might be the typical settlement of the harbor infrastructure, also as a consequence of recent interventions and maintenance works. In this regard, it is worth recalling that LOS displacement rates between −0.2 and −0.5 cm/year had been already found with ERS-1/2 and ENVISAT data in the period 1995-2009 [38].
The horizontal field ( Figure 8b) provides a further interesting insight, given that both the breakwaters and most of the old town (including the massive Castle of Charles V) fall within a wide pattern of eastward deformation, with V E ranging between +0.2 and +0.6 cm/year. Such pattern corroborates the evidence of +0.3 to +0.5 cm/year found with COSMO-SkyMed data for the period 2011-2014 [28]. Any slight difference between the two InSAR analyses could mainly be due to the different time periods surveyed, as well as the different techniques used to process the SAR data and the respective processing parameters that were set up, above all the reference point and its precise location within a known tectonically active region. Continuous geodetic monitoring data from two GPS stations installed in Crotone (CROT and KROT) indeed confirmed the presence of significant eastward horizontal motions, alongside northward ones [93]. The latter, however, are not captured in the P-SBAS results due to the known relatively poor visibility of north-south horizontal movements to the SAR sensor LOS, and the V N = 0 assumption that was made to estimate the vertical and east-west deformation field (see Section 3.3). In this case, such an assumption prevents any north-south displacement to be estimated, and also introduces an error in the other two components (vertical and east-west), especially at locations where a significant north-south displacement occurs [81]. Detailed date-by-date comparison of the P-SBAS and GPS monitoring data at the two stations would require the setting of a local reference frame for the GPS stations and the calibration of the P-SBAS results to the same reference system, as done in [28]. This would enable the quantification of the P-SBAS datasets accuracy against external data, though such analysis is beyond the scientific goals of this study.

Potential and Challenges of Big Data
The above results and geological interpretation demonstrate how the paradigm of InSAR studies is rapidly changing not only with the advent of SAR big data, but also owing to the accessibility of cloud/grid infrastructure and chains to process such huge datasets.
The nearly 300 images per single observation geometry that were processed to cover the study area of this paper well represent the typical big SAR stack that a satellite constellation such as Sentinel-1 is already able to provide, now that its two twin satellites have been collecting images globally for more than six years, with a site revisit of up to six days. Even if the site revisit was constrained to 12 days only (as it currently is in geographic areas outside Europe; [9]), the constellation would still collect 30 images per year. To put this number into context, this exceeds the typical minimum amount of C-band SAR images per single observation geometry that InSAR practitioners consider as the minimum requirement to carry out a meaningful and reliable PSI processing (i.e., 15-20 scenes; [5]), regardless of the total period of observation. The latter plays a crucial role. Relatively short time periods might be unable to capture the long-term behavior of the deformation processes of interest, and/or be biased by the occurrence of annual fluctuations or other short-term dynamics that could influence the estimation of deformation rates. Therefore, the use of Sentinel-1 imagery to investigate "silent" and subtly evolving geohazards that manifest over long periods and with trends that can be captured over years (such as the slow-moving landslides and coastal subsidence investigated in the province of Crotone), today poses the technical challenge of processing big SAR stacks with adequate computing resources.
Our experiment proves that nowadays, with the GEP it is possible to process a whole time series of 300 Sentinel-1 SAR images at full swath in less than three days, without the necessity for in-house high-performance computers and data-storage capacity (as an approximate figure with which to compare this time performance, the same processing in a standard desktop PC could take over a week or more, depending on its hardware configuration). The latter aspect is frequently overlooked when planning an InSAR investigation, but it should not be. A stack of 300 SLC Sentinel-1 dual-polarization images alone, prior to any pre-processing or InSAR analysis, requires approximately 2.4 terabytes of disk space, not to mention that these input data need to be sourced from a repository, and this step normally requires a certain amount of time given the conditions for data downloading from catalogues (e.g., a maximum of two simultaneous image downloads as currently allowed to users by the Copernicus Open Access Hub). On the other hand, it is reasonable to expect that a more direct connection to the data repository and a faster data retrieval, combined with a greater amount of computing resources made available to run the processing steps, can significantly improve the overall performance with total processing times shorter than the~3 days achieved in this experiment. Other P-SBAS studies carried out over the same and larger regions of interest across Italy [54,59] report processing time demands of less than one day, though they refer to shorter Sentinel-1 stacks (i.e.,~120-180) and a different cloud-computing infrastructure (i.e., AWS). Therefore, a one-to-one quantitative comparison of time demands cannot be made against these studies.
Nonetheless, there is no doubt that a processing run in the GEP saves time that can instead be invested by users to tailor the processing parameters, post-process and interpret the output InSAR products to address their scientific and/or technical goals (e.g., creating a hazard map, estimating landslide intensity, modeling the subsidence process). In this regard, in the GEP, the selection of parameters such as the number and acquisition date of the input SAR data, the study area (clip to AOI), the processing thresholds (e.g., temporal coherence), and the location of the reference point, is allowed through the geo-portal within the framework of a consolidated and standardized P-SBAS processing routine. In the case of Capo Colonna and Crotone, this functionality has definitely been advantageous compared to previous studies that exploited already processed and ready-to-use InSAR datasets [28,39,46]. At equal opportunity of reliance on trustable InSAR datasets produced by means of consolidated processing routines, in the present experiment, we could tailor the processing parameters to suit our purposes (e.g., reference point selection and respective implications, as discussed in Section 5.1.3).
On the other side, compared to those papers investigating land subsidence and regional tectonics in the study area based on InSAR products generated by the authors themselves [38,45,46], the processing times achieved thanks to the use of the GEP were extremely short on a relative scale, especially when considering that the Sentinel-1 stacks used in the present study were one order of magnitude longer (i.e., 300 vs. 40-60 images). A drawback is, of course, due to the current limitations in the range of SAR sensors that can be selected for P-SBAS processing in the GEP. At present, the P-SBAS on-demand service runs with Sentinel-1 data only, and it is also available for ERS-1/2 and ENVISAT SAR data archives via a link with the ESA G-POD platform [55].
Of course, such eased accessibility to reliable InSAR products does not solve the challenge of how to manage the post-processing and geospatial analysis of hundreds of thousands (or even millions) of coherent points. In the present research, given the intended demonstration character of this paper and the regional scale of the geohazard analysis (therefore, with no intention to deepen the geological interpretation beyond the correlation between the observed motions and the type of natural and/or anthropogenic process), the multitude of coherent points was creamed-off by narrowing down the geological interpretation to a relatively small region of 300 km 2 , and to its most relevant displacement patterns. This was done based on the initial LOS velocity values and the extracted vertical and horizontal displacement field, and then through geospatial integration with ancillary data (e.g., lithological map, landslide inventory) and the literature (e.g., evidence from previous InSAR studies). Such an approach is viable and effective when there is strong background knowledge of the geohazards affecting the study area. In other contexts, methods such as hotspot identification and clustering [94], time series analysis and estimation of trend deviations [15,17,95], statistical tests, and probabilistic approaches [16][17][18] are recommended.
Research trends in the field of InSAR are now aiming toward the development of processing chains and approaches to continuously update multi-interferogram processing outputs (as well as their post-processing and automated analysis) every time that one or more new SAR scenes are available for the study area. These approaches may provide key information for monitoring and early warning of the occurrence of anomalous trends and/or acceleration associated with precursors or incipient deformation processes, and thus input alerting and operational monitoring services (e.g., [25,96]). To achieve these goals, efficient approaches for data reprocessing after new data inputs become available are required to optimize processing loads and times (in the order of a few hours; [25]) for the updating of a multi-interferogram processing (e.g., via a sequential estimator; [97]). Currently, the P-SBAS service available in ESA's GEP does not include such functionality, and therefore the benefit of these approaches on its time performance cannot be assessed in this case.

Conclusions
In this study, six full years of Sentinel-1 SAR ascending and descending mode imagery were processed with the P-SBAS service available in ESA's GEP to generate the longest ever interferometric stacks so far used for a single study site and published in the specialist SAR remote sensing literature. The multitude of coherent points retrieved and the respective LOS displacement time series were then analyzed to develop use-cases of geological interpretation related to "silent" and subtle hazards-i.e., slow-moving landslides, erosion landforms, and coastal subsidence-affecting different locations spread across the province of Crotone in Calabria, southern Italy.
The results achieved make this research novel because they demonstrate that Sentinel-1 time series in the order of nearly 300 images (per observation geometry) can be processed in less than three days, owing to the high-performance computing infrastructure available in the GEP. Tens of thousands of coherent points can therefore be retrieved through P-SBAS, with a level of spatial density that is suitable to depict deformation patterns in rural and agricultural land-cover types as well. Furthermore, the very frequent displacement time series provided for each coherent point enable a robust identification of not only the longterm deformation trend of natural and anthropogenic hazard processes, but also of their short-term variations on a weekly basis, which previous analyses relying on shorter data stacks (or long stacks with less frequent observations) for the same study area could struggle to unveil. This comes out evidently in the characterization of the temporal evolution of land subsidence at the Capo Colonna promontory. Any investigations of shorter time spans, indeed, could be biased by the occurrence of seasonal or intra-annual variations, and therefore would not be able to provide a robust estimate of the long-term trend. On the other hand, long time spans with less frequent observations (e.g., on the monthly or bimonthly scale) may not be able to sense shorter scale non-linearities. Therefore, the geo-information provided by the 300 image-long datasets with weekly revisit acts as a key input for investigations of geological processes.
In this regard, the choice of a study area that was previously investigated with InSAR techniques proved effective to highlight the advantages brought by P-SBAS processing of such big SAR data, and thus helped to demonstrate the concept of how to develop applications with such big InSAR products for geohazard assessment and environmental monitoring purposes. Post-processing of P-SBAS results according to well-established methods included the vertical and east-west deformation field estimation, and the PSIbased matrix approach for landslide activity assessment.
In the case of Capo Colonna, the vertical and east-west deformation field estimation allowed constraining of the geometry of land subsidence and locating the center of the subsidence bowl, thus confirming and bringing in new knowledge compared to previous InSAR observations that were limited to analysis of LOS displacements, or only one of the two deformation components. In Crotone, the method helped to discriminate different sectors of the city (harbors, old town, inland settlement) and to detect horizontal motions that relate to the regional tectonics.
Finally, the exercise of updating the existing landslide inventory by means of the PSIbased matrix approach was intentionally meant as a simple, but effective proof-of-concept of the types of value-adding that users can undertake if they have access to the platform and become familiar with generating and exporting P-SBAS products to input them into their workflows for the generation of thematic products and development of downstream applications.

Informed Consent Statement: Not applicable.
Data Availability Statement: Sentinel-1 and Sentinel-2 data are openly available at the Copernicus Open Access Hub. NASA's SRTM data are openly available at the Earth Explorer data portal.