**Monitoring Land Subsidence Using Remote Sensing**

Editors

**Massimo Fabris Nicola Cenni Simone Fiaschi**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Massimo Fabris Department of Civil, Environmental and Architectural Engineering University of Padova Padova Italy

Nicola Cenni Department of Geosciences University of Padova Padova Italy

Simone Fiaschi UCD School of Earth Sciences University College Dublin Dublin Ireland

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Remote Sensing* (ISSN 2072-4292) (available at: https://www.mdpi.com/journal/remotesensing/ special issues/Monitoring Land Subsidence Using RS).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-1388-1 (Hbk) ISBN 978-3-0365-1387-4 (PDF)**

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## *Editorial* **Editorial for Special Issue "Monitoring Land Subsidence Using Remote Sensing"**

**Massimo Fabris, Nicola Cenni and Simone Fiaschi**


#### **1. Introduction**

Land subsidence is a geological hazard that affects several different communities around the world. The main consequences of subsidence can be related to environmental degradation, damage to buildings, and interruption of services [1]. The effects produced by the lowering of the ground level on building and infrastructure can be considered as a major problem in many countries. More than 150 cities all over the world have endured land subsidence with rates up to tens of centimeters per year ([2] and references therein).

Land subsidence can have both natural and anthropogenic origin: natural subsidence can be due to the compaction of lithological layers of the soil, the oxidation of peat, and geodynamic processes (e.g., tectonic-plate movements, volcanism [3]); anthropogenic subsidence derives mainly from the compaction of aquifers associated with groundwater/oil/natural gas extractions, drainage of organic soils, underground mining, hydrocompaction, sinkholes, stress provided by newly-built man-made structures, and thawing permafrost ([4,5] and references therein); the combination and coexistence of these factors have a strong negative impact on the territory [1]. The effects of this global problem are more evident along transitional environments, such as coastal areas, deltas, wetlands, and lagoons, which are becoming increasingly vulnerable to flooding, storm surges, salinization, and permanent inundation [6]. In these areas, the effects of subsidence are linked also to the retreat of coastlines and disappearance of emerged surfaces [7,8].

The ground surface movements due to anthropogenic activities have been deeply investigated, particularly in critical environments such as high-urbanized areas and coastal zones [9].

The monitoring activities allow to acquire useful information that can be used to prevent damage to buildings and infrastructures, plan more sustainable urban development, and mitigate the risk; the knowledge of the temporal and spatial distribution of the ground surface deformations is essential to delineate the areas most affected by subsidence and to understand the involved mechanisms [1].

In the past, only the traditional leveling technique was used for the monitoring of land subsidence; however, despite its high accuracy, its use has been reducing with time due to the high costs and limited number of points potentially measurable.

The monitoring of ground movements made great progress in the last decades with the development of Global Positioning System (GPS)–Global Navigation Satellite System (GNSS) [10] and Interferometric Synthetic Aperture Radar (InSAR) [11] technologies together with the use of different approaches, from analytical to 3D numerical, for the analysis of the involved physical processes [9]. Starting from the 2000s, further advances in the satellite-borne and in-situ technologies made the monitoring of Earth surface motions an easier and more common geodetic task [12].

**Citation:** Fabris, M.; Cenni, N.; Fiaschi, S. Editorial for Special Issue "Monitoring Land Subsidence Using Remote Sensing". *Remote Sens.* **2021**, *13*, 1771. https://doi.org/10.3390/ rs13091771

Received: 22 April 2021 Accepted: 28 April 2021 Published: 1 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The GPS/GNSS technique allows the 3D monitoring of points based on continuous or repeated surveys with lower costs in comparison with the leveling approach, but, generally, with low spatial resolution.

The Differential Interferometric Synthetic Aperture Radar (DInSAR) approach, which was developed from the InSAR technique, has the advantages of wide spatial coverage, high spatial resolution, all-weather working conditions, and cost effectiveness [13].

However, the standard DInSAR technique is strongly limited by the presence of atmospheric effects that reduce the accuracy, and by the lack of time-series [1,2]. These limitations were overcome by the development of multi-temporal DInSAR techniques such as the Permanent Scatterer Interferometry (PSI, which looks for point-like scatterers) and the Small Baseline Subset (SBAS, which looks for distributed scatterers within the cell resolution) that allow improved time-series analysis with millimetric accuracies [2,14]. Both methods are based on the analysis of stacks of SAR acquisitions from which are extracted suitable measure points. The PSI and SBAS approaches can be successfully applied over many different ground targets such as buildings, infrastructures, outcrops, base soils, low-vegetated area, etc. [11], and in several research fields, such as tectonics and volcanology [15], landslides [10,11,16], clay deposits deformations [5,17], and groundwater/oil/natural gas depletion [9,18].

The possibility to combine the results obtained from the analysis of SAR images acquired in both ascending and descending geometries allows the computation of the vertical and east–west components of movement, which is fundamental for data interpretation and for the integration with other techniques [18].

Despite the large advantages of the multi-temporal DInSAR methodology, the acquisition of ground-based data, mainly from leveling and GNSS techniques, is crucial for validation and integration of the satellite measurements. Nowadays, the land subsidence monitoring mostly relies on the integration between available and reliable InSAR data and sparse, but highly accurate, GNSS measurements, which are often compared to other geodetic observation methods [9–12,15,19,20].

This Special Issue consists of nine individual works that used different approaches and methodology for land subsidence applications. In the next section, each work is presented and the general contribution to this Special Issue is summarized.

#### **2. Overview of Contributions**

The papers included in this Special Issue cover a wide spectrum of applications related to the land subsidence monitoring in different areas of the world using remote sensing techniques integrated with ground-based data. The accepted works are presented here in order of publishing.

Zhou et al. [2] analyzed the overexploitation of groundwater in the Beijing–Tianjin– Hebei (BTH) area in China, in the 2012–2018 period, using the Interferometric Point Target Analysis (IPTA) with Small Baseline Subset (SBAS) technique. Authors used 126 RADARSAT-2 and 184 Sentinel-1 images acquired in ascending and descending orbits, respectively, to derive land subsidence rates in the study area: they validated the results using 72 leveling benchmarks and compared the measured vertical land motion with data of groundwater monitoring wells from 2012 to 2015. Moreover, the authors correlate the detected land subsidence rates with eleven subsidence features and land use types: the results showed serious vertical land motion with rates up to 131 mm/y. They demonstrate that the land subsidence changes are consistent with the seasonal trends of the groundwater level changes.

Gido et al. [5] used the Permanent Scatterer Interferometry (PSI) technique and historical leveling data to study the ground surface deformation of Gävle city in Sweden: two ascending and descending Sentinel-1A/B datasets (91 images in total), collected between January 2015 and May 2020, were processed and analyzed together with a long record of a leveling dataset (4 leveling lines), covering the period from 1974 to 2019. The authors performed the comparison between the obtained data at some locations showing a close

agreement between the subsidence rates extracted from precise leveling and PSI. Land subsidence rates were compared also with the geological information of the analyzed area: they suggest that the land subsidence (with maximum displacement rate that reaches up to −6 ± 0.46 mm/y in the LOS direction and only in localized deformation zones) occurred due to relatively weak subsurface layers (hazard zone reported as an artificial fill area) that either was affected by loading of new constructions or by hydro-compaction.

Doke et al. [15] performed an InSAR time series analysis of the Hakone Volcano (Japan) from 2006 to 2011 using the SBAS method and ALOS-PALSAR scenes (24 from the ascending and 22 from the descending orbits, respectively). The authors corrected the obtained InSAR displacements using the available GNSS data. The results showed highly localized subsidence (500 m in diameter, with maximum vertical rates of about 25 mm/y) to the west of Owakudani from 2006–2011. The authors suggested that the land subsidence was caused by a reservoir contraction at approximately 700 m above sea level. Based on the structure of the hot spring wells and their chemical components, they suggested that the contraction source can be considered as a reservoir containing hydrothermal fluids, which demonstrates the feasibility of the InSAR technique to monitor hydrothermal activity in shallow parts of volcanic areas.

Even et al. [18] applied InSAR using both PS and SBAS approaches to study the complex displacement field caused by convergence and operational pressure changes of the natural gas storage field at Epe (N–W Germany). The authors processed 86 and 118 SAR images, respectively, in ascending and descending orbits; they compared the InSAR results with leveling data acquired during three surveys between 2015 and 2017 (517 points in total) and ground water measurements at 97 locations. The authors combined separately the different components of the phase model (geometric orbit combination) for a better understanding of the phenomena that contribute to the displacement field; in addition, a method that allows to perform an orbit combination based on simplistic geomechanical modeling of the spatial displacement field was presented (Multi-Mogi approach). They demonstrated that the InSAR-derived displacements were in reasonable agreement with the leveling data taking into account the geometric orbit combination, and in good agreement with the Multi-Mogi approach; for the vertical components, the comparison with leveling data provided Root Mean Square (RMS) of 3.41 mm/y and 2.39 mm/y for the geometric orbit combination and for the Multi-Mogi approach, respectively.

Benetatos et al. [9] presented a multi-physics investigation of the ground movements related to the cyclical and seasonal injection and withdrawal of natural gas in/from a depleted reservoir located in the Po Plain area (Italy) using the Persistent Scatterer Pairs InSAR (PSP-InSAR) approach and GNSS data. The authors developed an integrated geological, fluid-flow, and geomechanical numerical modeling approach to reproduce the main geometrical and structural features of the involved formations. They processed 432 and 428 SAR scenes (from 2003) in both ascending and descending orbits, respectively, from RADARSAT-1/2 and Sentinel-1 satellites, and daily data from a continuous GNSS (CGNSS) station (from 2008), using the Network approach and the Precise Point Positioning data processing. They found (i) agreement between the InSAR and the GNSS results; (ii) gentle long-term subsidence trends; (iii) a strong correlation between the cumulative volumes curve of the gas storage and the historical series of the ground displacements above the reservoir, considering both the vertical and east–west planimetric components; and (iv) cyclical subsidence/uplift limited to the field area.

Grgi´c et al. [12] showed the conjoint analysis of vertical land motion of the Dubrovnik area (Croatia) using 75 ascending Sentinel-1 images from 2014 to 2020, continuous GNSS observations at Dubrovnik site obtained starting from 2000, differences of the sea-level change derived from all available satellite altimeter missions for the study area and tide gauge measurements in Dubrovnik starting from 1992. The data from the CGNSS station were used to correct the obtained InSAR ground motion rates and to reference the motions with respect to an absolute reference frame. They compared and analyzed trends obtained with the different techniques in the overlapping period, from 2014 to 2020: the results showed vertical land motion velocities in the order of some mm/y with only some limited areas characterized by rates exceeding −5 mm/y.

Sopata et al. [3] focused on describing vertical surface displacements related to seven mining-induced tremors in the Upper Silesian Coal Basin in the south of Poland using the standard DInSAR approach. The authors processed 15 interferograms of Sentinel-1A/B satellites from March to May 2017 in ascending orbit, which overlap with seven reported tremors of the rock mass of magnitude ML = 2.3–2.6. As a result that the obtained land subsidence isolines showed residual signal noise, the authors developed a procedure to eliminate the occurring irregularities interpolating the subsidence profile by means of eighth-degree orthogonal polynomials and manually entering the corrections to the surface distribution of isolines for each individual subsidence. Moreover, they evaluated the threats to building structures according to the classification used in mining areas. They found that the structures with resistance lower than the limit values of land subsidence speed can be exposed to higher damage risk.

Mancini et al. [21] developed a PSI-based workflow to process dual-orbit SAR observations with open-source tools complemented by the use of GNSS observations as constraints for the global reference frame and final accuracy assessment of the vertical and horizontal velocity maps. The authors investigated the land subsidence processes in the eastern sector of the Po Plain (Italy) and in particular in the metropolitan area of the Bologna city analyzing interferometric and GNSS observations acquired between 2015 and 2019. With respect to the InSAR analysis, they added a procedure to refer the LOS-projected velocities to the GNSS reference frame and an algorithm for decomposition analysis. The validation of the velocity maps through the comparison between the decomposed InSAR and GNSS annual rates provided differences at the millimeter level, which confirmed the substantial agreement between the PSI and GNSS measurements.

Cenni et al. [22] highlighted the spatial and temporal evolution of the land subsidence in the Po River Delta (PRD) area (Italy) analyzing time-series obtained from CGNSS stations using a moving window approach temporally overlapping with both the surveys of a new GNSS network (PODELNET—PO DELta NETwork) (46 non-permanent sites measured in 2016 and 2018), and InSAR data (SBAS processing of Sentinel-1A/B images from 2014 to 2017). The authors investigated the integration between these data: since the InSAR technology does not perform well in high vegetated areas or areas with high temporal variability, the PODELNET sites represent an important improvement for the monitoring of the land subsidence in the PRD. Moreover, they highlighted that an integrated monitoring system, combining GNSS and InSAR data, permit to overcome the limits of the two techniques. The obtained results suggested that the land subsidence velocities in the easternmost part of the area of interest (of about −10 mm/y) were characterized by values greater than the ones located in the western sectors (in the order of −5 mm/y), which can be linked to the different morphological characteristics of the subsoil in the PRD.

**Author Contributions:** The three authors contributed equally to all aspects of this Editorial. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The guest Editors would like to sincerely thank all the Authors who contributed to this Special Issue and the Reviewers who dedicated their time and provided the Authors with valuable and constructive recommendations. They would also like to thank the Editorial team of Remote Sensing for their support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Land Subsidence Response to Di**ff**erent Land Use Types and Water Resource Utilization in Beijing-Tianjin-Hebei, China**

#### **Chaofan Zhou, Huili Gong, Beibei Chen, Mingliang Gao, Qun Cao, Jin Cao, Li Duan, Junjie Zuo and Min Shi**


Received: 24 December 2019; Accepted: 28 January 2020; Published: 1 February 2020

**Abstract:** The long-term overexploitation of groundwater leads to serious land subsidence and threatens the safety of Beijing-Tianjin-Hebei (BTH). In this paper, an interferometric point target analysis (IPTA) with small baseline subset InSAR (SBAS-InSAR) technique was used to derive the land subsidence in a typical BTH area from 2012 to 2018 with 126 Radarsat-2 and 184 Sentinel-1 images. The analysis reveals that the average subsidence rate reached 118 mm/year from 2012 to 2018. Eleven subsidence features were identified: Shangzhuang, Beijing Airport, Jinzhan and Heizhuanghu in Beijing, Guangyang and Shengfang in Langfang, Wangqingtuo in Tianjin, Dongguang in Cangzhou, Jingxian and Zaoqiang in Hengshui and Julu in Xingtai. Comparing the different types of land use in subsidence feature areas, the results show that when the land-use type is relatively more complex and superimposed with residential, industrial and agricultural land, the land subsidence is relatively more significant. Moreover, the land subsidence development patterns are different in the BTH areas because of the different methods adopted for their water resource development and utilization, with an imbalance in their economic development levels. Finally, we found that the subsidence changes are consistent with groundwater level changes and there is a lag period between land subsidence and groundwater level changes of approximately two months in Beijing Airport, Jinzhan, Jingxian and Zaoqiang, of three months in Shangzhuang, Heizhuanghu, Guangyang, Wangqingtuo and Dongguang and of four months in Shengfang.

**Keywords:** land subsidence; IPTA; land-use; water resource utilization; groundwater change

#### **1. Introduction**

Land subsidence, with decreasing land elevation, has become an important environmental geological phenomenon, and a common urban geological disaster in many large cities. More than 150 cities all over the world have endured land subsidence with rates up to tens of centimeters per year. Affected areas include China [1–3], Canada [4,5], Mexico [6–8], Italy [9–13], the United States [14], and others. In China, approximately 100 cities have undergone land subsidence to different degrees, and the total area of subsidence with cumulative subsidence exceeding 200 mm has reached 79,000 km<sup>2</sup> [15–17]. Land subsidence in North China Plain (NCP), where Beijing-Tianjin-Hebei (BTH) is located, is most serious in China [18]. Land subsidence in NCP has become one of the most widely distributed geological hazards, especially in BTH [19]. The area of land subsidence in BTH accounts for approximately 90%, and the maximum cumulative subsidence accounts for 1.7 m, 3.4 m, and 2.6 m, respectively [20].

Compared with traditional earth observation methods, such as global positioning system (GPS), leveling, and layer-wise mark measurement, Interferometric Synthetic Aperture radar (InSAR) is a new earth observation technique used to obtain ground deformation information [21,22]. Synthetic Aperture Radar differential interferometry (DInSAR), which was developed from the InSAR technique, has the advantages of wide coverage, working on all-weather conditions, and high-precision monitoring of ground deformation. However, the DInSAR technique is limited by spatial and temporal decorrelation, atmospheric errors, orbit errors, and terrain errors. The Persistent scatterer InSAR (PS-InSAR) and Small Baseline Subset InSAR (SBAS-InSAR) techniques can overcome the limitations of DInSAR to some degree by analyzing the pixels that retain stable scattering characteristics in a time series to obtain ground deformation information [22–27]. Moreover, the SBAS-InSAR technique can reduce the decorrelation influence by recognizing good coherence pixels without strong backscatter and then obtaining time series deformation [16,28].

Due to the continuous development of the cities, BTH urban agglomeration has formed the largest groundwater funnel in the world. Approximately 70% of the water per year used in the region comes from groundwater [29,30]. Nearly 20 billion m3 of groundwater is exploited per year, which forms the largest groundwater depression cone in the world, with an area of more than 70,000 km<sup>2</sup> [31,32]. The overexploitation of groundwater leads to the development of serious regional land subsidence [33–36]. By June 2014, the area of the BTH region with annual subsidence exceeding 50 mm reached 4569 km2. Meanwhile, land use change affects groundwater balance to a certain extent and leads to the evolution of land subsidence [37]. However, previous studies on land subsidence in the BTH region have mostly focused on the acquisition of subsidence information and the analysis of subsidence characteristics in small areas. Studies considering land subsidence over a large area and its response to different land use types and water resource utilization are relatively rare.

In this paper, the Interferometric Point Target Analysis (IPTA) with SBAS-InSAR technique was used to derive the land subsidence in the BTH region from 2012 to 2018 with 126 Radarsat-2 and 184 Sentinel-1 images. Then, the development of land subsidence in BTH was preliminarily investigated. On this basis, we analyzed the spatial and temporal evolution characteristics of typical subsidence features (Section 4). Finally, we discussed the land subsidence response to different land use types (Section 5.1) and water resource utilization (Section 5.2) in typical subsidence features.

#### **2. Study Area and Data Source**

#### *2.1. Study Area*

The BTH includes Beijing, Tianjin, and the typical region of Hebei. The land area is approximately 2.18\*10<sup>5</sup> km2, with a total population of approximately 110 million people. The automobile industry, electronics industry, machinery industry, and metallurgical industry are the main industries in this region, forming the main high-tech and heavy industrial base of the country. BTH is located in the northern North China Plain, which is the most extensive land subsidence area in China.

Land subsidence in the Beijing Plain first appeared in 1935; before 1952, the maximum land subsidence was only 58 mm. However, since the 1950s, land subsidence developed rapidly because of the development of industry and agriculture in the Beijing Plain area caused by the increase of groundwater exploitation and the significant decrease in groundwater level [38,39]. By 2015, the cumulative subsidence had reached 1.4 m and the area of land subsidence exceeding 100 mm had

reached 4000 km2. The areas affected by land subsidence comprise more than 60% of the Beijing Plain and form six land subsidence funnels [40–44].

In 1932, land subsidence first occurred in the Tianjin area. The development of land subsidence in Tianjin is divided into four main stages [45]. The first and second stages occurred from 1923 to 1957 and 1958 to 1966, respectively. Due to the low intensity of groundwater exploitation, the subsidence reached only 7.1 to 12 mm/year and 30 to 40 mm/year, respectively. However, from 1966 to 1985, the subsidence reached 80 to 100 mm/year because the amount of groundwater exploitation increased. After 1985, because the reduction in groundwater exploitation, land subsidence in the urban area was controlled at approximately 15 mm/year. Through 2015, the average subsidence was 26 mm, and the maximum subsidence was 117 mm [46]. No obvious land subsidence occurred in most parts of the north, but the land subsidence situation in most parts of the south is still grim.

Land subsidence in Hebei appeared later than that in Beijing and Tianjin. From 1950 to 1975, the land subsidence rate was less than 10 mm/year. From 1975 to 1985, land subsidence continuously increased due to the exploitation of deep groundwater in the central and eastern plains, reaching 18 to 104 mm/year [47]. Since 1986, land subsidence has developed rapidly, except for relief in Cangzhou [48]. Land subsidence occurred mainly in the middle and eastern Hebei Plain area. By 2009, there were 14 land subsidence centers, including Cangzhou, Renqiu, Suning, Hejian, Xianxian, Dongguang, Hengshui, Nangong, Pingxiang, Fengnan, Tanghai, Langfang, Baoding and Handan. By the end of 2009, the cumulative land subsidence of the Cangzhou subsidence center had reached 2580 mm, and the area with cumulative land subsidence greater than 1000 to 2000 mm had exceeded 11.3 km2 [49].

#### *2.2. Data Source*

In this paper, we applied 126 Radarsat-2 and 184 Sentinel-1 images to derive land subsidence information in the BTH area. Radarsat-2 is an Earth observation satellite which was successfully launched on 14 December 2007, with a revisiting time of 24 days. Radarsat-2 operates in the C-band with the HH polarization mode and a wavelength of 5.6 cm. Sentinel-1 as the first satellite developed by the European Commission (EC) and the European Space Agency (ESA) for the Copernicus global earth observation project which launched in April 2014. This satellite carries a C-band SAR, with a revisit period of 12 days. For this paper, we chose an Interferometric Wide swath (IW) imaging mode with a medium resolution (5 m × 20 m). The location and information of the Radarsat-2 and Sentinel-1 images are shown in Figure 1a and Table 1. The track number of Radarsat-2 and Sentinel-1 is 26 for descending track and 142 for the ascending track, respectively. The external DEM data we selected are the Shuttle Radar Topography Mission DEM (SRTM DEM) data with 90 m resolution.

**Table 1.** Information of the Radarsat-2 and Sentinel-1 images. R2 and S1 are the abbreviations of Radarsat-2 and Sentinel-1, respectively. The numbers in S1 are the frame numbers of the master image. "Numbers" indicates the numbers of SAR images.


Furthermore, we chose the groundwater monitoring wells near the subsidence features, mentioned in Section 5.2 to compare subsidence changes with groundwater level changes. Because we only had permission to obtain the groundwater level before 2015, only the subsidence obtained by Radarsat-2 from 2012 to 2015 is included in this analysis. The information of groundwater monitoring wells is shown in Table 2.

**Figure 1.** (**a**) The location of typical Beijing-Tianjin-Hebei (BTH) area and the coverage of synthetic aperture radar (SAR) images. (**b**) The location of the North China Plain (NCP). The dark blue line is the middle route of the South-to-North Water Diversion Project (MSWDP). (**c**) The location of BTH in China.


**Table 2.** Information of the groundwater monitoring wells; the monitoring wells are named with the subsidence features. The location of groundwater monitoring wells is shown in Figure 1a.

#### **3. Methods**

#### *3.1. IPTA*

IPTA is a method developed by the GAMMA Company to acquire time series deformation from point targets that have stable spectral characteristics or high backscattering characteristics for single-look complex (SLC) data. Although the IPTA method can accurately obtain the results of nonlinear deformation information in a small area, its ability is often limited when processing large-scale data. However, the SBAS-InSAR technique can reduce the effects of spatio–temporal decorrelation by generating interferograms pairs based on a spatio–temporal baseline and the Doppler centroid frequency shift threshold [22]. Considering the respective advantages of the SBAS-InSAR and IPTA methods, this paper adds SBAS-InSAR into the IPTA method to process SAR images. The specifics steps of the algorithm are as follows:

(1) Master image selection

The master image is selected on the basis of the temporal and spatial baseline, and the Doppler centroid frequency, and its Doppler center is as close as possible to the average Doppler center of the other SAR images. In this paper, the master images for the three Radarsat-2 data and Sentinel-1 data were taken on 23 April 2014 and 25 June 2017, respectively.

(2) SAR images registration

All SAR images are registered using the same radar coordinates. N interferograms are formed based on N + 1 scenes of SAR data in the SLC format, where one is chosen as a master image, and the other images are registered to the master image. This registration includes two steps, which include rough registration based on an image parameter file and precise registration based on the least square method.

(3) Small baseline differential interferograms generation

We selected interferometric image pairs with spatial and temporal baseline shorter than 300 m and 300 days, respectively. The final interferometric image pairs are showed in Figure 2. The final interference pairs for R2-a, R2-b and R2-c are 205, 284 and 233, respectively, and the final interference pairs for S1-116, S1-121 and S1-126 are 177, 174, 183, respectively. The external digital elevation model (DEM) data of SRTM DEM with a 90 m resolution is used to remove the topographic phase signature. The differential interferometric phase is as follows:

$$
\varphi\_{\rm dint} = \varphi\_{\rm def} + \varphi\_{\rm topo} + \varphi\_{\rm atm} + \varphi\_{\rm orbit} + \varphi\_{\rm noise} \tag{1}
$$

where ϕ*dint* is the differential interferometric phase, ϕ*de f* is the deformation phase including the line deformation and non-line deformation phases, ϕ*topo* is the topographic phase, ϕ*atm* is the atmospheric phase, ϕ*orbit* is the orbital phase, and ϕ*noise* is the noise phase.

(4) Interferometric point target extraction

Coherent point target extraction is the foundation and premise of interferometric point analysis technology. The point target extraction method we used is the coherence coefficient method which refers to setting standards to extract points with coherence information [50]. The coherence coefficient was set to 0.3 in this paper.

(5) Interferometric point analysis

Based on the extracted interferometric point targets, we performed a regression analysis to acquire the linear deformation and residual phases. The residual phase includes the atmospheric, non-linear deformation, and noise phases [51]. The non-linear phase shows a high correlation in the spatial domain, while the atmospheric delay shows a high correlation in the spatial domain and low correlation in the temporal domain. Noise shows a low correlation in both the spatial and temporal domains. First, the low-pass filtering in the spatial domain is used to suppress the noise phase from the residual. Secondly, in the temporal domain, the low-pass filtering is again used for spatially filtered residual phases to identify the non-linear phase. Thirdly, we superpose the non-linear phase to the linear

phase to acquire the complete deformation phase. Finally, phase unwrapping is used to acquire the deformation with point targets.

**Figure 2.** (**a**–**c**) is the interferograms for Radarsat-2a, Radarsat-2b, Radarsat-2c data and (**d**–**f**) is the interferograms for Sentinel-1a, Sentinel-1b and Sentinel-1c data, respectively.

#### *3.2. Merging the InSAR Monitoring Results*

(1) Merging the InSAR results in space

In this study, we only consider the difference of deformation in point targets, caused by a different selections of reference points [41]. We first acquired the deformation for each SAR image with the IPTA technique. Secondly, we extracted the same interferometric point targets in overlapping regions. Thirdly, the adjacent strip data results are adjusted to splice different strip data results. The calculation method for the difference is shown in Equation (2).

$$\overline{a} = \frac{\sum\_{i=1}^{n} (a\_{1i} - a\_{2i})}{n} \tag{2}$$

where *n* represents the number of point targets in the overlapping area, and *a* represents the average difference between two strip data results. For the Radarsat-2 data, first, we let *a*1*<sup>i</sup>* be the subsidence at *i* point targets of Radarsat-2a, and let *a*2*<sup>i</sup>* be the subsidence at the *i* point targets of Radarsat-2b. Then, we use Equation (2) to merge the InSAR results of Radarsat-2a with those of Radarsat-2b. Next, let *a*1*<sup>i</sup>* be the subsidence at *i* point targets of the merged result of Radarsat-2a and Radarsat-2b and let *a*2*<sup>i</sup>* be the subsidence at the *i* point targets of Radarsat-2c. Then, we use the Equation (2) to merge the InSAR

results of Radarsat-2a and Radarsat-2b with those of Radarsat-2c. Afterwards, we merge the InSAR results of S1-116, S1-121, and S1-126 using the same steps.

(2) Merging the InSAR results in a time-series

First, the InSAR results derived by the Radarsat-2 and Sentinel-1 images are transformed from line of sight (LOS) to the vertical direction by Equation (3) to ensure the InSAR results have a common direction:

$$d\_{\rm ul} = \frac{d\_{\rm LOS}}{\cos \theta} \tag{3}$$

where θ is the central incident angle, and *dLOS* is the deformation in the LOS direction. Secondly, the InSAR results of Rasarsat-2 were resampled based on the InSAR results of Sentinnel-1. Thirdly, the mean difference between the InSAR results of Rasarsat-2 and Sentinel-1 is carried out in overlapping observation periods from January 2016 to October 2016. Then, the corresponding time-series InSAR result for Sentinnel-1 is resampled based on the InSAR result of Rasarsat-2 with the previously-calculated difference. Finally, the long time-series subsidence from 2012 to 2018 in BTH is derived.

#### *3.3. Impulse Responses Function*

The impulse response function (IRF) method is the corresponding innovation accounting method in the vector autoregressive (VAR) model which was proposed by Sims in 1980 [52]. The IRF is used to measure the impact of a standard deviation shock from the random disturbance term (innovation) on the current and future values of variables. This model can intuitively describe the dynamic interaction between variables and their effects [53,54]. In this paper, the IRF method was employed to estimate the dynamic effect of groundwater level changes on land subsidence changes. The IRF curve is drawn to analyze the impact of groundwater level changes to land subsidence changes. The period corresponding to the peak value of the curve is the lag time of land subsidence. It is necessary to test the stability of the time series data before the IRF because this method is used for stationary data. In our case, an augmented Dickey-Fuller (ADF) test method was used for test, and the confidence coefficient was 5%. The test results are shown in Table 3. For the non-stationary time series values in the original time series data, a first-order difference is carried out, and the ADF test is carried out again for the sequence after this difference to ensure that the sequences used for the IRF analysis are all stationary sequences. The p-value is used to test a statistical hypothesis. The D(p-value) refers to the p-value of an ADF test based on the difference of the time series data. In Table 3, when the p-value or D(p-value) is less than 0.05, the time series data have passed the ADF test.


**Table 3.** The ADF test results of subsidence and groundwater level change in subsidence features from (a) to (j). The sub and Gl in table are subsidence and groundwater level, respectively.

#### **4. Results**

Figure 3 shows the distribution of land subsidence in a part of the BTH plain, which includes Radarsat-2 and Sentinel-1 data overlapped with the average land subsidence rate from 2012 to 2018. The land subsidence in BTH is very serious, and eleven subsidence features have been identified in Beijing, Langfang, Tianjin, Baoding, Hengshui, and Xingtai. We will investigate these subsidence features in the following sections.

**Figure 3.** The average land subsidence rate from 2012 to 2018. The three white dotted frames are the main subsidence area in BTH, which we will discuss in detail in Figures 4–10.

#### *4.1. Land Subsidence in Beijing*

In Figure 4, we find that the spatial distribution of the land subsidence rates in Beijing is quite different. The annual subsidence increased from 2012 to 2013, during which the maximum subsidence increased from 120 mm to 148 mm. The maximum subsidence decreased from 148 mm in 2016 to 106 mm in 2018. Furthermore, the area of land subsidence exceeding 50 mm (Figure 5a) is calculated and four subsidence features (Figure 5b) are identified. From 2012 to 2014, the area of land subsidence over 50 mm shows an increasing trend from 308 km2 to 374 km2 even though the maximum subsidence decreased from 148 in 2013 to 131 in 2014. After 2014, the area of land subsidence over 50 mm and the annual subsidence especially show a decreasing trend. As seen in Figure 5b, the subsidence in the Jinzhan and Heizhuanghu area shows a similar trend, which increased from 2012 to 2013 and decreased from 2013 to 2018. Significantly, this subsidence was maintained at approximately 40 mm from 2012 to 2018 in the Beijing Airport area, so we should pay attention to its development.

**Figure 4.** The annual land subsidence in the typical Beijing region. The subsidence in (**a**–**d**) is based on results from the Radarsat-2 data and (**e**–**g**) show the results from the Sentinel-1 data, respectively. The four red circles are the locations of the subsidence features.

**Figure 5.** (**a**) The time-series subsidence changes and area of subsidence exceeding 50 mm from 2012 to 2018 in the typical Beijing region. (**b**) The annual land subsidence change from 2012 to 2018 for the four subsidence features.

**Figure 6.** The annual land subsidence in typical Langfang and Tianjin regions. The subsidence observed in (**a**–**d**) is based on results from the Radarsat-2 data, while (**e**–**g**) show the results from the Sentinel-1 data. The three red circles are the locations of the subsidence features.

**Figure 7.** (**a**) The time-series subsidence change and area of subsidence greater than 50 mm from 2012 to 2018 in the typical Langfang and Tianjin region. (**b**) The annual land subsidence changes from 2012 to 2018 for three subsidence features.

**Figure 8.** The annual land subsidence in typical Baoding, Hengshui and Xingtai regions. The subsidence in (**a**–**d**) is based on results from the Radarsat-2 data and (**e**–**g**) show results from the Sentinel-1 data. The four red circles are the locations of the subsidence features.

**Figure 9.** (**a**) Time-series subsidence change and area of subsidence greater than 50 mm from 2012 to 2018 in the typical Baoding, Hengshui, Cangzhou and Xinngtai regions. (**b**) The annual land subsidence changes from 2012 to 2018 for four subsidence features.

**Figure 10.** Comparison between InSAR-derived land subsidence and leveling-derived observations at 70 leveling benchmarks, whose locations are shown in Figure 1a with levelings-II. There are 35 leveling benchmarks from 2012 to 2013 in (**a**) and 2017 to 2018 in (**b**), respectively.

#### *4.2. Land Subsidence in Langfang, Tianjin*

In Tianjin and Langfang, the maximum subsidence is relatively severe, with 145 mm compared to the subsidence in the Beijing area. In this typical region, three subsidence features were identified: Guangyang and Shengfang in Langfang and Wanqingtuo in Tianjin. As seen in Figure 6, we found that the land subsidence trend in Guangyang significantly decreases after 2015. However, the land subsidence in Shengfang and Wangqingtuo shows a decreasing trend, although that trend is not obvious. Furthermore, the time-series characteristics of land subsidence and the area exceeding 50 mm in this region are investigated. For 2012 to 2018, the area of land subsidence greater than 50 mm shows a decreasing trend from 344 km<sup>2</sup> to 147 km2. The area of land subsidence greater than 50 mm dropped by almost half from 2015 to 2018, despite a rebound in 2017. From the three subsidence features (Figure 7b), we found that in Guangyang, land subsidence decreased obviously from 60 mm in 2012 to 23 mm in 2018, and the distribution of land subsidence also decreased (Figure 6). However, in Shengfang, the land subsidence increased from 101 mm in 2012 to 105 mm in 2013 and decreased from 105 in 2013 to 77 in 2018. Nevertheless, although the land subsidence decreased, but the subsidence range in space did not decrease (Figure 6). In Wangqingtuo, the land subsidence was largely maintained at around 70 mm, except for 2014 (86 mm). Further, there was basically no change in the land subsidence range in space.

#### *4.3. Land Subsidence in Heibei (Baoding, Hengshui, Xingtai)*

The last subsidence bowl and four subsidence features occur in Baoding, Hengshui, Canzhou, and Xingtai of Heibei province. The maximum subsidence is 131 mm in this area, which is relatively lower compared to the other two areas (Sections 4.1 and 4.2). From Figure 8, we found that the subsidence in this area slowed from 2012 to 2018. Meanwhile, we figure out the area of land subsidence greater than 50 mm, as shown in Figure 9a. The result is similar to that for Beijing; the area of land subsidence greater than 50 mm decreased from 2012 to 2015 and from 2017 to 2018, increased from 2016 to 2017, with the same trend of subsidence. As seen in Figure 9b, we found that the land subsidence in Dongguang, Jingxian, Zaoqiang and Julu experienced an overall decreasing trend from 2012 to 2018, especially after 2015.

#### *4.4. Accuracy Assessment of InSAR Results vs Leveling Data*

To evaluate the accuracy between the InSAR-derived subsidence from the Radarsat-2 and Sentinel-1 data and the leveling data, we converted the InSAR results into vertical displacements using a trigonometric equation and neglected the horizontal component of movement. In this paper, there are total of 72 leveling benchmarks: 35 from 2012 to 2013 to verify the subsidence derived from

the Radarsat-2 images, 35 from 2017 to 2018 to verify the subsidence derived from the Sentinel-1 data, and two from 2012 to 2018 to verify the time series subsidence. The locations of these leveling benchmarks can be found in Figure 1. The precise subsidence was measured at leveling benchmarks from 2012 to 2013 and 2017 to 2018, and we selected the point targets within 200 m of the leveling benchmarks. The results of this comparison are shown in Figure 10. The differences between the two measuring methods vary from 1 mm to 23 mm and from 1 to 20 mm; the standard deviation is 6 mm and 5 mm, respectively.

The results of the two leveling benchmarks from 2012 to 2018 were used to further evaluate the accuracy of InSAR-derived time series subsidence and leveling results (Figure 11). We found that the time series subsidence of two methods is in good agreement. The mean square errors (MSE) was 8 mm and 7 mm, respectively. Meanwhile, the linear regression analysis between InSAR and leveling results in 2012 to 2013 and 2017 to 2018 were shown that the correlation coefficient (R2, at a 95% confidence level) reached 0.97 and 0.95, which indicated that InSAR results were highly correlated with leveling values.

**Figure 11.** (**a**) and (**b**) Comparison between time series subsidence derived by InSAR and leveling results in two levelings, respectively, whose locations are shown in Figure 1 with levelings-I.

#### **5. Discussion**

#### *5.1. Land Subsidence Response to Di*ff*erent Land Use Types*

As seen in Figure 12 and Table 4, the land-use in the subsidence features of Shangzhuang, Beijing airport, Jinzhan, Heizhuanghu, Wangqingtuo and Shengfang is relatively complex. The maximum subsidence in Jinzhan, Heizhuanghu, Wangqingtuo, and Shengfang areas exceeds 100 mm/year, with 117, 105, 120 and 101 mm/year, respectively. In the Jinzhan and Heizhuanghu areas, the dense residential land and high population density in the area lead to a considerable water demand. Meanwhile, these areas are located in the groundwater funnel area of Beijing, which leads to large land subsidence. Since the MSWDP began operation in 2015, the urban water supply has been alleviated, and land subsidence has been slowed to some degree. As for Wangqingtuo and Shenfang areas, the dense population as well as agricultural and industrial lands there all result in an increasing demand of water and further aggravation of land subsidence. A large number of water-consuming plants and crops are planted in this region, and 85% of the well water is used for residents' lives and agricultural production [55]. Unlike the Wangqingtuo area, Shengfang mainly focuses on industry, mainly the steel industry, glass industry, and sheet metal industry. The steel industry, in particular, involves high water consumption industry. There are three large steel mills in the region, with a steelmaking capacity of 10 million tons. Industries with high water consumption may lead to the

serious development of land subsidence in this area. Similarly to Shengfang, in the Guangyang area, the main industrial categories are machinery manufacturing, building materials, electronic appliances, food, seasoning, clothing, glass products, printing, and chemical. In addition, the agricultural water demand is also large (winter wheat, maize, peanut, and cotton are mainly planted here). The crop area of Jinzhou town of Guangyang is one of the top 1000 crop sown areas in China. Dongguang, Jingxian, Zaoqiang, and Julu are the main grain production areas in Hebei. The main crops in the region are wheat and maize, which feature the largest total water consumption in BTH, with the areas of 680, 1293, 753, and 607 million km2, respectively [56]. The annual groundwater overdraft for crop irrigation ranges from 3 billion to 5.5 billion m3 [57]. A large amount of groundwater overexploitation results in serious land subsidence in the region. The serious land subsidence could affect the urban public facilities, transportation and buildings, such as the damage of urban underground pipe network, the instability of railway subgrade and cracks of buildings. The response characteristics of groundwater change and land subsidence are discussed in Section 5.2.

**Figure 12.** Remote sensing images in eleven subsidence features areas. (**a**–**k**) is the land use in Shangzhuang, Beijing Airport, Jinzhan, Heizhuanghu, Guangyang, Wangqingtuo, Shengfang, Dongguang, Jingxian, Zaoqiang and Julu.


**Table 4.** Main types of land subsidence in eleven subsidence features areas. The subsidence range represents the range of the average subsidence from 2012 to 2018 in each area.

#### *5.2. Land Subsidence Response to Di*ff*erent Water Use*

As seen in Figure 13, we found that the Beijing and Hebei supplied water mainly from the groundwater, while Tianjin mainly supplied surface water. In Beijing and Hebei, the total groundwater supply amount was 2.04 and 15.13 to 1.63 and 10.61 billion m3 from 2012 to 2018, and the total surface water supply amount was 0.8 and 4.13 to 1.23 and 7.04 billion m<sup>3</sup> from 2012 to 2018, respectively. However, in Tianjin, the total groundwater and surface water supply amount were 0.55 to 0.44 and 1.6 to 1.95 billion m<sup>3</sup> from 2012 to 2018, respectively. From the point of view of water consumption, Beijing needs more water for residential use, while Tianjin and Hebei need more water for agricultural use, especially Hebei, which used more than 10 billion m3 water for agricultural use from 2012 to 2018. The water for residential use in BTH showed an increasing trend, which ranged from 1.6 to 1.84 billion m3 in Beijing, 0.5 to 0.74 billion m<sup>3</sup> in Tianjin, and 2.34 to 2.78 billion m<sup>3</sup> in Hebei, from 2012 to 2018, with an increase in population. The population of Beijing, Tianjin, and Hebei increased from 20.69, 14.13, and 72.88 million to 21.54, 15.6, and 75.56 million, respectively. Due to major water transfer projects such as the MSWDP, the surface water supply in BTH shows an increasing trend, while the groundwater supply shows a decreasing trend from 2012 to 2018. Further investigations revealed that Beijing and Tianjin received 34.9 and 33.8 million m3 water from MSWDP as surface water for the urban water supply. In summary, after 2015, the exploitation of underground water decreased and the supply of surface water increased, which slowed the land subsidence in BTH.

**Figure 13.** Water resources utilization in BTH from 2012 to 2018. (**a**–**c**) is Beijing, Tianjin and Hebei, respectively. Total surface water supply refers to the water intake of surface water engineering, which are water storage, water extraction and water diversion. Total groundwater supply refers to the exploitation of water wells. Total industrial water refers to the water used by industrial and mining enterprises in manufacturing, processing, cooling, air conditioning, purification, washing and other aspects in the production process. Total agricultural water use includes farmland irrigation water, forest and fruit land irrigation water, grassland irrigation water, fish-pond recharge water and livestock and poultry water. Total residential water use includes urban domestic water and rural domestic water.

In Figure 14, the groundwater level shows a drop with a significant fluctuation trend, especially in wells h to j. The groundwater level is the lowest around May to August in wells h to j, possibly because these wells are located in Hengshui of the Hebei area, which is one of the national grain bases. The main crop in the area is winter wheat [57]. Because the irrigation period of this crop is February and March, and the precipitation is lower at this time, the groundwater exploitation volume is large, and the groundwater level drops significantly. After June, due to an increase in precipitation, the exploitation of groundwater decreased and the water level gradually rises. Furthermore, the comparison results of land subsidence changes and groundwater level changes reveal that the subsidence changes are consistent with the groundwater level changes, showing a seasonal trend. Meanwhile, we found that there is a time lag between land subsidence and groundwater level.

Furthermore, we applied the IRF method to investigate the lag time between subsidence change and groundwater level change. The ADF test and IRF results are shown in Table 3 and Figure 14, respectively. It can be seen from Figure 15 that groundwater level change has an obvious effect on land subsidence at each subsidence feature, and this effect has a lag effect. For the subsidence features of Shangzhuang, Heizhuanghu, Guangyang, Wangqingtuo, and the Dongguang area, the groundwater level changes in the third period (the third months) have the greatest impact on land subsidence, which means that there is a lag period of approximately three months between land subsidence and groundwater level change. However, for the subsidence features of Beijing Airport, Jinzhan, Jingxian and Zaoqiang, the groundwater level change in the second period (the second months) have the greatest impact on land subsidence, which means that there is a lag period of approximately two months between land subsidence and groundwater level change. Moreover, for the Shengfang subsidence feature, the groundwater level changes in the fourth period (the fourth months) have the greatest impact on land subsidence, which means that there is a lag period of approximately four months between land subsidence and groundwater level change.

**Figure 14.** Comparison between subsidence changes and groundwater levels in different groundwater monitoring wells. (**a**–**j**) The corresponding subsidence features from (**a**) to (**j**), which were previously discussed. The locations of the groundwater monitoring wells are shown in Figure 1.

**Figure 15.** Orthogonal impulse response from groundwater level change to subsidence in subsidence features of (**a**) Shangzhuang Beijing Airport, Jinzhan and Heizhuanghu, (**b**) Guangyang, Wangqingtuo and Shengfang, (**c**) Dongguang, Jingxian and Zaoqiang. The horizontal axis represents the lag order of the impact (month), and the vertical axis represents the response degree.

#### **6. Conclusions**

Because of the rapid growth of population and the development of industry and agriculture in BTH, the regional water demand is huge. However, the over-exploitation of groundwater in the BTH is obvious due to the shortage of water resources, which leads to the serious land subsidence in BTH. The serious subsidence in BTH resulted in the elevation loss, damage of roads, bridges, underground pipes and other municipal facilities, urban flooding and hydrops. In this study, the IPTA with SBAS-InSAR method was optimized to investigate land subsidence based on 126 Radarsat-2 and 184 Sentinel-1 images from 2012 to 2018 in BTH. We compared the InSAR-derived results with leveling-derived results from 2012 to 2013 and 2017 to 2018. The correlation coefficients (at a 95% confidence level) all exceed 0.9, indicating that our result is reliable. The land subsidence in BTH is very serious, with the maximum subsidence rate exceeding 120 mm/year. Eleven subsidence features in three subsidence areas were identified to be distributed in Beijing, Tianjin and Langfang, Baoding, Hengshui, and Xingtai, with maximum subsidence rates of 148, 145 and 131 mm/year, respectively. The area of subsidence greater than 50 mm in these areas was shown a decreasing trend after 2015. Meanwhile, the analysis of time-series subsidence among eleven features also shows a decreasing trend starting from 2015.

Furthermore, we discussed the relationship between land subsidence and different land-use types in subsidence features. We found that when the land-use type in an area is residential land or overlaps with agricultural and industrial land, land subsidence is relatively serious. For example, in Jinzhan, Heizhuanghu, Wangqingtuo, and Shengfang area, the maximum subsidence exceeds 100 mm/year. The distribution of different land use types leads to different water use structures. There is a large number of residential areas in Jinzhan and Heizhuang, with a large population and a correspondingly large water demand. In Guangyang, Wangqingtuo, and Shengfang, many high-water consumption industries and agriculture increase the regional water demand. However, a large amount of agricultural land is distributed in the Dongguang, Jingxian, Zaoqiang and Julu areas, and most of that land is cultivated wheat. Due to the lack of precipitation during the irrigation period of wheat, a large amount of groundwater exploitation leads to land subsidence in the region.

The relationship between land subsidence changes and groundwater level changes in the subsidence features was examined by using ten groundwater monitoring wells. We determined that the land subsidence changes are consistent with the groundwater level changes and show a seasonal trend. Nevertheless, there is a time lag between land subsidence changes and groundwater levels. The IRF method was used to investigate the time lag between subsidence changes and groundwater level changes. Our analysis shows that the lag time is two months in Beijing Airport, Jinzhan, Jingxian, and Zaoqiang, three months in Shangzhuang, Heizhuanghu, Guangyang, Wangqingtuo, and Dongguang, and four months in Shengfang.

**Author Contributions:** Conceptualization, B.C. and H.G.; methodology, M.G. and Q.C.; software, M.G., Q.C. and L.D.; validation, J.C., L.D. and M.S.; formal analysis, C.Z.; investigation, C.Z. and J.Z.; data curation, Q.C., J.C. and L.D.; writing—original draft preparation, C.Z.; writing—review and editing, J.Z. and M.S.; project administration, B.C. and H.G.; funding acquisition, B.C. and H.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China (No. 41930109/ D010702, 41771455/ D010702), China Postdoctoral Science Foundation (No. 2018M641407), Beijing Outstanding Young Scientist Program (No. BJJWZYJH01201910028032), Natural Science Foundation of Beijing Municipality (No. 8182013), Beijing Youth Top Talent Project, a program of Beijing Scholars and National "Double-Class" Construction of University Projects.

**Acknowledgments:** We thank the National Bureau of Statistics of China for the data of water consumption. The provision of IPTA developed by the GAMMA Company and the free R package "vars" created by Bernhard Pfaff and Matthieu Stigler are gratefully acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Localized Subsidence Zones in Gävle City Detected by Sentinel-1 PSI and Leveling Data**

#### **Nureldin A. A. Gido, Mohammad Bagherbandi and Faramarz Nilfouroushan**


Received: 12 May 2020; Accepted: 12 August 2020; Published: 14 August 2020

**Abstract:** Among different sets of constraints and hazards that have to be considered in the management of cities and land use, land surface subsidence is one of the important issues that can lead to many problems, and its economic consequences cannot be ignored. In this study, the ground surface deformation of Gävle city in Sweden is investigated using the Persistent Scatterer Interferometry (PSI) technique as well as analyzing the historical leveling data. The PSI technique is used to map the location of hazard zones and their ongoing subsidence rate. Two ascending and descending Sentinel-1 datasets, collected between January 2015 and May 2020, covering the Gävle city, were processed and analyzed. In addition, a long record of a leveling dataset, covering the period from 1974 to 2019, was used to detect the rate of subsidence in some locations which were not reported before. Our PSI analysis reveals that the center of Gävle is relatively stable with minor deformation ranged between −2 ± 0.5 mm/yr to +2 ± 0.5 mm/yr in vertical and east–west components. However, the land surface toward the northeast of the city is relatively subsiding with a higher annual rate of up to −6 ± 0.46 mm/yr. The comparison at sparse locations shows a close agreement between the subsidence rates obtained from precise leveling and PSI results. The regional quaternary deposits map was overlaid with PSI results and it shows the subsidence areas are mostly located in zones where the subsurface layer is marked by artificial fill materials. The knowledge of the spatio-temporal extents of land surface subsidence for undergoing urban areas can help to develop and establish models to mitigate hazards associated with such land settlement.

**Keywords:** deformation; geology; Gävle; InSAR; PSI; precise levelling; sentinel-1; subsidence; Sweden

#### **1. Introduction**

The management of cities and land use are challenging activities, where different sets of constraints, hazards, and preventive decisions have to be considered. Among the hazards, land surface subsidence and its subsequent effects on building and infrastructure can be considered as a major geo-hazard in many cities around the world. The principal causes of such deformation can be attributed to an aquifer-system compaction associated with groundwater/oil depletion, drainage of organic soils, underground mining, hydrocompaction, natural compaction, sinkholes, stress ended by a load of construction and foundation type, and thawing permafrost [1–4].

Globally, numerous cities and regions undergoing development are characterized by different rates of land movement phenomena, such as in Hong Kong [5], London [6], Las Vegas [7], and Tehran [3]. Therefore, knowledge of large-scale land surface subsidence could be an important factor in the long-term planning and managing of cities.

In recent years, Synthetic Aperture Radar interferometry (InSAR) technology has proved its capability in providing an accurate large cover land surface deformation monitoring information from SAR images [8–10], when compared for instance, with a traditional GNSS (Global Navigation Satellite System) and repeated precise leveling methods. InSAR and its differential method (DInSAR) use the phase difference between radar images acquired at different times and geometries, and allow the generation of digital elevation models and measurements of the centimeter-scale earth surface movements [11–13]. The DInSAR method is mostly used when measuring large deformations induced, for instance, by earthquakes or volcano activities. For small deformation, the method has considerable limitations, where the most important one is the atmospheric phase contribution, which can be addressed by applying the so-called Multi-temporal InSAR (MT-InSAR) technique [13–16]. The method looks through multi-temporal SAR data, even with long time intervals, for targets with relatively stable scattering properties and good coherence. Two other techniques were developed, which explored two different kinds of radar scattering targets. The Small Baseline Subset (SBAS), which looks for a distributed scatterer within the resolution cell [17], while the Permanent Scatterers Interferometry (PSI) looks for a point-like scatterer [13,18,19].

The PSI method exploits multiple SAR images obtained over the same area, and appropriate data processing and analysis procedures are used to separate the displacement phase from the other phase components [14,20]. Deformation time series and velocity estimation at measurement points (Permanent Scatterers:PS) are the main outcomes of the PSI analysis over the studied area. A large number of available images would result in a better PSI analysis e.g., 15–20 images for C-band satellite sensors [21]. The PSI method is applicable in many areas such as (i) urban, peri-urban, and build areas [17,22], (ii) geophysics, e.g., tectonics and volcanology [17,23], (iii) landslides [16,24–27], i.e., surface ground deformation monitoring that is due to subsurface clay deposits [28,29] and groundwater/oil depletion [3,30].

Numerous studies have been carried out using different SAR datasets, e.g., TerraSAR-X, ENVISAT, and Sentinel-1 images with SBAS and PSI methods to investigate ground surface deformation associated with different causes (e.g., groundwater depletion, landslides, surface geology, etc.). For example, Haghshenas and Motagh [3] presented the results of (InSAR) analysis of Tehran using different SAR data between 2003 and 2017, where they could identify three distinct subsidence features attributed to groundwater extraction. Zhou et al. [31] assessed the spatial–temporal analysis of land subsidence in Beijing using a small baseline subset (SBAS) InSAR based on 47 TerraSAR-X SAR images from 2010 to 2015, where distinct variations of the land subsidence were found in the study regions. Terrain deformation in the metropolitan area of Sydney city was investigated by Ng et al. [29] using 49 sets of C-band ENVISAT radar images acquired between 2003 and 2010 and the PSI technique. The results reveal the relative stability of the city with major deformation ranged between −3 mm/yr to +3 mm/yr. In another study, Fryksten and Nilfouroushan [28] mapped the ongoing deformation in Uppsala City using PSI analysis and Sentinel-1 radar datasets covering the period of 2015 to 2019. Their results show that the city was undergoing significant subsidence in some areas that were dominated by postglacial clay, reaching to −6 mm/yr in the Line of Sight (LOS) direction. Furthermore, Hu et al. [5] generated surface deformation map for Shenzhen, China, and Hong Kong using PSI and SBAS methods in addition to GNSS observations. Significant subsidence was revealed in the land reclamation area, in the metro construction area, and in the building with a shallow foundation. Recently, Goorabi et al. [32] studied the land subsidence and its relationship with geological and geomorphological in Isfahan metropolitan using the sentinel-1A InSAR observation and PSI technique employed by SARPROZ software. More similar studies can be read in, e.g., Keren [33], Gernhardt, and Bamler [34], Sousa et al. [35], Lan et al. [36], Dehghani et al. [37], Béjar-Pizarro et al. [38], Casagli et al. [39], and Barra et al. [40].

Since towns and cities became the focus for regeneration and developments, up-to-date ground information (e.g., related to characteristics of land subsidence and ground stability) for the urban environment became very essential, where the properties for an anthropogenic and natural superficial deposits need to be defined and addressed. The quality of such information is of great importance for different purposes, such as urban spatial planning and development activities, environmental management, geo-hazard monitoring, and assessment efforts. For instance, the PanGeo project (https://www.eurogeosurveys.org/projects/pangeo/), which was a geo-hazards European project, used satellite radar data integrated with local geological information to provide geo-hazard information for the largest 52 cities across Europe, including Stockholm and Gothenburg in Sweden. Therefore, the integration of the radar satellite data with the local in situ data (e.g., geology and precise leveling data) would be of great importance for the geo-hazard mapping.

Generally, land subsidence impacts in urban regions are quite numerous and can be classified into infrastructural, environmental, economic, and social impacts. In Gävle, which is one of the old coastal growing cities in Sweden, the subsidence impacts can be seen in forms of cracking and damage of buildings and infrastructure, especially in the northern part of the city, which negatively will affect the rehabilitation and maintenance cost and the quality of the living environment. For instance, according to the latest surveys and reports made by the Gävle municipality (www.gavle.se), some streets in the northern part of the city (i.e., Norra Strandgatan) were repaired due to the subsidence problems. Furthermore, visible related problems probably due to subsidence of the buildings in the Gävle Strand area, have been reported by, e.g., WSP and SWECO Companies. Additionally, subsidence for 13 benchmark points in the area of Alderhomen, around Gävle beach, ranging between 4 and 12 mm, have been reported by Lindén and Rosenqvist [41], using leveling data collected between 2006 and 2011. Hedberg and Ottekrall [42], have reported similar problem for three different buildings located in the city center using leveling data between 1985 and 2011.

The city of Gävle crosses the Gävleån River that formed a large delta at the outlet in the Gulf of Gävle. The sub-surface geology map of the city shows areas marked with artificial fill that have been used for urban development activities. According to the Department of the Environment, Transport, and the Regions, and Environment Agency in the UK (DETR/EA, 1998) "the fill material at any site is commonly an admixture of organic, chemical and inert material which can lead to serious problems in the built environment". Thus, the ongoing and expected development of the city according to the municipality of Gävle, along the coastal line, within the formed delta and the artificial fill area, is prone to subsidence hazard and needs further investigation with, e.g., satellite remote sensing (InSAR) for deformation mapping. According to the Gävle municipality, there are plans for numerous large construction projects; e.g., Näringen, which is one of the largest urban transformations in the municipality and targeted to become one of Europe's most sustainable districts with 6000 new houses and space of 450,000 m<sup>2</sup> of gross area businesses. The targeted area is located around and within the identified geo-hazard zone areas according to the previous studies. Moreover, the availability of long precise leveling observation records that last for about 45 years in monitoring the stability of some locations in the city, in addition to other information from remote sensing data would help and motivate us to further study and map the localized deformation in the city of Gävle.

The aim of this study is to generate a precise deformation map of Gävle city, to detect and highlight the hazard zones, which relatively deform with a higher rate in this city, and correlate the results with subsurface geology. In this study, Sentinel-1 satellite radar datasets provided by the European Space Agency (ESA) have been processed and analyzed using the PSI technique, which is essentially an improved version of the standard PS-InSAR algorithm, developed by Perissin and Wang [43] and implemented in the SARPROZ software. For this purpose, two freely available datasets of Sentinel-1 radar images in ascending and descending geometries were used, that covered more than a five-year period starting from January 2015 to May 2020. Available long records of precise leveling data, back to early 70s, observed and processed by the municipality and the University of Gävle, were utilized to examine the obtained PSI results for the Sentinel-1 observation period as well as to analyze the long-term stability of some locations in the city which were not reported before.

#### **2. Study Area**

The study area covers about 24 km2 of Gävle city which is located by the Baltic Sea in east-central Sweden, at about 60 degrees north and 17 degrees east (see Figure 1). Gävle city, which is the capital of Gävleborg County is the oldest in the historical Norrland, with a population of about 102,000 in 2019, which makes it the 13th largest city in Sweden. According to the quaternary deposit map of the study area (Figure 2), provided by the Geological Survey of Sweden (SGU), the surface geology of the city mainly consists of post-glacial sand which is dominant and covers most of the city center and extends toward the west. A glaciofluvial sedimental strip crosses the city from southwest to the northeast. Glacial clay and Post-glacial clay layers located toward the south, the east, and partially to the north of the city. In addition, some areas to the northeast are of artificial fill material. The two banks of the outlet part of Gävleån River (i.e., about 6 km across the city) that crosses the city from the west to the east toward the Baltic Sea, consist of artificial fill material and have been constructed. The fill areas were part of the original route of the river canal and extend toward the north from the riverbank. Furthermore, the outer areas of the city are almost sandy till layers type. According to the SGU (https://www.sgu.se/) the estimated soil depth in the city ranges between 3 to more than 50 m in some areas along the coastline, and between 10 and 20 m in the artificial fill areas. Figure 1 shows the location map of the city, urban zone, and the area covered by leveling measurements (the red square in Figure 1b).

**Figure 1.** Shows the study area of Gävle city (**a**) location map (**b**) map of the urban zone and (**c**) locations of the four buildings (B-1, B-2, B-3, and B-4) which have been measured and monitored using precise leveling between 1974 and 2019. Coordinate system: SWEREF 99 TM.

**Figure 2.** The quaternary deposit map of the study area © Geological Survey of Sweden (SGU). Coordinate system: SWEREF 99 TM.

#### **3. Data and Methods**

#### *3.1. PSI Method*

According to Crosetto et al. [19], the comprehensive DInSAR phase components can be read as below, if a digital elevation model (DEM) of the imaged scene is available (i.e., the topographic term can be simulated and subtracted from the interferometric phase Δϕint yielding the so-called DInSAR phase components Δϕ*D*−int):

$$
\Delta q\_{\rm D-int} = \Delta q\_{\rm Disp} + \Delta q\_{\rm Top-res} + \Delta q\_{\rm atm-S} - \Delta q\_{\rm atm-M} + \Delta q\_{\rm orb-S} - \Delta q\_{\rm orb-M} + \Delta q\_{\rm noise} + 2\pi n \tag{1}
$$

where Δϕ*Disp* is the displacement phase term in the LOS direction, Δϕ*Top*−*res* represents the residual topographic error components, Δϕ*atm* is the atmospheric phase components at an acquisition time of each image (Master and Slave), and it can be reduced by applying stacking technique [44], Δϕ*orb* is the orbital error phase components of each image (Master and Slave), and it can be estimated using precise orbits [45], Δϕ*noise* is the phase noise which can be minimized with filtering and multi-look techniques [45,46]. The last term of Equation (1) represent the general 2π ambiguity (integer value *n*) is a result of errors in the phase unwrapping process, which is necessary to solve for this ambiguity, i.e., the DInSAR phases are bounded in the range (−π, π) [11,47]. The goal of any DInSAR technique is to estimate the displacement phase components Δϕ*Disp* from the Δϕ*D*−int, which requires its separation from the other components in Equation (1).

To overcome the limitations of the deformation monitoring using the DInSAR (i.e., temporal and geometrical decorrelation and atmospheric disturbances), Ferretti et al. [13,14] developed the first PSI technique based only on a subset of pixel analysis. It is called permanent scatterers (PS) and assumed to possess stable phase values during the time series of the whole used datasets of the same area. From these datasets, a single master is selected considering a high coherence in all formed

interferograms over the whole datasets. Temporal sampling of the data should be regular to some extent. Prior to the PS estimation, a DEM of the area should be used to remove the topographic phase components [19]. Accuracy of about 20 m is sufficient for the used DEM [14].

The PS is the strong reflecting objects that are dominant in a cell and it can be identified based on the amplitude stability index and the temporal coherence [48]. Stable pixel selection depends on the use of amplitude dispersion index [49], which can be computed as in Ferretti et al. [50]:

$$D\_A = \frac{\sigma\_A}{m\_A} \tag{2}$$

where σ*<sup>A</sup>* and *mA* are the standard deviations and the mean of the amplitude, respectively. Therefore, a pixel with high reflection has a low *DA*. Temporal coherence explains the mean difference between the observed and modeled phase for each pixel and interferogram [9]. Over the selected PS points, the atmospheric phase components should be estimated for every interferogram, where the interpolated map is called the atmospheric phase screen (APS). According to Colesanti et al. [51], at least 25 images are needed in addition to a PS density of about 5 to 10 PS/km2 to achieve proper estimation and removal of the APS effects. Furthermore, a deformation model has to be introduced prior to the velocity estimation e.g., linear, nonlinear, and seasonal models [52].

It is worth mentioning here that, the PSI technique [13,14] was followed by several contributions and new approaches such as PSI procedure for crustal displacement in non-urban environments [23], the adaptation of the GPS LAMBDA method to PSI [53], the SqueeSAR algorithm, which is an extension of the original PS-InSAR algorithm [54], and the PSI approach introduced by Perissin and Wang [43], which is called the QPS algorithm. The PSI method by Perissin and Wang (i.e., QPS) exploits partially coherent targets, and it has been incorporated in the SARPROZ InSAR processor software which is used in this study.

#### *3.2. Sentinel-1 Data and Analysis*

In this study, freely available C-band Sentinel-1 A and B data from the European Space Agency (ESA) were downloaded (https://scihub.copernicus.eu/dhus/#/home) and used to detect and quantify the ground surface deformation of Gävle city, using the PSI method [18,51]. Two datasets consisting of 41 and 50 Single Look Complex (SLC) images on ascending and descending geometries were processed and analyzed over the study area covering periods from January 2015 to May 2020 and from June 2015 to May 2020, respectively (see Table 1). The temporal resolution for the acquired images was about one month, with some exceptions due to the unavailability of some images with certain characteristics or some images were removed due to the winter season which may reduce image coherence.


**Table 1.** Details of the Sentinel-1 A and B datasets used for the PSI time series analysis and their properties.

For the data processing, SARPROZ software [49] was used, which could determine the LOS rate and its associated uncertainty for each individual PS point. The SLC images were co-registered to a single master. Then, the 50 m grid local DEM, provided by the National Land Survey of Sweden (Lantmäteriet), and the precise orbits for each image were used to remove the topographic and flat earth components, respectively. Reference points were selected among the selected Persistent Scatter Candidates (PSCs), which are relatively unaffected by deformation. For the PS point selections, a quality factor called amplitude stability index is used and defined as 1 − *DA* in the software.

The first and second sets of data i.e., the 41 ascending and 50 descending images, were used to analyze the land surface of the study area. We used the images dated by 2016-11-06 (for ascending) and 2017-11-13 (for descending) as master images, respectively, according to certain characteristics including the length of the baseline, where a high stack coherence is achieved.

The APS effect was estimated and removed in the location of the PSCs, which were selected based on amplitude stability indexes of 0.7 and 0.8 as quality parameters for the ascending and descending datasets, respectively. Temporal coherence of 0.60 was used to select all persistent scatterers (PS) points. For each dataset, the line of sight (LOS) velocity and displacement time series were estimated utilizing a linear model. To ensure reliable results, the two reference points (located in Sätra, Figures 3 and 4) with zero displacements, were selected in the same very stable area and away from the known deformed zones. Furthermore, the estimated LOS displacement velocities were decomposed to the east–west and vertical components by combining ascending and descending results using the following Equations (Equations (3) and (4)) presented by Milillo et al. [55] and implemented and used in SARPROZ software.

$$d\_{\rm East-West} = \frac{1}{2} \left[ \frac{d\_{\rm Dcs}}{\sin(\theta\_{\rm Dcs})} - \frac{d\_{\rm Asc}}{\sin(\theta\_{\rm Asc})} \right] \tag{3}$$

$$d\_{Vcr} = \frac{1}{2} \left[ \frac{d\_{\rm Des}}{\cos(\theta\_{\rm Des})} - \frac{d\_{\rm Asc}}{\cos(\theta\_{\rm Asc})} \right] \tag{4}$$

where θ and *d* are the incidence angle and the estimated displacement velocity, respectively, in the ascending (*Asc*) and descending (*Des*) geometry. AD pair (i.e., Ascending Descending pair) method was implemented by using 10 m as a maximum offset distance (in planar and height) between the selected pair.

**Figure 3.** Line of Sight (LOS) displacement rates of all generated ascending PS points in Gävle city relative to the reference point (pink color). Area-1 shows the maximum displacement zone. A1–A5 refer to the five different small areas. The observation period is from January 2015 to May 2020. Coordinate system: SWEREF 99.

**Figure 4.** LOS displacement rates of all generated descending PS points in Gävle city relative to the reference point (pink color). Area-1 shows the maximum displacement zone. A1–A5 refer to the five different small areas. The observation period is from June 2015 to May 2020. Coordinate system: SWEREF 99.

#### *3.3. Levelling Data and Analysis*

A long record of leveling measurements has been conducted and processed by the municipality of the city of Gävle and partially by the University of Gävle, covering the period from 1974 to 2019. The data were used to examine the obtained PSI results as well as to analyze the ground stability of some locations in the study area. The record consists of four short separated leveling lines ranging between 100 and 500 m and has been established by the municipality to monitor four different buildings located in the center of the city (see Figure 1c) within the vicinity of the Gävleån River. Each line consists of several metal pegs that mounted in the building foundation and associated with a certain benchmark (BM). A double run leveling procedure has been performed using different types of leveling instruments including Leica DNA03, NA2, and NA720. The measurements are referenced to the official national height system RH 2000 in Sweden. Periods of observation records for each building varies and

ranges between 6 years and 34 years. For instance, the levelling records for Building-1 ranged from 1985 to 2019, Building-2 ranged between 2000 and 2019, Building-3 ranged between 1976 and 1982, and Building-4 between 1974 and 1988 (Figure 1c). In addition, the temporal resolution for repeated leveling observations also varies from months to several years (this will be discussed in Figure 6). By using regression analysis, the annual rates of vertical changes have been computed at each peg.

#### *3.4. Geological Data and Analysis*

According to the National Geotechnical Institute "Statens Geotekniska Institute" (SGI) assessment report of Gävle municipality [56], the dominant soil along the coastal stretch of the Gävle municipality is almost blocky and large-block moraine, while areas with sand and fine sediments, such as clay and silt, are found in areas near to watercourses and small lakes as well as in sea bays (cf. Figure 2). Therefore, based on the soil condition, the report highlighted five areas in the region, attributed to the hazard of landslide, race, and erosion [56]. The identified five areas almost consist of silt and/or clay as well as fill that probably underpinned of fine sediment and adjacent to the coastal lakes' watercourses. The fill term in this study is used to describe ground that has been formed by human activities in which natural soils or sometimes man-made materials are deposited, in contrast to natural soils which has its origin in geological processes [57]. At the outlet of the Gulf of Gävle, the Gävleån River, which crosses the city, formed a large delta, which has been covered using fill material accordingly and used for urban development. In addition, the associated small river arms are also channeled or have been filled.

By combining the geological information of Gävle city with the resulted PSI deformation map, it is possible to assess the hazard of land movement in the city. The combined ascending and descending PSI results (i.e., vertical components) and their rate of movement were analyzed based on the provided 1:25000–1:100000 quaternary deposit map by the SGU (Figure 2). The PS points were overlaid on the geological map where the rate of displacement can be evaluated accordingly (cf. Figure 7).

#### **4. Results and Discussions**

In this section, we report our results of the estimated land subsidence rate over the city of Gävle and the relation between the obtained PSI results with the available leveling data and the geological information.

#### *4.1. PSI Results*

The two Sentinel-1 datasets which consist of ascending and descending images were analyzed over the city, where the PSI technique was applied to detect relative ground deformation, using the processing criteria in Section 3.2. By using 0.6 as a temporal coherence mask, a total of 7265 and 9314 PS points were obtained using ascending and descending geometries, respectively, by employing SARPROZ software.

Generally, the obtained ascending and descending results show good agreement in the LOS direction relative to selected references, which confirm the stability of the two references that have been selected in the same location (see Figures 3 and 4). It is worth noting that, the selection of the reference points has to be very careful due to its effects on all PS points. The selected two reference points (ID 5776 and ID 1976) have zero velocity and 0.94 and 0.99 temporal coherence, respectively, and are located in a very stable part of the city (Sätra area) that consist of types of sandy till layers.

Figure 3 shows all generated ascending PS points, the main deformed areas (Area-1 in the white square), and five different small areas (i.e., A1–A5) where the PS rates have been averaged to generate time series for the selected five areas using 5 to 10 PS points (see Figure 5a), in addition to the location of the reference point. Point selection for averaging is based on the closeness of the points to each other (i.e., within 100 m radius) and having the minimum coherence of 0.75. The identified hazard zone Area-1, which is located toward the northeast of the city center, shows the maximum negative LOS movements in the city, based on the estimated PS points in the area which reaches up to −6 mm/yr. Some sparse areas and points showed a notable subsidence rate ranging between −1 mm/yr and

−3 mm/yr, e.g., in the railway lines (area A5) to the northwest of Area-1 (Figure 3). The rest of the land surface in the city is almost relatively stable, where the maximum movement is ranged between about −2 mm/yr and +2 mm/yr. The small and sparse positive LOS movements may be attributed to ground heave (i.e., when the soil is saturated by water or accumulated due to construction activities).

**Figure 5.** Time series of the selected five areas (i.e., A1–A5), see Figures 3 and 4, where the permanent scatterers (PS) rates have been averaged for (**a**) the ascending, (**b**) the descending, and the references ID 5776 and ID 1976, respectively. The period is between January 2015–May 2020 and June 2015–May 2020.

Figure 4 shows all the generated descending PS points and the area with maximum subsidence rate (Area-1), and similarly the five different small areas (i.e., A1–A5), in addition to the location of the reference point. In both cases (i.e., ascending and descending results), the density of the obtained PS results is relatively poor (about one PS for each 2500 m2), based on the medium spatial resolution of the SLC Sentinel-1 radar data which is 5 m × 25 m. In contrast, more dense PS points would be expected when using high-resolution images e.g., X-band images, which offer a resolution up to 1.1 m × 0.6 m as for TerraSAR-X [58].

Figure 5 shows the time series and the accumulative LOS relative movement in millimeter for the five ascending and descending small areas (i.e., A1–A5), and the reference points ID 5776 and ID 1976, respectively. There are a lack of data during periods of 31 October 2015 to 15 June 2016 (about 7 months) and between 12 January 2018 to 18 April 2018 for the descending dataset, in addition to different periods in the ascending set, ranged between 2 and 5 months (Figure 5), despite the short revisit time of Sentinel-1 sensors which is 6–12 days. The results show a continuous subsiding rate with a notable deviation in e.g., January, and February which can be attributed to the winter season. The reference points (ID 5776 and ID 1976) demonstrate a very stable and consistent time series with a minor deviation in some months e.g., February.

#### *4.2. Leveling Data And the PSI Results*

By utilizing the leveling records associated with the four buildings (Figure 1c), the annual rate of change at each leveling peg has been computed using the regression analysis method. Figure 6 illustrates the time series and the computed rate of change at a selected peg for each building. The selection of the peg point is based on its closeness to the generated PS points on or very near to the building to facilitate the validation process. It is worth noting that, the leveling records of the Buildings 1 and 2 (Figure 6a,b) extend for a long period (1985 to 2019) and overlap with the satellite radar data (2015 to 2020). In contrast, the leveling records of the Buildings 3 and 4 (Figure 6c,d) ended at the years 1982 and 1988, respectively, i.e., there was no overlap between the leveling records and the SAR data. However, even that old leveling records for Buildings 3 and 4 in that period was used here to infer if the present-day subsidence rate from PSI is comparable with the one obtained from old leveling records at different time period.

**Figure 6.** The computed annual displacement rate in mm/yr on the selected leveling pegs, relative to four different benchmarks in the four selected buildings (**a**–**d**), according to Figure 1c.

For validation purposes, the available leveling information data were used to confirm the obtained PSI results, where the generated ascending and descending PS points on or close to the buildings and the leveling pegs were considered.

The LOS movement rates of all generated PS points from the PSI analysis are relative to one stable reference point, ID 5776 for ascending and ID 1976 for descending, while the computed displacement rates of the metal pegs, mounted in the buildings are referenced to specific benchmarks (BM), and these benchmarks could be moving with respect to the stable reference of the PSI analysis. To overcome this relative problem, the trend rate in the nearest PS points to each peg in the building will be computed relative to the nearest PS point to the leveling benchmark. Then, a comparison between the computed leveling rate (i.e., between the BM and the peg) and the relative rate between the PS points that represent the BM and peg (i.e., between the two PS points near to the BM and to the peg) can be performed. It is worth noting also, that the PSI results are in the LOS direction, while the leveling rate results in the vertical direction. To compensate for such differences, the PSI results (i.e., the displacement rate and the accumulative displacement) have been converted to vertical and horizontal (only possible for east–west) movements from the LOS direction by combining the ascending and descending datasets according to Milillo et al. [54] (see Figures 7 and 8 and Equations (3) and (4)). Table 2 shows the track type i.e., ascending and descending pair (AD), point number (ID), the relative vertical displacement rate in mm/yr, the vertical accumulative displacement in mm, and the coherence of the generated combined PS points on or near to the four buildings using about five years' radar data. In addition, the corresponding computed annual rate of the leveling measurements (Pre. Lev) in mm/yr for the four buildings, using different records ranged between 1974 and 2019, have been presented.

**Figure 7.** The vertical displacement rate of the combined ascending and descending PS points overlaid on the quaternary deposits map of the city of Gävle. Negative values denote for subsidence, while the positive denotes for uplift. Base map: Quaternary deposit © Geological Survey of Sweden (SGU). Coordinate system: SWEREF 99. Unit: mm/yr.

**Figure 8.** The east–west displacement rate of the combined ascending and descending PS points overlaid on the Topographic map of the city of Gävle. Negative values denote for west movement while the positive denotes for east movement. Coordinate system: SWEREF 99. Unit: mm/yr.

**Table 2.** Comparison between the computed precise leveling rates (Pre. Lev) using four different leveling records relative to four different BMs, and the relative vertical rate of the generated combined PS points in the four validation buildings using the ascending and descending combination of the 41 and 50 SAR images, respectively, collected between 2015 and 2020.


According to Table 2, the relative vertical PS results of point (ID 498) on Building-1 shows close agreement with the computed leveling rate at Point 7 (nearest combined PS point to the leveling point) i.e., −0.9 mm/yr vs. −1.2 mm/yr, respectively. In Building-2, the relative vertical results of point (ID 430) shows very close agreement in terms of the rate and the accumulative displacement with the leveling data at Point 6. In Building-3, despite the relativity short and old record of the leveling data (from 1976

to 1982) when compared to radar data (from 2015 to 2020), a lesser agreement is revealed in terms of the displacement rate between the leveling records at point 6 and the relative vertical combined PSI result at point (ID 396), which is −2.0 mm/yr and −0.6 mm/yr, respectively. In Building-4, the vertical relative combined result shows almost complete stability (i.e., zero rates) with a very minor cumulative displacement for the 5 years (−1.6 mm), while the historical leveling record (from 1974 to 1988) shows relative displacement rate of −1.8 mm/yr and accumulative displacement of −31.0 mm for a total of 14 years (i.e., there is no agreement). The result reveals that Building-4 was experiencing a high relative subsidence rate during the years after construction phase, due probably to the load and foundation type, and reached complete stability during the last years.

The above comparison between the obtained relative PSI results in the vertical direction, using combined ascending and descending data for about 5 years, and the precise leveling information record for different periods at the four buildings reveals close agreement in Buildings 1 and 2.

#### *4.3. Subsidence and Subsurface Geology E*ff*ect*

By correlating the geological information of Gävle city with the resulted combined (ascending and descending) PSI deformation map of the land surface, it is possible to assess the vertical land movement in the city. The combined measured PS points and their rate of the vertical movement were analyzed based on the provided 1:25000–1:100000 quaternary deposit map by the SGU (Figure 2). The combined PS points were overlaid on the geological map where the high rate of displacements can be illustrated accordingly (Figure 7). Moreover, the resultant horizontal movement in the east–west direction has been overlaid on the topographic map of the city (Figure 8).

Based on our results, the Area-1 (see Figures 3, 4 and 7) shows the PS points with maximum displacement rates that reach up to −6 mm/yr. It is clearly visible from the corresponding geological map that the subsiding area belongs to the artificial fill zone that covers the banks of the outlet of Gävleån River and extends toward the north [59]. Some sparse PS points with a relatively high rate of subsidence are located on areas marked with clay to silt deposits (pale yellow solid zone in Figure 7) in the southeast of the city. The land surface of the central part of the city is relatively stable and it is almost located on the post-glacial sand deposits. However, some sparse PS points, mostly along and around the river, show a westward movement about 1–2 mm/yr relative to the selected reference point (Figure 8). This can be further investigated in future by analyzing the higher resolution radar satellite images, e.g., TerraSAR-X, which provides deformation maps with denser PS points.

Field investigation reveals that most of the construction buildings in areas A1 and A4, those located in the artificial fill Area-1, are of light type (e.g., warehouses and workshops). Meanwhile, area A2, which is the newest build area in the city, contains multi-story buildings. Recently, visible related problems probably due to subsidence of the buildings in area A2 has been reported to construction consultancy companies. Additionally, it is worth noting that, the fill area to the northeast of Area-1, which contains a large number of huge storage tanks, shows very stable PS results, with minor vertical movement of order of −1.0 to −2.9 mm/yr concentrated in the railway line in the vicinity of the oil storage tanks area, while the storage tanks area shows vertical movements of order less than −1.0 mm/yr. The storage tank foundations are probably built with special care to the subsurface geology type in the area. Furthermore, uneven subsidence evidence has been reported in the church building in the city center (Heliga Trefaldighetskyrkan) to the east of Building-1, where the very visible cracks can be attributed to different reasons such as building foundation settlement or age-related settlement or the close location to Gävleån River. However, due to the resolution of the PSI results, it was not possible to get enough PS points on this building that could be used for evaluation.

#### **5. Conclusions**

In this paper, we have carried out a comprehensive study to generate an accurate deformation map of Gävle city, to detect and highlight the hazard zones. For this purpose, two datasets of Sentinel-1 radar images in ascending and descending geometries, provided by the ESA, and covered the period of

January 2015 to May 2020, were processed and analyzed using the PSI method and combined ascending descending datasets. Available long records of precise leveling data were utilized to validate the PSI results. Furthermore, a correlation analysis was performed between the PSI results and the available geological maps for the study area.

Our results show that land surface in Gävle city is relatively stable with exceptions of some localized subsiding zones. There is no obvious pattern of deformation observed for the whole or a big part of the city. Only localized deformation zones were detected by our analysis that have been subsiding in the last five years. The correlation of our results with subsurface geology types, i.e., the artificial fill, suggests that the subsidence occurs due to relatively weak subsurface layers which either are affected by loading of new constructions or due to hydrocompaction. A further analysis is suggested to analyze the hydrogeological and geotechnical measurements and correlate them with InSAR results to better understand the ongoing land deformation process in this city.

The detected subsiding zones show a maximum displacement rate that reaches up to −6 ± 0.46 mm/yr in the LOS direction. To the northwest of Area-1 (i.e., in Gävle railway station's marshalling yard), the area A5 is experiencing average subsidence ranging from −1 ± 0.47 mm/yr to −3 ± 0.53 mm/yr. The rest of the city, including the city center, is almost relatively stable with minor deformation rates ranged between −2 ± 0.5 mm/yr to +2 ± 0.5 mm/yr in vertical and east–west components. A comparison between the computed precise leveling rates using the available four different leveling records relative to four different BMs, and the relative vertical rate of the generated PS points in the four validation buildings, shows close agreement in Buildings 1 and 2. According to the correlation analysis between the obtained PSI results and the geological information of the city, the achieved result is highly correlated with the quaternary deposit map, e.g., in the detected hazard zone that is reported as an artificial fill area (see Figure 7). This conclusion can be considered by the Gävle municipality for their future projects, e.g., in Gävle harbor.

Knowledge of the spatial and temporal extents of land surface subsidence for cities that are vulnerable and undergoing development can help to develop and establish measures to mitigate hazards associated with such land settlement.

**Author Contributions:** N.A.A.G. wrote the manuscript and processed and analyzed the data. M.B. designed the study and helped to write and prepare the manuscript. He also contributed to historical and precise leveling data processing. F.N. contributed to supervision and writing the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** We are immensely grateful to J. Gust. Richert Stiftelse foundation (Project number 2019-00485) for their financial support for this study.

**Acknowledgments:** The authors are grateful for the free use of the Sentinel-1 data provided by the European Space Agency (ESA), and to the SARPROZ Company for providing a free license for SARPROZ software for this study. We are also immensely grateful to the municipality of Gävle for providing us with the long record of leveling data. We are grateful to four anonymous reviewers who provided very useful and constructive reviews.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Technical Note*

## **Very Local Subsidence Near the Hot Spring Region in Hakone Volcano, Japan, Inferred from InSAR Time Series Analysis of ALOS**/**PALSAR Data**

#### **Ryosuke Doke, George Kikugawa and Kazuhiro Itadera**

Hot Springs Research Institute of Kanagawa Prefecture, 586 Iriuda, Odawara 250-0031, Japan; kikugawa@onken.odawara.kanagawa.jp (G.K.); itadera@onken.odawara.kanagawa.jp (K.I.) **\*** Correspondence: r-doke@onken.odawara.kanagawa.jp; Tel.: +81-465-23-3588

Received: 28 July 2020; Accepted: 31 August 2020; Published: 1 September 2020

**Abstract:** Monitoring of surface displacement by satellite-based interferometric synthetic aperture radar (InSAR) analysis is an effective method for detecting land subsidence in areas where routes of leveling measurements are undeveloped, such as mountainous areas. In particular, InSAR-based monitoring around well-developed hot spring resorts, such as those in Japan, is useful for conserving hot spring resources. Hakone Volcano is one of the major hot spring resorts in Japan, and many hot spring wells have been developed in the Owakudani fumarole area, where a small phreatic eruption occurred in 2015. In this study, we performed an InSAR time series analysis using the small baseline subset (SBAS) method and ALOS/PALSAR scenes of the Hakone Volcano to monitor surface displacements around the volcano. The results of the SBAS-InSAR time series analysis show highly localized subsidence to the west of Owakudani from 2006–2011 when the ALOS/PALSAR satellite was operated. The area of subsidence was approximately 500 m in diameter, and the peak rate of subsidence was approximately 25 mm/year. Modeling using a point pressure source suggested that the subsidence was caused by a contraction at approximately 700 m above sea level (about 300 m below the ground surface). The rate of this contraction was estimated to be 1.04 <sup>×</sup> 104 m3/year. Hot spring water is collected from a nearby well at almost the same depth as the contraction source, and its main dissolved ion component is chloride ions, suggesting that the hydrothermal fluids are supplied from deep within the volcano. The land subsidence suggests that the fumarole activity is attenuating due to a decrease in the supply of hydrothermal fluids from deeper areas.

**Keywords:** InSAR; ALOS/PALSAR; time series analysis; SBAS; land subsidence; hot spring water; hydrothermal fluids; Hakone Volcano

#### **1. Introduction**

Land subsidence due to the over-extraction of groundwater needs to be carefully monitored as it can lead to the destruction of buildings and infrastructure in urban areas. In Japan, land subsidence is monitored by leveling [1]. Although leveling measurements can be used to precisely determine vertical displacements, they can only be used to monitor displacements if benchmarks have been installed; if local subsidence occurs in areas away from benchmarks, it may be overlooked. Consequently, conventional leveling measurements are not very effective for detecting land subsidence in mountainous areas where leveling routes are not well developed. There are numerous volcanoes in Japan, many of which have hot springs and are popular tourist destinations. Large-scale resort development has occurred in these mountainous areas, which has meant that extensive drilling for new sources of hot spring water has occurred. Therefore, it is important to monitor land subsidence in mountainous areas with many large spa resorts, not only to maintain the infrastructure, but also to conserve hot spring resources.

Interferometric synthetic aperture radar (InSAR) analysis, using satellite-based sensors, is one of the most effective tools for detecting unknown land subsidence [2]. While InSAR is less precise than conventional leveling measurements, the technique does not require a network of ground-based observations. In addition, InSAR analysis is well suited for application in mountainous areas, where there is no existing observation network. SAR satellites were launched by governments and private companies, primarily to monitor the natural and physical characteristics of the Earth, including vegetation, natural resources, and crustal deformation for disaster monitoring. Various wavelengths are used for land surface observations, ranging from the L-band to the X-band, depending on specific requirements. InSAR analysis using short-wavelength sensors (e.g., X-band and C-band) generally does not provide good coherence in forested areas, which makes it challenging to monitor ground movement. Conversely, InSAR analyses using long-wavelength (L-band) sensors can achieve good coherence, even in forested areas, enabling users to monitor ground deformation. It is therefore considered that monitoring by L-band sensors will play an important role in clarifying ground deformation in mountainous areas in places such as forested area of Japan.

The Hakone Volcano (Figure 1), located in the plate collision zone of the Izu Peninsula with Japan mainland [3], has been active for approximately 400,000 years and is one of several active volcanoes in Japan. Hakone is a caldera volcano with a complex of stratovolcanoes of basaltic to andesitic composition forming the caldera rim (≥230 ka) (Figure 2) [4]. Within the last 37,000 years, andesitic eruptive activity has formed a central cone of the volcano [4]. The latest magmatic eruption occurred about 3000 years ago on the northern slope of the central cone and caused a sector collapse that formed an amphitheater [4]. Owakudani, the largest fumarole area of the volcano, is located at the eastern margin of the amphitheater (Figure 3). In the Owakudani, hot spring water is artificially produced from the self-discharged volcanic steam that comes up from borehole (steam wells), which are several hundred meters in depth. Total heat flux of natural fumaroles and steam wells added together was estimated to be approximately 15 to 20 MW from 1983 to 2012, and it dropped in recent years [5]. The recent activity of the volcano started around the year 2000, and a small-scale phreatic eruption occurred in Owakudani in 2015 [5]. Hakone Volcano is located approximately 100 km west of Tokyo (Figure 1), the capital of Japan. The region surrounding the volcano has a well-developed infrastructure, comprising railways, roads, hotels, and other recreational amenities. It is a major hot-spring tourism destination in Japan that is visited by approximately 20 million tourists a year.

**Figure 1.** Map showing the area around Hakone Volcano. The topographical map is based on 50 m mesh height data released by the Geospatial Information Authority of Japan (GSI). The red and blue rectangles show the area included in the scenes of the ascending and descending orbits, respectively. Red points show the location of ground control points used for the analysis.

**Figure 2.** Geological map around Hakone Volcano (modified from the Seamless Geological map of Japan [6]). The original geological map was simplified based on ages and types of rocks. Contour lines in which intervals are 100 m were extracted from 10 m mesh height data released by the GSI. The area enclosed by the rectangle indicates the focus area of this study (Figure 3).

**Figure 3.** Location map of the focused area of this study. The base map is a true-color image captured by ALOS/AVNIR-2 on 10 November 2006. Contour lines in which intervals are 50 m were extracted from 10 m mesh height data released by the GSI.

Based on seismic activity and InSAR analyses, the occurrence of the phreatic eruptions of Hakone Volcano and its associated seismic swarm activities are thought to be related to hydrothermal activity in the shallow parts of the volcano [7,8]. It is considered that this activity is closely related to the formation of hot spring water [9]. It is therefore considered important to clarify the reasons for ground deformation in this area, not only for preventing a volcanic disaster, but also to better understand the formation mechanism of hot springs in Hakone Volcano and the conservation of hot spring resources.

In this study, we performed an InSAR time series analysis using the small baseline subset (SBAS) method [10], and data obtained from the Phased Array type L-band Synthetic Aperture Radar (PALSAR) aboard the Advanced Land Observing Satellite (ALOS). We detected highly localized subsidences in the fumarolic area, which is the principal hot spring area of Hakone.

#### **2. Materials and Methods**

We used data from the ALOS/PALSAR satellite launched by the Japan Aerospace Exploration Agency (JAXA), which has a wavelength of 23.6 cm (L-band). The satellite always observes the right side of the azimuth direction, and its revisit interval is 46 days [11]. ALOS/PALSAR was operated from 2006 to 2011, and its successor, ALOS-2/PALSAR-2, has been in operation since 2014. In this study, we obtained 24 scenes from Path 407-Frame 690 observed from the ascending orbit (observation from the western sky) and 22 scenes from Path 59-Frame 2910 observed from the descending orbit (observation from the eastern sky) which covers Hakone Volcano (Figure 1 and Table 1). The off-nadir angles for the scenes of both the ascending and descending orbits were 34.3◦. The ALOS-2/PALSAR-2 data were not used in this study because the data acquisition interval was long, and the data available was insufficient for the analysis.


**Table 1.** Observation dates of scenes used for the interferometric synthetic aperture radar (InSAR) analysis in this study.

<sup>1</sup> Data excluded from the analysis. <sup>2</sup> Data used as the super primary scene.

The ENVI+SARscape software program (Sarmap Sa, Caslano, Switzerland) was used for the SBAS-InSAR time series analysis in this study. Given that coherence is generally worse when the perpendicular baselines (Bperp) between ALOS/PALSAR orbits are more than 1000–1500 m, the analysis was conducted by extracting the interference pairs whose Bperp was less than approximately 1000 m. The pairs with weak interference were then removed visually, and 21 scenes-64 pairs in the ascending orbit and 21 scenes-57 pairs in the descending orbit were used for the analysis (Figure 4). Three of the scenes in the ascending orbit were excluded because of the large Bperp between the other scenes (Table 1). One of the scenes in the descending orbit was excluded because it was observed after the Tohoku-Oki Earthquake (M9.0) of 11 March 2011, and was affected by the earthquake displacement (Table 1).

To improve the signal-to-noise ratio, we averaged over approximately 15 × 15 m for each pixel in the azimuth and range directions. To remove the influence of topography in the initial interferogram, we used 10 m mesh elevation data obtained from the Geospatial Information Authority of Japan (GSI) and the geoid model of the Earth Gravitational Model 2008 [12]. The adaptive interferogram filtering algorithm [13] was applied to reduce the noise, and interferograms were unwrapped using the minimum cost flow approach [14]. During the time series inversion analysis, we used the global navigation satellite system (GNSS) of the GSI for establishing ground control points (GCPs) (Figure 1). The InSAR displacements were corrected with the displacements of the GCPs. For the GCPs, we used daily coordinates of GNSS sites published by the GSI and filtered these coordinate values in a 10-day time window. The velocity distributions obtained from the InSAR time series analysis were converted into a geographic coordinate system with an approximately 25 × 25 m mesh. Then, the velocity maps were obtained by masking pixels with a height estimation error of more than 5 m or a velocity estimation error of more than 8 mm/year. For the descending orbit, the data include Mount Fuji, which is the highest mountain in Japan (3776 m) and forms an independent peak with no continuity with the surrounding mountains. Including Mount Fuji in the scene causes a large discrepancy between the

displacements of GCPs and pixels in the InSAR time series analysis. However, when the western part of the scene, including Mount Fuji, was cropped, the displacements estimated by InSAR and the GCPs almost agreed. Therefore, we employed the result that excluded Mount Fuji.

**Figure 4.** Temporal and spatial baselines for the SBAS-InSAR analysis of ALOS/PALSAR data from the (**a**) ascending orbit, and the (**b**) descending orbit. Red points show the super primary scenes used for the analysis, which the software selected as the scene with the highest number of connections to other scenes.

#### **3. Results**

#### *3.1. Ascending Orbit*

Figure 5 shows the analysis results obtained for the ascending orbit (observations from the western sky). The colors of the figures show the satellite line-of-sight (LOS) velocity; the positive values shown in red indicate the LOS velocity toward the satellite, while the negative values shown in blue indicate the LOS velocity away from the satellite. Since the LOS velocities are corrected by GNSS displacements, these are not relative, but absolute velocities. A LOS velocity of approximately 4 mm/year toward the satellite was observed over the entire area (Figure 5a). A maximum LOS velocity of about 23 mm/year away from the satellite was observed to the west of Owakudani (Figure 5b), the largest fumarole area in the Hakone Volcano.

**Figure 5.** Distribution of line-of-sight (LOS) velocity in the ascending orbit for the (**a**) entire region; (**b**) area of interest shown in the black rectangle (inset) in (**a**). The base map of (**a**) is a topographical map based on 10 m mesh height data released by the GSI, and its color scale is the same as Figure 1. Intervals of contour lines in (**b**) are 50 m in height.

#### *3.2. Descending Orbit*

Figure 6 shows the results of the analysis for the descending orbit (observation from the eastern sky). A LOS velocity of about 4.5–5.8 mm/year away from the satellite was observed over almost the entire area (Figure 6a). In the areas around Hakone Volcano and Mount Fuji, LOS velocities toward the satellite (about 0–7 mm/year) were observed (Figure 6a). A maximum LOS velocity of about 16 mm/year away from the satellite was observed to the west of Owakudani (Figure 6b).

**Figure 6.** Distribution of LOS velocity in the descending orbit for the (**a**) entire region; (**b**) area of interest shown in the black rectangle (inset) in (**a**). The base map of (**a**) is a topographical map based on 10 m mesh height data released by the GSI, and its color scale is the same as Figure 1. Intervals of contour lines in (**b**) are 50 m in height.

#### *3.3. 2.5-D Analysis*

The LOS velocities obtained for ascending and descending orbits were transformed to the velocities of the quasi-eastward and quasi-upward components by applying 2.5-D analysis [15]. Quasi-eastward and quasi-upward velocities (*VQE*, *VQU*) are estimated from LOS velocity in ascending and descending orbits (*VA*, *VD*), and these unit vectors (<sup>→</sup> *A*, → *D*) by following equations.

A normal vector <sup>→</sup> *N* (x, y, z) of the LOS plane created by both LOS vectors are estimated as follows:

$$
\overrightarrow{N} = \overrightarrow{A} \times \overrightarrow{D} \tag{1}
$$

A line vector where the LOS plane and horizontal plane intersect is represented as <sup>→</sup> *L* (y, −x, 0), and its magnitude is as follows:

$$
\left| \overrightarrow{\tilde{L}} \right| = \sqrt{\left(\mathfrak{x}^2 + \mathfrak{y}^2\right)} \tag{2}
$$

Angles (θ*A*, <sup>θ</sup>*D*) between <sup>→</sup> *<sup>L</sup>* and LOS vectors (<sup>→</sup> *A*, → *D*) are estimated as follows:

$$\begin{vmatrix} \stackrel{\rightarrow}{L} \end{vmatrix} \cos \theta\_A = \stackrel{\rightarrow}{A} \stackrel{\rightarrow}{\cdot L} \tag{3}$$

$$\begin{aligned} \left| \overrightarrow{L} \right| \cos \theta\_D &= \overrightarrow{D} \cdot \overrightarrow{L} \end{aligned} \tag{4}$$

Quasi-eastward and quasi-upward velocities (*VQE VQU*) are estimated as follows:

$$\left(\begin{array}{cccc}V\_A & V\_D \\ \end{array}\right) = \left(\begin{array}{cccc}V\_{QE} & V\_{QII}\end{array}\right) \left(\begin{array}{cccc}\cos\theta\_A & \cos\theta\_D\\ \sin\theta\_A & \sin\theta\_D\end{array}\right) \tag{5}$$

The strike and dip of the LOS plane varied slightly from point to point, but these were approximately 269◦ and 85◦, respectively.

The quasi-eastward velocity shows that almost the entire area was displaced westward (Figure 7a). This westward displacement is considered to be a regional tectonic displacement caused by the Pacific Plate, which is subducting beneath the northeastern Japan [16,17]. On the other hand, the areas near Mount Fuji and Hakone Volcano showed eastward velocities (Figure 7a,b). Since expansions of volcanic edifices were observed by GNSS around Mount Fuji and Hakone Volcano during the study period [18,19], the eastward velocities might be consistent with the GNSS displacements.

**Figure 7.** Distribution of quasi-eastward velocity over the (**a**) entire region; (**b**) area of interest shown in the black rectangle (inset) in (**a**). The base map of (**a**) is a topographical map based on 10 m mesh height data released by the GSI, and its color scale is the same as Figure 1. Intervals of contour lines in (**b**) are 50 m in height.

The quasi-upward velocity findings showed that the areas around Mount Fuji, Hakone Volcano, and the Izu Peninsula were uplifted, and the Ashigara Plain, located to the east of Hakone Volcano, subsided slightly (Figure 8a). These displacements are thought to be caused by tectonic deformations by the Izu Peninsula colliding with the main part of the Japanese Islands [3] and volcanic activities [18,19]. To the west of Owakudani, where volcanic rock debris-avalanche deposits are distributed on the surface (Figure 2), a maximum subsidence rate of about 25 mm/year was observed within an area with a diameter of about 500 m (Figure 8b). The displacement appears to have progressed at a constant rate over the observation period (Figure 9). Moreover, regarding the quasi-eastward component, both sides of the center of the subsidence area were displaced inwards; that is, the western side was displaced eastward, and the eastern side was displaced westward (Figure 7b). This result indicates that contractional displacement was dominant to the west of Owakudani.

**Figure 8.** Distribution of quasi-upward velocity over the (**a**) entire region; (**b**) area of interest shown in the black rectangle (inset) in (**a**). The base map of (**a**) is a topographical map based on 10 m mesh height data released by the GSI, and its color scale is the same as Figure 1. Intervals of contour lines in (**b**) are 50 m in height.

**Figure 9.** Time variation of displacement in pixels where the maximum velocities were observed to the west of Owakudani. Displacements were obtained by subtracting the displacement of an arbitrarily selected reference point near Sengokuhara (the green points in Figures 5b and 6b) to remove the displacement around the Hakone Volcano.

#### *3.4. Modeling*

Since the contractional displacements are dominant, we used a model involving a point pressure source in an elastic half-space body [20] to explain the subsidence to the west of Owakudani. A modeling tool in ENVI+SARscape, which implements inversion of the Levenberg–Marquardt least-squares approach [21], was used to estimate the volume variation and 3-D position of the pressure source. During the inversion analysis, the influence of topography was compensated for by the elevation data. Moreover, the offset parameters of the InSAR images were also estimated to consider the LOS velocity around the Hakone Volcano. After the optimum values of the model parameters were obtained by the

inversion, the variations in the residual root mean square (RMS) values were obtained by forward modeling with varying parameters (i.e., altitude and volume variation) to validate the sensitivities of each parameter.

The modeling results show that the subsidence can be explained by a point pressure source contracting at a rate of 1.04 <sup>×</sup> 104 m3/year at the height of about 700 m above sea level (about 300 m below the ground surface) just under the area of subsidence (Table 2 and Figure 10). Based on the residual RMS values (Figure 11) estimated by the forward modeling, the parameters of the model are considered well-constrained, especially with respect to altitude. For example, a residual RMS value of 4 mm or less would be achieved at an altitude of 620–750 m and a volume variation of <sup>−</sup>1.4 <sup>×</sup> 10<sup>4</sup> to <sup>−</sup>6.0 <sup>×</sup> <sup>10</sup><sup>3</sup> m3/year. However, it may be inappropriate to apply a point pressure source model, because the estimated model has a very shallow underground distribution, and it may not be deformed as an elastic body. Modeling that takes into account the physical properties of the ground is a future issue. At a minimum, the modeling results suggest that the cause of subsidence may be due to very shallow contraction, possibly several hundred meters below the ground surface.

**Figure 10.** Results of the inversion analysis for estimating the point pressure source beneath the ground surface to the west of Owakudani: (**a**) observed velocity in the ascending orbit; (**b**) simulated velocity in the ascending orbit; (**c**) observed velocity in the descending orbit; (**d**) simulated velocity in the descending orbit. Contour lines in (**a**,**c**), in which intervals are 25 m, were extracted from 10 m mesh height data released by the GSI. The yellow "x" marks the location of the point pressure source, the parameters of which are shown in Table 2. Estimated offsets for the ascending and descending orbits are +0.15 and +2.43 mm/year, respectively. The color scales used for the velocity are the same as those in Figures 5 and 6.

**Table 2.** Estimated parameters of the point pressure source beneath the area of subsidence to the west of Owakudani.

**Figure 11.** Residual root mean square (RMS) values estimated by forward modeling and varying (**a**) altitude and (**b**) volume variation. Altitude and volume variation changes every 50 m and 2000 m3/year, respectively. During forward modeling, the other parameters were held at optimum values (Table 2).

#### **4. Discussion**

Hot spring wells have been developed in the Owakudani region for decades, and many of these wells are located to the east of the subsidence area. Most of the hot spring wells in Owakudani are steam wells; hot water for bathing is artificially produced by adding local groundwater to the self-discharged volcanic steam that comes up from boreholes, which are several hundred meters in depth. The hot water produced by steam wells in Owakudani is acidic with high concentrations of chloride ions, which are derived from a magma activity of about 10 km in depth [22,23].

On the other hand, in the western region of Owakudani, most of the hot spring water is drawn from the wells. These hot spring waters contain primarily sulfate and bicarbonate, which are naturally produced by volcanic gas added to groundwater distributed near the surface, except for well MO39 (Figures 12 and 13). The hot spring well MO39, which is closest to the area of subsidence, is a steam well that produces acidic spring water from self-discharged volcanic steam with chloride ions as its main dissolved chemical component (Figure 13).

(**b**)

**Figure 12.** Characteristics of the hot spring wells to the west of Owakudani as (**a**) a horizontal distribution of hot spring wells (orange-colored points) on the map of quasi-upward velocity; (**b**) cross-section indicated by A-A' in (**a**). Contour lines in (**a**), in which intervals are 100 m, were extracted from 10 m mesh height data released by the GSI. MO (Motohakone) and SE (Sengokuhara) indicate the numbers of hot spring wells managed by the Odawara Health and Welfare Office. Gray and red colors in (**b**) show the extent of cased and uncased parts of the borehole, respectively. MO4 is a natural hot spring at the surface.

**Figure 13.** Triangular plot of major anions in hot spring water from wells to the west of Owakudani. The numbers in the triangular plot indicate the composition (%) of each ion. The colored circles show the chloride ion content in the hot spring water. The numbers of the hot spring wells and their locations are shown in Figure 12.

Figure 12 shows the depth and characteristics of the hot spring wells to the west of Owakudani. To prevent cold water from mixing with the hot spring water, casing pipes are typically installed from the top of the well to a certain depth in the wells. The casing pipes allow for the collection of the hot spring water or steam directly from the depth of the uncased part of the borehole, where the hot spring aquifers are exposed. Based on the structure of the wells to the west of Owakudani, hot spring water, or steam is collected from around 600–800 m above sea level (Figure 12). These data indicate that the hot spring water or steam in this region is distributed at a depth of several hundred meters below the ground surface. Interestingly, the contraction source estimated by this study is also located at a depth of several hundred meters, which is almost coincident with the hot spring source. Therefore, the subsidence observed in the area is considered to be due to the contraction of the reservoir from which the hydrothermal fluids are being supplied from below.

So, what has caused this reservoir contraction? The Owakudani fumaroles, which is located a few hundred meters east of the area of subsidence, are more active than the fumaroles to the west of Owakudani (Figure 14). In Owakudani, an area 200 m in diameter was locally uplifted during volcanic unrest in 2015, and a phreatic eruption occurred in the vicinity of the uplifted area [7,24]. This local uplifting was explained by expansion sources approximately 80–150 m below the ground surface [7,24]. This expansion suggests the existence of a reservoir that supplies hydrothermal fluids and causes the fumarolic activity on the surface of Owakudani. The hydrothermal fluids with highly concentrated chloride ions are considered to be directly supplied to this reservoir through cracks from the deeper parts of the volcano [7].

**Figure 14.** False-color image around the Owakudani fumarole area captured by ALOS/AVNIR-2 on 31 January 2007. The red color in the image indicates areas of vegetation, and in the Owakudani fumarole area, there is no vegetation due to geothermal and fumarolic activity.

Conversely, there are no such active fumaroles to the west of Owakudani, where the subsidence was observed (Figure 14), although there are small fumarolic activities that make a discharge. However, in 1872, the French engineer Jean-Paul Ishidor Vidal is on record as having measured the temperature of the steam in the western area of Owakudani as 103 ◦C [5,25,26]. This implies that the area of subsidence was once a fumarole zone with heated steam. Moreover, chloride ion that is contained as a major anion in hot spring water of MO39 (Figure 13), which is the closest to the area of subsidence, is implying the participation of the hydrothermal fluids supplied directly from the deeper parts of the volcano [22]. This suggests that the observed subsidence is associated with attenuation of fumarole activity in response to a decrease in the supply of hydrothermal fluids from below. It is also possible that this attenuation may have arisen due to closure of the fluid path from deeper regions, such as a crack. Another possible cause of subsidence may be excessive pumping of hot spring water. However, MO39, the well closest to the subsidence area, is a steam well and uses self-discharged steam, not artificially pumping hot spring water. It is not considered that forced dissipation of pore water pressure, such as land subsidence resulting from excessive groundwater pumping, has occurred. To accurately clarify these factors, it will be necessary to continuously monitor the discharge volume and temperature of the hot spring water.

In this study, we used SAR data acquired by the ALOS/PALSAR satellite from 2006 to 2011. After the termination of ALOS/PALSAR in 2011, the ALOS-2/PALSAR-2 satellite was launched. We consider it essential to continue monitoring the surface displacement in the area to the west of Owakudani to clarify whether the subsidence is continuing or not; this will form part of future studies.

#### **5. Conclusions**

Based on an SBAS-InSAR time series analysis, we observed highly localized subsidence in an area to the west of Owakudani in the Hakone Volcano, Japan, from 2006 to 2011 using ALOS/PALSAR data. The area of subsidence was approximately 500 m in diameter, and the peak rate of subsidence was approximately 25 mm/year. The results of point pressure source modeling suggested that the cause of the observed subsidence was due to reservoir contraction at approximately 700 m above sea level (about 300 m below the ground surface). The rate of contraction was estimated to be 1.04 <sup>×</sup> <sup>10</sup><sup>4</sup> m3/year. Based on the structure of the hot spring wells and their chemical components, the contraction source is considered to be a reservoir containing hydrothermal fluids. Moreover, the subsidence suggests that the attenuation of fumarole activity is occurring due to a decrease in hydrothermal fluid supply.

In addition, the findings of this study demonstrated the application of InSAR surface displacement measurements to monitor hydrothermal activity in shallow parts of volcanic areas. Monitoring surface displacement is considered useful, not only for evaluating the potential occurrence of phreatic eruptions but also for clarifying the production mechanism of hot spring water in volcanic regions.

**Author Contributions:** Conceptualization, R.D.; methodology, R.D.; validation, R.D.; formal analysis, R.D.; investigation, R.D.; data curation, R.D. and G.K.; writing—original draft preparation, R.D.; writing—revision and editing, R.D., G.K. and K.I.; visualization, R.D. and K.I.; funding acquisition, R.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by JSPS KAKENHI, grant number 19K04005.

**Acknowledgments:** We thank Teruyuki Kato for his useful comments and encouragement. We also thank editors and three anonymous reviewers for their constructive comments. ALOS/PALSAR and AVNIR-2 data were provided by JAXA via the Coordinating Committee for the Prediction of Volcanic Eruption as part of the project "ALOS Domestic Demonstration on Disaster Management Application" of the Volcano Working Group. The original ALOS/PALSAR and AVNIR-2 data belong to JAXA/Ministry of Economy, Trade and Industry (METI). The daily coordinates of GNSS sites were obtained through GSI. We used Generic Mapping Tools [27] to draw the figures.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Complex Surface Displacements above the Storage Cavern Field at Epe, NW-Germany, Observed by Multi-Temporal SAR-Interferometry**

#### **Markus Even, Malte Westerhaus and Verena Simon**


Received: 11 September 2020; Accepted: 9 October 2020; Published: 14 October 2020

**Abstract:** The storage cavern field at Epe has been brined out of a salt deposit belonging to the lower Rhine salt flat, which extends under the surface of the North German lowlands and part of the Netherlands. Cavern convergence and operational pressure changes cause surface displacements that have been studied for this work with the help of SAR interferometry (InSAR) using distributed and persistent scatterers. Vertical and East-West movements have been determined based on Sentinel-1 data from ascending and descending orbit. Simple geophysical modeling is used to support InSAR processing and helps to interpret the observations. In particular, an approach is presented that allows to relate the deposit pressures with the observed surface displacements. Seasonal movements occurring over a fen situated over the western part of the storage site further complicate the analysis. Findings are validated with ground truth from levelling and groundwater level measurements.

**Keywords:** satellite geodesy; radar interferometry; persistent scatterers; distributed scatterers; orbit combination; Sentinel-1; cavern field; salt deposit; geophysical modeling

#### **1. Introduction**

The storage cavern field at Epe has been brined out of a salt deposit belonging to the lower Rhine salt flat, which extends under the surface of the North German lowlands and part of the Netherlands. Near Epe the deposit has a thickness between 200 and 400 m and the top of the salt layer (abbreviated as top salt throughout the text) lies in an average depth of 1000 m. The currently 114 caverns are used for brine production and for storage of natural gas, helium and crude oil by in total 8 companies, which follow independent operating strategies. Cavern convergence, i.e., the long-term shrinking of the caverns caused by deformation of the salt, and pressure changes due to the injection and withdrawal of gas cause surface displacements. Mining-caused effects are monitored regularly with levelling, ground water measurements, and other techniques. For our study the potential of SAR interferometry (InSAR) for monitoring nonlinear movements over the storage site has been investigated.

InSAR is a technique that provides valuable information for various monitoring situations and offers its own capacities that complement longer established techniques as levelling and GNSS. Levelling gives high precision measurements of elevation at a moderate number of selected positions with a temporal sampling that often counts in years and GNSS allows to obtain high precision 3D-positioning with dense temporal sampling but with stations usually many kilometers apart from one another. InSAR, on the other hand, provides a high number of displacement measurements in the line of sight (LOS) of the sensor, but the scattering properties of the earth's surface determine the quality of the signal. Nevertheless, usually hundreds of measurements per square kilometer can be utilized. In addition, modern satellite systems allow for a high repeat rate, e.g., six days for Sentinel 1,

thereby providing a high spatio-temporal sampling. With InSAR, phenomena can be studied that cannot be observed with levelling or GNSS alone.

There are several cases where InSAR techniques were successfully used for monitoring of gas reservoirs or storing sites (mostly in combination with geophysical modelling). A particular well investigated example is the in Salah gas field in Algeria [1–7], where natural gas is extracted from a 20-m-thick layer of sandstone via a set of horizontal wells. The natural gas contains excess CO2, which is separated and reinjected at the flanks of the gas field. Geophysical modelling at In Salah making use of the excellent spatio-temporal sampling of InSAR displacement measurements made it possible to retrieve valuable information on the reservoir. It allowed to map reservoir volume changes, to deduce flow properties (e.g., permeability), helped understanding the role of the fault and fracture system in CO2 propagation, and was complementary for the calibration of reservoir models. Similarly, InSAR has been used to calibrate a 3–D fluid-dynamic model and develop a 3–D transversally isotropic geomechanical model for the "Lombardia" gas field in northern Italy [8]. The latter has been successfully used to reproduce the vertical and horizontal cyclic displacements over the storage site. This allows prediction of reservoir behavior under increased pressure (there is great economical interest to increase the working gas volume as much as possible). Another case is the Roswinkel gas field in the northeast of the Netherlands, a severely faulted anticlinal structure, constituting up to 30 reservoir compartments. [9] demonstrates that a carefully executed inversion combining modelling with InSAR can help to identify possibly un-depleted compartments in the reservoir. The authors of [10] have investigated the ability of efficient global optimization to reduce the parameter uncertainties usually affecting geomechanical modeling for the Tengiz giant oil field, Kazakhstan. Efficient global optimization is used to identify the parameter set that minimizes the difference in land displacements obtained from InSAR and 3D geomechanical modeling. In [11] InSAR was used for establishing a risk map for the Solotvyno salt mine area in Ukraine displaying risks of sinkholes and landslides related with mining activities. A recent study [12] investigates the temporal evolution of displacement over the gas storage sites at Lussagnet and Izaute in southwestern France with help of InSAR, data of the Soil Moisture Ocean Salinity (SMOS) satellite and simple modeling. The finding is a linear superposition of several signals: linear, pressure induced, ground temperature induced and soil moisture induced (clay swelling).

For sites with porous storage media as In Salah [13], Berlin-Spandau [14], the "Lombardia" gas field [8], Lussagnet and Izaute [12] or the secret one in the case study [15] the geomechanic response can be described as elastic: displacement is almost proportional to reservoir pressure and displays the same pronounced seasonal behavior. Under different geological circumstances another phenomenon can be observed. [16] investigates cavern integrity at Bryan Mound, where crude petroleum is stored in caverns situated in a salt dome. They observe that pressure peaks in the caverns correspond to peaks of displacement occurring after a time lag of 24 days. The general appearance of the displacement is that of a strongly smoothed and shifted version of the pressure curve. At Epe presumably the same effect is visible. In Section 2.2.1, a temporal model for displacement with pressure changes (pressure response) is derived that relates cavern pressure with observed displacement and is based on the theory for visco-elastic behavior of a Kelvin-Voigt body (for the theory of Kelvin-Voigt bodies see, for instance, [17]). It will be shown that this approach potentially can facilitate InSAR processing and enhance monitoring of cavern storage sites situated in salt deposits.

Further aspects of monitoring with InSAR are density of sampling and determining 2D- or 3D-displacement. Density of sampling depends mainly on the backscattering properties of the earth's surface, acquisition parameters (e.g., sensor, mode, geometry, temporal sampling), and algorithms used for processing. There are two main categories of scatterers that provide usable information for InSAR displacement analysis: persistent scatterers (PS) and distributed scatterers (DS). The former are predominately associated with man-made structures, as most PS are caused by trihedral structures or poles. The latter are often characterized by Gaussian scattering and need to be composed of a sufficiently large number of ground resolution cells that share the same scattering behavior in order to be exploited with help of statistical methods. A prominent algorithm capable of processing jointly PS and DS is SqueeSARTM. SqueeSARTM was used for many of the works cited above [6,7,10–12,16]. The basic idea is pre-processing the DS and using the obtained signals as if they were PS signals in any PS algorithm [6,18,19]. The present work likewise makes use of DS pre-processing [20] and combines it with the software package StaMPS v3. 3b [21,22]. To this purpose, StaMPS was modified in order to allow joint processing of filtered DS and PS and to support unwrapping (Section 2.2.3) with a phase model composed of linear trend, pressure response and a seasonal component that is caused by ground water level changes (Section 2.2.1). Furthermore, in order to avoid leakage of the displacement signal to the spatially correlated noise term, the iterative estimation scheme of StaMPS was refined (Section 2.2.3). With help of these modifications StaMPS processing could be improved successfully for dealing with the challenging displacement field at Epe. Finally, the results presented in this study confirm the validity of the approach to joint processing of PS and DS recently introduced in [20]. The additional information provided by DS allowed to obtain a more complete picture of the displacement field. In particular, it is interesting that it was possible to analyse the displacements over a fen, where a groundwater driven seasonal signal superposed that of the cavern field.

Determining 2D- or 3D-displacements in vertical, east-west and in case, north-south direction from InSAR line of sight-displacements is fundamental for interpretation and integration with other data. As all SAR satellites fly approximately in north-south direction 3D-displacements from InSAR alone can only be obtained near the poles (e.g., at the Henrietta Nesmith Glacier in northern Ellesmere Island, Canada [23]). Hence other strategies have been developed (for a review on resolving three-dimensional surface displacements from InSAR measurements see [24]). Often InSAR data from ascending and descending orbits are combined with assumptions on the physical nature of the deformation. Ref. [25] as well as [26] applied a surface parallel flow assumption to estimate the three-component velocity field for glaciers in Greenland. The authors of [27] based their method on the assumption that landslides in Central–Southern Italy move along the steepest slope. Another possibility is the combination with GPS data. This approach was e.g., applied by [28] to investigate the 1999 Izmit (Turkey) earthquake. In [29], PS-InSAR, GNSS and levelling were combined to estimate small surface displacements in the Upper Rhine Graben. Strong displacements in flight direction can also be estimated from SAR data with speckle tracking [30,31] or using interferograms of sub-apertures [32]. The authors of [33] based their algorithm on the hypothesis that the horizontal displacement is proportional to the tilt using InSAR data of only a single imaging geometry. On the other hand, it is point out in [34] that using two or more imaging geometries provides more robust results and that the combination of ascending and descending orbits is frequently possible. As flight directions deviate roughly 15◦ from north-south any movements in this direction also project to line of sight of the SAR sensor. However, in case of moderate north-south movements and small incidence angles this contribution is relatively small [35]. This is sometimes taken as rationale for determining vertical and east-west displacements using only ascending and descending displacements. With the acquisition geometry of the ascending and descending data of Sentinel-1 used for this work the contribution of the north-south movement to LOS is approximately 15%. Finally, it has to be remarked that in case only LOS displacements from one orbit are available the practice to project them to the vertical assuming wrongly the absence of horizontal displacements does not merely result in false magnitudes but also distorts the spatial pattern of displacement [36].

In this study, a basic method of orbit combination and another one supported by a simplistic geophysical model were applied in order to obtain 2D-displacements. For the basic method the north-south component was handled as if it were zero. The geophysical model predicts the LOS effect of NS displacements. It assumes that caverns act as spherical pressure or volume sources embedded in an elastic halfspace, often called "Mogi" sources in honor of [37], and that the spatial pattern of surface movements results from the superposition of the corresponding displacements. This Multi-Mogi model is used here to describe either the parameters of the linear component of the displacement model or of the pressure response. A novelty of the orbit combination implemented for this study is that the different components of the phase model are combined separately. This allows for a better understanding of the phenomena that contribute to the displacement field.

In Section 2 data and methods are described. Methods of modelling for the support of unwrapping are explained in Section 2.2.1: modelling of the pressure response and of seasonal movements with ground water fluctuations. Section 2.2.2. introduces modelling of the spatial pattern of the linear component of displacement and of the pressure response with help of a Multi-Mogi model. InSAR methods are presented in Section 2.2.3: DS pre-processing and joint processing of DS and PS with a modified StaMPS; support of unwrapping with help of the phase model. Orbit combination is presented in Section 2.2.4. Section 3 presents the results. They are validated versus levelling and ground water measurements. Finally, Sections 4 and 5 give a synopsis and conclusions.

#### **2. Materials and Methods**

#### *2.1. Data*

The SAR data used for this work are from Sentinel-1 in interferometric wide swath mode. Acquisitions from two orbits in ascending and descending flight direction were used. Those from ascending orbit were taken between 3 February 2015 and 7 March 2018 (86 acquisitions, incidence angle 35.21◦), those from descending orbit between 5 February 2015 and 21 March 2018 (118 acquisitions, incidence angle 34.36◦). 5 interferograms from ascending orbit and 3 from descending orbit were discarded because of strong noise.

At Epe 114 caverns have been created by brining out and are used for storage of natural gas, helium, brine and petrol (cp. Figure 1). For each cavern, among other information, northing, easting, volume, salt top, salt bottom, stored medium and operator were provided by the project partner Salzgewinnungsgesellschaft Westfalen (SGW). From the internet site of the Aggregated Gas Storage Inventory (AGSI) [38] historic data with filling levels of gas caverns at Epe are available. In addition, SGW provided levelling data and ground water measurements (GWM). Levelling data were acquired during three measurement campaigns in the years 2015–2017 and comprise 548 points. The positions of measurements are shown in Figure 1. Their colors correspond to linear displacement rates estimated from these measurements. GWM have been taken at 97 locations. Twenty four of them with daily logging are likewise displayed in Figure 1.

**Figure 1.** Locations of levelling measurements (markers colored according to linear trends), ground water measurements (GWM) with daily logging and of caverns. The background map is from Open Street Map.

#### *2.2. Methodology*

Modelling was used for two purposes. First, modelling was used to enhance unwrapping. Significant gradients of the displacement posed a challenge to the unwrapping algorithm of StaMPS and lead to unsatisfactory results. To improve on this, unwrapping was supported by use of a 3D phase model (Section 2.2.1). Second, modelling the spatial displacement pattern helped to obtain more accurate results from orbit combination. To this end, a model comprising multiple Mogi sources in elastic half-space was established in order to describe the spatial pattern of surface movements. With its help, LOS observations of linear convergence and of pressure response from both orbits could be combined resulting in a more accurate estimation of vertical displacement than with a pointwise calculation of vertical and east-west component of displacement. The model is introduced in Section 2.2.2 and its use for orbit combination is explained in Section 2.2.4.

#### 2.2.1. Modelling for Support of Unwrapping

Because of significant temporal gradients of the observed displacement, unwrapping was supported by use of a 3D phase model. The basis for the model is the observation that displacement above the deposit has three main components: 1. linear subsidence caused by convergence of the caverns; 2. pronounced nonlinear movement in response to pressure changes due to the injection and withdrawal of gas; 3. seasonal displacement in peat areas due to water level changes. For each pixel the coefficients of a linear combination of the three components are estimated from the data. For the pressure response and the seasonal movement temporal models were gained from preliminary InSAR results and in case of the pressure response also geophysical considerations. The proceeding is described below. Details on the algorithmic use are given in Section 2.2.3.

In order to define the temporal model for phase variations with pressure changes (pressure response) two steps were performed. For the first step points supposedly displaying a pressure induced phase constituent were selected manually from a preliminary StaMPS result. For both orbits the spatial average (median) over the selected points was calculated, a linear trend was fitted to the median series and subtracted. The detrended median series of the two orbits were arranged to one time series (depicted in Figure 2a) without further scaling. Note that the incidence angles are almost the same such that the vertical displacement projected to line of sight has the same magnitude for both orbits. The east-west component of displacement projects to the lines of sight of the two orbits with opposite signs. Potentially this could cause a split in the combined curve, but there is no clear evidence of such an effect. This is presumably because the projection of the east-west component is small as the spatial median is taken over the area of strongest displacement: horizontal movement of the surface is towards the maximum of subsidence and displacements to the east balance out with those to the west in average. To this combined time series, a model based on cavern pressures is fitted in the second step.

**Figure 2.** (**a**) Model for pressure response with combination of median time series of both orbits selected from preliminary InSAR results and average filling level *fcavern* as in equation (10); (**b**) Model for seasonal displacement together with combination of median time series of both orbits selected from preliminary InSAR results.

Since we did not have detailed information about the cavern pressures of the Epe cavern field, as only averaged pressure data from one operator have been available, we use information from the AGSI [38] on filling levels of gas caverns as a surrogate for pressure. From the AGSI website, filling level time series that reach back sufficiently far are available for two companies and natural gases of types H and L. These comprise 59 of the 76 gas filled caverns. Filling levels are broken down corresponding to company and content and are not available for individual caverns.

A comparison of the time series of filling levels of the caverns with the average InSAR signal derived in the first step shows no simple match. A time shift of approximately 100 days is observed between operational changes of the gas volumes inside the caverns and the supposedly pressure driven variations in the InSAR time series (Figure 2a). We attribute this delay to the visco-elastic behavior of the salt layer wherein the caverns are embedded [39]. To account for this, we employed a simple Kelvin-Voigt-body, that consists of a parallel arrangement of an elastic element with Young's modulus *E* and a viscous element with viscosity η. The constitutive equation for this body (cp. (13.29), p. 690 in [17]) that relates stress σ with strain ε is

$$
\sigma = E \cdot \varepsilon + \eta \cdot \dot{\varepsilon} \tag{1}
$$

Replacing stress σ with the pressure difference *pcavern*(*t*) − *pcavern*(*t*0) and assuming ε(*t*0) = 0 the equation is solved by the delayed pressure function

$$E \cdot \varepsilon = p\_{\rm hop}(t) = \int\_{t\_0}^{t} \dot{p}\_{\rm curv}(s) \cdot \left(1 - e^{-\frac{E}{\eta}(t-s)}\right) ds \tag{2}$$

The subscript indicates that *ptop* is the pressure at salt top from where it is propagated elastically to the surface. To obtain a model curve, it has been assumed

$$p\_{\text{current}} = \alpha\_{1,H} \cdot f\_{1,H} + \alpha\_{1,L} \cdot f\_{1,L} + \alpha\_{2,H} \cdot f\_{2,H} + \alpha\_{2,L} \cdot f\_{2,L} \tag{3}$$

with filling levels *f*∗,<sup>∗</sup> of natural gases of types H and L in the caverns of the two companies for which suitable data have been available as above. The coefficients α∗,<sup>∗</sup> are unknown. Optimizing of retardation time <sup>η</sup> *<sup>E</sup>* and coefficients α∗,<sup>∗</sup> constrained to be positive leads to a model curve *mpres*(*t*) that fits the average signal gained from InSAR reasonably well. *mpres*(*t*) describes the response of the observed phases with respect to variations of the filling levels of the caverns. It differs from *ptop*(*t*) by a factor that scales between pressure and InSAR phases. This factor is estimated in Section 2.2.2 using the physics that are implemented in the Mogi model.

As retardation time 84 days was determined from optimization. Assuming an elastic modulus *<sup>E</sup>* = <sup>25</sup>·10<sup>9</sup> Pa, a viscosity of <sup>η</sup><sup>ˆ</sup> = 1.81·10<sup>17</sup> Pa·<sup>s</sup> results, which both are plausible values for a salt body (see [39], Tables 5.4, 5.8, 5.9 and 5.18 and references given there; The values of *E* given in [39] for different types of salt vary between 10·109 Pa and 37·109 Pa). Figure 2a shows the InSAR signal, the fitted temporal model and the linear combination of the time series of filling levels using the coefficients obtained by optimization. The quality of the fit is good up to the last two months. The reason why it deteriorates during this period is not known to us. A possible cause could be that the assumption that filling levels of a subset of caverns are able to describe the pressures of all gas filled caverns is no longer adequate.

The model for seasonal displacement was as well generated with help of a preliminary StaMPS result. In this case, points displaying significant seasonal phase variations were manually selected. A phase model *a* + *b*·*t* + *c*·*mpres*(*t*) being a linear combination of constant, linear trend and pressure response was fitted to the time series of each of these points and subtracted in order to obtain residual seasonal signals that do not contain contributions of the other components. In order to be insensitive to unwrapping errors and outliers in general, the spatial average (median) over all residual time

series was calculated. That way average seasonal phase time series for both orbits were obtained. Combined, they form a new time series with values for all acquisition times of both orbits. This time series still contains some outliers caused by unwrapping errors in some interferograms that affected whole groups of points. Therefore, robust locally weighted quadratic regression (RLOESS) in a moving window of length 15 acquisitions was applied in order to prevent that the outliers influence the model significantly. This means that inside the moving window a quadratic polynomial is fitted to the time series using robust iterative reweighting. In order to preserve the pronounced shape of the signal the regression was performed separately inside the four intervals defined by the three peaks. Figure 2b shows the InSAR signal and the fitted seasonal model.

It has to be remarked, that we did not use observed time series of groundwater measurements to model the seasonal phase constituent. The reason is that well level variations are often representative for local areas only, and even well levels in immediate vicinity may give significantly differing results. Additionally, a time delay between the response of the surface layers to water input and the water level of a certain well has to be expected. To account for that, diffusion models should have been implemented. We did not do this because we don't have enough information about the governing hydraulic parameters that are supposed to exhibit strong spatial variations. In Section 3.1, we compare INSAR-derived seasonal displacements over a fen in the western part of the cavern field with water level variations of the surrounding wells (cf. Figure 1). Although the levels partly show marked differences to our seasonal model, the general agreement is good which indicates that the model is plausible.

#### 2.2.2. Modelling Used for Orbit Combination

For orbit combination, a simplistic model for the spatial pattern of surface displacements is employed that will be called Multi-Mogi model in this paper. It assumes that each cavern is surrounded by a spherical salt mantel with a thickness of 75 m. The uppermost point of the salt sphere coincides with top salt at the position of a certain cavern. The salt spheres act as pressure or volume sources embedded in an elastic halfspace, often called "Mogi" sources in honor of [37]. The pressure at the outside of the salt spheres is given by *ptop*(*t*), which is related to the pressure variation in the interior of the caverns according to equations (1) and (2). Thus, the model employs an elastic part, governed by the Mogi approach, and a visco-elastic component governed by the Kelvin-Voigt body. We further assume that the surface deformation pattern results from the linear superposition of the contributions of each cavern. For simplicity, we use the term "cavern" to denote the composite pressure source in the following.

The elastic contribution of a cavern situated at (*xc*, *yc*) in depth *dc* at coordinates (*x*, *y*) to surface displacement is proportional to

$$M\_{\mathbf{c}}(\mathbf{x}, \mathbf{y}) = \frac{2\left(1 - \nu^2\right) a\_{\mathbf{c}}^3}{E} \begin{pmatrix} (\mathbf{x} - \mathbf{x}\_{\mathbf{c}})/R^3 \\ (\mathbf{y} - \mathbf{y}\_{\mathbf{c}})/R^3 \\ d\_{\mathbf{c}}/R^3 \end{pmatrix} . \tag{4}$$

where *E* is the elastic module and ν is the Poisson ratio in the elastic halfspace,

$$R = \sqrt{(\mathbf{x} - \mathbf{x}\_c)^2 + (y - y\_c)^2 + d\_c^{-2}} \tag{5}$$

And *dc* = *sc* − *ac*, where *sc* denotes top salt (having negative sign), and *ac* = 75*m* + *rc* is the radius of a sphere of salt with a spherical cavern of radius *rc* in its center (cp. Figure 3). *rc* is the "virtual" radius of the cavern calculated from the cavern volume under the assumption that it is spherical. In reality, caverns are irregularly formed and rather of cylindrical or ellipsoidal shape. The assumption that the top of the spherical cavern is located 75 m below top salt is somewhat arbitrary, as we don't know the exact vertical positions of the caverns in the salt layer. Positions, values of top salt and

volumes of the caverns were provided by the project partners. For the elastic half space an elastic modulus *<sup>E</sup>* = 30·109 Pa and a Poisson ration <sup>ν</sup> = 0.25 have been chosen, because values for the cap rock at Epe are not known to us. As a consequence of the big diversity of values of *E* for cap rock, this number is even less reliable as that assumed for the salt. As a final remark, it should be pointed out that the uncertainty regarding the assumed values of Young's modulus *E* has no influence on the results from InSAR as they always occur as *E* times estimated parameter, where the estimated parameter is not interpreted on its own. The Mogi model is used here to describe two components of the displacement field: either for the parameters of the linear component of the displacement model or of the pressure response. In the case of the linear component it is assumed that all caverns converge with the same rate. This means only one parameter *plin* is assumed to be unknown (describing volume change with time *t*) and has to be determined such that the projections to LOS of the displacement

$$D\_{\rm lin}(\mathbf{x}, y, t) = t \cdot p\_{\rm lin^\*} \sum\_{c \in \mathcal{C}} M\_{\rm c}(\mathbf{x}, y) \tag{6}$$

approximate optimally the linear model component for both orbits and all points found by InSAR processing (*C* denotes the set of all caverns). According to the Mogi [37] formalism, *plin* is related to a continuous volume loss *Vlin* of the spherical pressure source:

$$p\_{\rm lin} = \frac{V\_{\rm lin} \cdot E}{a\_c^3 \cdot \pi \cdot 2 (1 + \nu)}\tag{7}$$

**Figure 3.** Sketch of the geometric arrangement of a cavern as used for modelling (*dc* depth of cavern center, *sc* top salt, *ac* radius of salt sphere, *rc* cavern radius).

Optimization gives annual volumetric convergence rates of the order of 1%, which are quite realistic. In case of the pressure response, only gas containing caverns are assumed to contribute and displacement is proportional to the pressure on the surface of the salt sphere, which is obtained from cavern pressure according to Equation (2). The pressure in all these caverns is assumed to be the same at all times. Displacement at coordinates (*x*, *y*) at time *t* as described by the Multi-Mogi model is obtained as

$$D\_{\rm pres}(\mathbf{x}, \mathbf{y}, t) = p\_{\rm top}(t) \cdot \sum\_{\mathbf{c} \in \mathcal{C}(\rm pres)} M\_{\mathbf{c}}(\mathbf{x}, \mathbf{y}) = m\_{\rm pres}(t) \cdot p\_{\rm pres} \cdot \sum\_{\mathbf{c} \in \mathcal{C}(\rm pres)} M\_{\mathbf{c}}(\mathbf{x}, \mathbf{y}) \tag{8}$$

As above *mpres* is the model for the pressure response and *C*(*pres*) denotes the set of all gas filled caverns. The parameter *ppres* is determined such that the projections to LOS of *Dpres* approximate optimally the component of the pressure response for both orbits and all points found by InSAR processing.

The parameter *ppres* estimated for the pressure response allows to guess the difference Δ*Pmax* of the maximal cavern pressure to the pressure *P*<sup>0</sup> corresponding to filling level zero (unknown to us). To this end we calculate the pressure difference

$$
\Delta p\_{\rm top} = p\_{\rm top}(t\_{\rm max}) - p\_{\rm top}(t\_{\rm min}) = p\_{\rm pre} \cdot \left( m\_{\rm pre}(t\_{\rm max}) - m\_{\rm pre}(t\_{\rm min}) \right) \tag{9}
$$

effecting elastically the maximal displacement [37], i.e., the difference between maximum and minimum of the displacement signal. The cavern pressure is proportional to the amount of gas in the cavern according to the Van der Waals equation for constant temperature. Temperatures in the depth of the caverns can be assumed to be constant on average. If it is assumed that the weighted average *fcavern* of known filling levels (with sum of weights equal one) is a good approximation of the actual filling levels, then the cavern pressure can be expressed as

$$p\_{\text{caverv}} = \Delta P\_{\text{max}} \cdot f\_{\text{caverv}} + P\_0 = \Delta P\_{\text{max}} \cdot \frac{\pounds\_{1,H} \cdot f\_{1,H} + \pounds\_{1,L} \cdot f\_{1,L} + \pounds\_{2,H} \cdot f\_{2,H} + \pounds\_{2,L} \cdot f\_{2,L}}{\pounds\_{1,H} + \pounds\_{1,L} + \pounds\_{2,H} + \pounds\_{2,L}} + P\_0 \tag{10}$$

and the left-hand side of Equation (9) can be calculated by evaluating Equation (2) at *tmax* and *tmin*:

$$\Delta p\_{ttop} = \Delta P\_{max} \cdot \left\{ \int\_{t\_0}^{t\_{max}} f\_{carrun}(s) \cdot \left(1 - e^{-\frac{E}{\eta}(t-s)}\right) ds - \int\_{t\_0}^{t\_{min}} f\_{carrun}(s) \cdot \left(1 - e^{-\frac{E}{\eta}(t-s)}\right) ds \right\} \tag{11}$$

From (9) and (11) the guess of <sup>Δ</sup>*Pmax* can be calculated. Anticipating the result, <sup>Δ</sup>*Pmax* = 5.3·106 Pa is obtained, which probably is at the lower end of the real pressure values. Nevertheless, the order of magnitude is reasonable which indicates that the simplistic model is physically plausible.

#### 2.2.3. InSAR Methodology

The InSAR processing runs through three principal steps: 1. Coregistration and interferogram formation with ESAs Sentinel-1 Toolbox; 2. DS pre-processing following the paradigm of SqueeSAR; 3. Joint processing of DS and PS with a modified and augmented version of StaMPS 3.3b1. Coregistration and interferogram formation are standard Sentinel Application Platform (SNAP) functionality and are not further discussed here.

In order to jointly process DS and PS with StaMPS two principal algorithmic changes have been performed: 1. The introduction of DS pre-processing as a new component, which estimates for each pixel a complex signal that in case of good quality can be used like that of a PS in an arbitrary PS algorithm. 2. The modification of the pixel selection step of StaMPS, where is decided which points are kept for further processing. This comprises assessing the quality of DS and PS signals and also deciding if the DS or PS signal shall be used in case both are of good quality.

The basic concept for pre-processing of DS stems from the SqueeSAR paper [18] and consists of grouping for each pixel a statistical homogeneous neighborhood, estimating the covariance matrix and DS signal estimation. Since the introduction of SqueeSAR, a multitude of approaches has been developed to perform these tasks e.g., [40–46], for a review see [19]. For the results presented in this work the following methods were applied: The approach of [44] was used for the grouping step. It consists of signal transformation, outlier removal based on the adjusted boxplot, and application of the one-sample t-test. The one-sample t-test was applied with significance level 0.99. From each pixel in the neighborhood up to eight outliers were allowed to be removed (the investigated stacks comprise 86 and 118 acquisitions). Pixels with more outliers were discarded. The covariance was estimated as a

*A*

3D-version of the sample covariance matrix (outlying values are set to zero). The square roots of its diagonal values serve as amplitudes of the DS signal. The phases of the DS signal were estimated with phase triangulation coherence maximization [19,47], which is for small neighborhoods more reliable than the maximum likelihood estimator derived in [48].

A distinctive feature of the approach used for this work is that the decision between alternative use of DS or PS signal is made at an early stage before the pixel selection step of StaMPS. The motivation for this proceeding derives from the design of the pixel selection step: those pixels are retained that possess small residuals compared to a background phase obtained by spatial filtering. As the reliability of the background phase depends on the availability of high quality pixels, it is plausible to combine DS and PS right from the start. The drawback is that this means to abstain from the use of the residual phases as a characteristic of quality for deciding between the use of DS or PS signals. Therefore, in order to make this decision, a newly developed criterion based on the difference between amplitude dispersion of the DS signal and of the PS signal is applied [20]. The approach of spatially filtering DS and PS separately and then deciding based on the residuals has been employed by [40]. A comparison between both approaches has not yet been done. Figure 4a gives a flowchart of the approach used for this work. The characteristics used are amplitude dispersion *DPS <sup>A</sup>* for the PS signal, amplitude dispersion *DDS <sup>A</sup>* for the DS signal, the difference <sup>Δ</sup>*DADPS <sup>A</sup>* <sup>−</sup> *DDS <sup>A</sup>* of amplitude dispersions, the number of pixels #Ω of the neighborhood of the DS, phase triangulation coherence γ*pt* and the temporal coherence <sup>γ</sup>*sel* calculated from the residuals with respect to the background phase. <sup>T</sup><sup>∗</sup> denotes thresholds for the different characteristics (for more details see [20]). For this study the settings were T*DPS A* = 0.45, T*DDS* = 0.45, TΔ*DA* = 0.0, T#<sup>Ω</sup> = 20, Tγ*pt* = 0.45, and Tγ*sel* = 0.8.

**Figure 4.** Flowcharts for joint processing of PS and DS and for orbit combination. Symbols are explained in the text.

Furthermore, the iteratively repeated steps of phase unwrapping, estimation of spatially correlated look angle error (scla), and estimation of spatially correlated noise (scn) of StaMPS had to be modified because of the strong gradients of the deformation field and leakage of signal into the scn term. These changes concern the arrangement of corrections for already estimated error terms done in the algorithmic steps depending on the iteration step (Table 1 gives an overview). In particular, as the most significant modification the estimation of a phase model that is either estimated during the unwrapping step from wrapped phase (this is explained in the next paragraph) or during the scn estimation step from unwrapped phase (indicated in Table 1 by w or uw) is newly introduced. In Table 1 the last four columns correspond to terms that are estimated during processing. The spatially uncorrelated look angle error (sucla) has been estimated earlier after joint spatial filtering from the wrapped phase. The other three are estimated during iterations. An entry in any of these columns indicates that the phase is corrected for this term before the algorithmic step is performed. The particularities of the scheme are based on the following considerations:



**Table 1.** Overview of corrections done during iterative estimation.

<sup>1</sup> spatially uncorrelated look angle error, <sup>2</sup> spatially correlated look angle error, <sup>3</sup> spatially correlated noise.

As said above, the most significant modification is the introduction of the model. During iteration 2 it is estimated from the wrapped phase. Before unwrapping, already known error terms from the first iteration are subtracted from phase. Between this preparative step and the call of the unwrapping algorithm, this estimation step has been newly inserted. It performs a parameter space search for a phase model. Those parameters are selected which maximize the temporal coherence of the residuals. The phase model for the gas deposit assumes a linear combination of linear trend, pressure response and seasonal movement and is subtracted from the wrapped phase before unwrapping starts. After the unwrapping algorithm has terminated, the modeled phase is added to its result.

In order to make this step reliable several improvements have been added:


#### 2.2.4. Orbit Combination

As DS and PS lie scattered and do not coincide for different orbits the two outputs of StaMPS have to be interpolated to common positions for the purpose of orbit combination. Likewise, times of acquisition do not agree, which makes temporal interpolation necessary. Unfortunately, interpolation and filtering tend to attenuate the estimated signal, particularly in the presence of noise. In this situation, the a priori knowledge about the deformation process that is reflected in the Multi-Mogi model described in Section 2.2.2 might potentially be used to preserve better the underlying signal. Hence orbit combination was implemented in two ways that are explained in the sequel. Before doing so, the principal steps of orbit combination are described:


This basic version of orbit combination uses the displacement time series from StaMPS. In order to generate the results in the next section, a more complex approach has been adopted (Figure 4b gives a flowchart). A parametric linear model consisting of constant, linear trend, pressure response and seasonal movement was fitted to the displacement time series. This way they were split in a modeled part and an unmodeled residual. Steps 3 to 5 were modified accordingly: (a) for the unmodeled residual they have been performed as before; (b) for the modeled part kriging filtering has been applied to the parameters and the transformation from LOS to vertical and east-west component has been performed separately for them, which has the advantage to preserve the decomposition in the different components thus allowing to investigate them separately. The algorithm used for kriging filtering of the parameters is different for two of them. For all parameters the assumption of a constant drift is clearly not valid. Therefore, the use of ordinary kriging is not recommendable, while using universal kriging is a potential alternative. To apply universal kriging a model for the drift is required. We do not have one for all parameters. The first parameter is not relevant for displacement analysis and appears to represent not removed master atmosphere. Hence it does not need extra care and ordinary kriging has been applied. The second and third parameter are related to the geophysical behavior of the deposit. For these two parameters models for the drift have been constructed with help of the Multi-Mogi model and universal kriging has been performed. The fourth parameter has been filtered preliminary with ordinary kriging despite the named objection, as a model for the drift has not yet been devised. We will refer to this method as "geometric" for brevity.

Alternative to universal kriging filtering for the second and third parameter and transforming the result to vertical and horizontal coordinates a second approach was tested. For each of the two displacement parameters the parameters of the Multi-Mogi model were estimated by optimizing for both orbits at the same time the fit to the LOS displacement parameter values at their original, orbit dependent positions (cp. Section 2.2.1). In case of the linear component the value of *plin* is determined that minimizes


Here *Rasc* and *Rdesc* are the sets of pixels from the two InSAR results for which the estimated model parameters *pasc lin* (*x*, *<sup>y</sup>*) or *<sup>p</sup>desc lin* (*x*, *<sup>y</sup>*) are bigger than some threshold. *<sup>A</sup>asc*(*x*, *<sup>y</sup>*) and *<sup>A</sup>desc*(*x*, *<sup>y</sup>*) are the vectors that project displacements according to the Multi-Mogi model to line of sight of each of the orbits. The value of *ppres* is determined analogously. With the obtained parameters of the Multi-Mogi model the displacement parameters in vertical and horizontal direction at arbitrary positions can be calculated. We will refer to this method as "Multi-Mogi".

Furthermore, it has to be remarked that the procedure of estimating the parameters of the Multi-Mogi model differed for the two parameters. *plin* describes the subsidence caused by convergence of the caverns. Therefore, all caverns contribute to the deformation, although not necessarily to the same degree. On the other hand, *ppres* describes the magnitude of the response to pressure changes in gas filled caverns. It is assumed that caverns filled with liquids do not contribute to this component of the displacement signal and consequentially these caverns are not included in the model in this case. Finally, points with big value of the parameter of the seasonal component have been masked out in order to prevent any adverse effects on the estimation because of correlation of the model for pressure induced displacement with the seasonal model.

#### **3. Results**

Data from both orbits have been processed with DS-preprocessing and the augmented StaMPS algorithm. A four dimensional linear phase model has been fitted to each time series of the results.

#### *3.1. Joint Processing of DS and PS*

Figure 5 displays the model parameters estimated for both orbits. On the left hand side results from ascending orbit are shown and on the right hand side those of descending orbit. From top to bottom linear trend, pressure response and seasonal signal are given. The color code for the linear trend expresses mm/y and negative values indicate subsidence. For the other two parameters the maximal displacement in mm is depicted, i.e., the difference between maximum and minimum of the signal. Significant linear trends can be observed over the whole cavern field. Apparently all caverns suffer from convergence and contribute to the signal. The other two components of the displacement signal express differently in the western part and in the eastern part. In the western part, the majority of liquid filled caverns is situated and a fen can be found that is displayed by a brownish color in the background map (cp. Figure 1). In the eastern part, gas filled caverns dominate (cp. Figure 1). As one

would expect, surface movements modeled by the pressure response occur over the eastern part of the cavern field, where gas is stored. Seasonal movements, on the other hand, are restricted to a fen in the western part of the cavern field. We demonstrate below that the seasonal movements are presumably caused by groundwater level changes. The presence of horizontal movements can be recognized for linear trends and, less distinctly, for displacements with pressure changes, e.g., the areas of strongest displacement appear to be different. The points with significant seasonal signal in the fen are all DS and no information was available in this area from a reference result generated with StaMPS using only PS. It is noteworthy that these points extend the usable information on convergence and pressure response to this area (cp. also Figures 12 and 13). The related case of using InSAR to determine displacements of pasture on drained peat soils was demonstrated to be a challenging task [51,52], but the authors predicted that the 6-day repeat pass of Sentinel-1 might improve the situation. The use of Sentinel-1 and the presence of structures caused by peat cutting might have contributed to a successful evaluation of the DS signals in this work. Figure 6 displays as indicator of pixel quality the root mean square (RMS) of the residual phases for each point for both orbits. Residual phase means here phase of the final result minus modelled phase. The residual phase comprises non-compensated errors (e.g., white noise) and the not modelled signal. The latter presumably gives a partial explanation for the high values of RMS over the fen, because the actual highly irregular seasonal signal supposedly deviates from its model.

Figure 7 displays two time series of levels of groundwater measurements (GWM) taken in the fen area together with the mean detrended displacement signals of nearby DS of both orbits (the InSAR signals have been scaled to match the groundwater levels). A good correlation between groundwater levels and surface movements detected by InSAR is evident. However, a small time delay between InSAR and GWM is visible that might be explained by the time the water needs to seep to the location of GWM after rain.

#### *3.2. Orbit Combination*

In this section results obtained after orbit combination are presented. As explained in Section 2.2.4 orbit combination was performed by two different methods to which is referred as geometric or Multi-Mogi in the following. In order to assess the quality of the different approaches orbit combinations have been done at those coordinates, where levelling data are available. This is different from the purpose of obtaining displacement fields as presented below, where the combination is done for those cells of a raster that contains information from both orbits. As a characteristic of quality the standard deviation σ between InSAR time series and levelling has been calculated at each coordinate. To do so, some assumptions were made. Levelling measurements were performed each year between middle of May and middle of July. It is not possible to assign a precise point of time to measurements. Hence σ was calculated by interpolating the InSAR displacement time series to the 15th of June of each year. Furthermore, data from some of the levelling locations were not used. These either comprise only two values or are not consistent with nearby measurements. For the latter a comparison with InSAR results does not make sense as this requires interpolation to the positions of the levelling points and cannot capture local peculiarities. This way 517 of 548 levelling points were retained and used for the comparison. Figure 8 shows the color coded values of σ for the two investigated approaches. The results on the left-hand side use geometric orbit combination, those on the right-hand side the Multi-Mogi approach. Both comparisons show deviations between InSAR result and levelling at two regions, where the approximately north-south running levelling line leaves the subsidence bowl. Deviations for those using the Multi-Mogi approach are less severe. Two points with large deviation between the geometric approach and levelling stand out, while Multi-Mogi is not sensitive to estimation errors affecting single points by design and provides good results. This is confirmed by the histograms of differences between InSAR and levelling in Figure 9. The distinct negative bias observed for the geometric approach is a bit surprising as it means that the deformation measured by InSAR is stronger than that measured by levelling. Partially, it can be explained by the spatial arrangement of PS and Ds on the one hand and levelling points on the other. The larger portion of the points showing significant

differences are situated along the road south of the cavern field. North of the road is the area of strongest subsidence, which is well covered by InSAR, while south of the road an area without InSAR points abuts. Hence, kriging uses predominately points with subsidence rates higher than measured by levelling to interpolate InSAR results to levelling positions. Another source of error is the ignorance of the north component of displacement for the geometric approach. To get an idea of the order of its consequences this error has been calculated based on the north component estimated with Multi-Mogi (Figure 10). It has 0.08 mm/y bias and 1.3 mm/y standard deviation. From this it appears that it does not significantly influence the bias and merely increases the standard deviation notably.

**Figure 5.** Estimated parameters in line of sight. On the left from ascending orbit, on the right from descending orbit. From first row to third row: linear rate, pressure response, seasonal displacement. The brownish color surrounding the points with strongest seasonal displacement displays the fen.

**Figure 6.** Root Mean Square (RMS) of residual phase. Left-hand side from ascending orbit, right-hand side from descending orbit.

**Figure 7.** Examples of ground water levels together with scaled mean over InSAR time series of nearby points from each orbit.

**Figure 8.** Standard deviations of differences between InSAR results and levelling. Left-hand side from geometric orbit combination, right-hand side Multi-Mogi approach.

Figure 11 shows the vertical displacement for some time series for each of the approaches together with the model curves estimated and levelling. The model curves are a combination of three components: linear subsidence, pressure induced and seasonal displacement. The given value of σ is calculated for Multi-Mogi and is intended to convey a visual impression of the significance of the deviation that goes along with the value of σ. Figure 11d is an example where the non-modeled signal contains valuable information at locations, where the model does not fit perfectly to the time series.

There are also cases where a larger deviation possibly is due to inaccuracies of levelling. Figure 11e shows a suspected case.

**Figure 9.** Histograms and scatter plots of differences between InSAR results and levelling. Left-hand side from geometric orbit combination, right-hand side Multi-Mogi approach.

**Figure 10.** Error of vertical displacement because of ignorance of north-south component.

Figure 12 shows the linear vertical displacement rate estimated from the time series generated with the Multi-Mogi approach, where data have been interpolated to a raster. Linear rates for levelling points are displayed as crosses in the same color scale. Both data sets blend nicely, what demonstrates again the good agreement between results from Multi-Mogi and levelling measurements. It is clearly seen that the DS located at the fen carry relevant information about the convergence process of the caverns and thereby contribute to a spatially dense displacement field that could not have been obtained from PS alone. However, Multi-Mogi does not capture as well the horizontal displacement field as it does the vertical as indicates a comparison with the geometric approach. The first row of Figure 13 shows linear displacement rates in vertical and east-west direction estimated from the time series generated with geometric orbit combination. The second row shows the corresponding results generated with the Multi-Mogi approach. The horizontal displacement extends far into the surroundings of the deposit and appears blurred in comparison with the result of the geometric combination.

**Figure 11.** Examples of vertical displacement time series for selected points.

**Figure 12.** Linear displacement rate transformed to vertical with Multi-Mogi approach (bullets) together with linear rates from levelling (crosses).

**Figure 13.** Linear displacement rates transformed to vertical and east-west direction. First row from geometric orbit combination, second row with Multi-Mogi approach.

Finally, Figure 14 displays for geometric orbit combination the vertical and east-west component of the parameters for the pressure response and the seasonal signal. The pressure signal which is related to deeply situated sources, has a large spatial dimension in both components. In contrast to this, the seasonal effect is restricted to the vertical component in the immediate vicinity of the fen, which indicates that the water-driven deformation takes place only in the uppermost soil layers.

**Figure 14.** Parameters for pressure response and seasonal displacement transformed to vertical and east-west direction from geometric orbit combination.

#### **4. Discussion**

The results presented in Section 3 demonstrate that multi-temporal InSAR techniques using Sentinel-1 data allow obtaining valuable information about complex surface displacement fields as that of the underground storage at Epe. The InSAR derived displacements are in reasonable agreement with levelling data for geometric orbit combination and in good agreement for the Multi-Mogi approach. They complement levelling with information from additional locations and deliver horizontal movements. In particular, a high temporal sampling is possible that is difficult to achieve with levelling. The successful analysis has been made possible by combining StaMPS with DS pre-processing and introducing a displacement model that improved unwrapping as well as orbit combination. DS allow to extract information also in areas, where no PS are available, in particular over the fen in the western part of the cavern field. In this way a more complete description of the displacement field is obtained.

The displacement model, in addition to its benefits for processing, helps in interpreting the observations by decomposing the displacement signal in physically founded constituents. Although direct access to absolute pressure data from each cavern was not given, the high correlations of InSAR time series with the modeled pressure response based on filling levels of AGSI as a pressure surrogate and on assuming propagation through a visco-elastic medium, suggests that the proposed causal relationship exists and operational pressure changes in the caverns are accompanied by surface

displacements of up to ± 6 mm. The presented model driven approach facilitates a clear separation of the pressure related displacements from ground water induced surface movements.

The InSAR based results also exhibit some shortcomings. The results of the geometric approach for the vertical and east-west component of displacement suffer from ignorance of the north-south component but still are reasonable. The observed bias between InSAR and levelling can be explained partially by the required interpolation. On the other hand, the Multi-Mogi approach captures vertical displacement well. However, it does overestimate the magnitude of horizontal displacement. This shortcoming presumably can be accounted for by the following two arguments:


Future research might help to resolve these shortcomings. A better modeling of the geophysical realities potentially can provide good estimation of vertical as well as of horizontal displacement. Finally, repeating this study with data from a longer period would serve to validate the approach presented here.

#### **5. Conclusions**

In this study the complex displacement field over the cavern storage site at Epe, NW-Germany, was investigated. The temporal behavior of each point could be modeled as linear combination of three components: linear subsidence caused by cavern convergence, displacement responding to pressure changes in the gas filled caverns and a seasonal movement over a fen that shows a good correlation with groundwater measurements. In order to achieve the presented results, simplistic geomechanical modelling of the temporal signal, improvements of InSAR processing and component-wise orbit combination were introduced.


In addition, a method that allows to perform an orbit combination based on simplistic geomechanical modelling of the spatial displacement field was presented. It treats caverns as Mogi sources, whose volume changes because of convergence lead to linear subsidence and whose pressure changes result in the pressure induced displacement. The model assumes for each of these two components one unknown parameter that is estimated such that the projections to line of sight for both orbits optimally approximate the orbit-wise estimations. We named this approach Multi Mogi approach.

The vertical components of the estimated displacements were validated with data from 517 levelling points available over the storage site. The agreement of linear trends of vertical displacement from geometric orbit combination and levelling was reasonable (root mean square 3.41 mm/y). The agreement of linear trends of vertical displacement from the Multi Mogi approach and levelling was very good (root mean square 2.39 mm/y). Unfortunately, the Multi Mogi approach overestimated the east-west component of displacement. This may be due to insufficient knowledge of physical parameters or imperfections of the model itself and deserves further investigation.

**Author Contributions:** Conceptualization, M.W.; methodology, M.E. and M.W.; software, M.E.; formal analysis, M.E.; preliminary analysis, V.S.; data curation, M.E. and V.S.; writing—original draft preparation, M.E.; writing—review and editing, M.W; visualization, M.E; project administration, M.W; funding acquisition, M.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** The present study was financed by BMBF via the project SUBI. We thank BMBF for making our work possible.

**Acknowledgments:** We are grateful for the support and data we have received from our project partners from SGW, in particular Stefan Meyer and Bernd Feldhaus, and from Uniper, in particular Bela Szöcs, Katrin Vosbeck, and Tobias Rudolph. We thank Andy Hooper for providing StaMPS as open source and ESA for the Sentinel-1 data and SNAP (Sentinel Application Platform).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Multidisciplinary Analysis of Ground Movements: An Underground Gas Storage Case Study**

#### **Christoforos Benetatos, Giulia Codegone, Carmela Ferraro, Andrea Mantegazzi, Vera Rocca, Giorgio Tango and Francesco Trillo**


Received: 23 September 2020; Accepted: 20 October 2020; Published: 23 October 2020

**Abstract:** The paper presents a multi-physics investigation of the ground movements related to the cyclical and seasonal injection and withdrawal of natural gas in/from a depleted reservoir located in the Po Plain area, Italy. Interferometric Synthetic Aperture Radar (InSAr) data (from 2003) and Global Navigation Satellite System (GNSS) data (from 2008) provided a full and coherent panorama of almost two decades of ground movement in the monitored area (more extended than the field boundary). The analysis of the acquired millimetric-scale movements together with the detailed geological analysis, both at reservoir and at regional scale, represents the focal point for understanding the investigated phenomena. Based on this information, a fully integrated and multidisciplinary geological, fluid-flow and geomechanical numerical modeling approach was developed to reproduce the main geometrical and structural features of the involved formations together with the poromechanics processes induced by the storage operations. The main achievement of the adopted methodology is a deep knowledge of the system and the involved processes, which is mandatory for the safety of the urbanized areas and the effective management of the underground resources.

**Keywords:** underground gas storage; InSAR; GNSS; ground movement; subsidence monitoring; integrated numerical simulation

#### **1. Introduction**

The ground surface movements due to anthropogenic activities, such as underground fluid production/injection, have been deeply investigated by the scientific literature [1–4], among the others, particularly in potentially critical environments such as high-urbanized areas and costal zones. In the Italian panorama, a large number of hydrocarbon fields were discovered and produced since the early 1950s onward in the Po Plain area (e.g., [5]); due to its high degree of urbanization, the area has been the focus of numerous studies [6–11], among others. The monitoring, investigation and forecasting approaches adopted for the subsidence analyses have shown a constant and continuous improvement of the technologies and methodologies to fulfill the higher and higher standards of the system safety and the social acceptance. The ground movement monitoring has improved from Global Navigation Satellite System (GNSS) to Interferometric Synthetic Aperture Radar (InSAR) technologies and the analysis of the involved multiple physical processes has improved from analytical to 3D numerical approaches.

The investigation of ground movements induced by the underground storage of natural gas (UGS) represents a specific branch of the Oil and Gas industry. During UGS activity, natural gas is

stored within subsurface geological formations (i.e., depleted hydrocarbon reservoirs, aquifers, salt caverns). The UGS is a worldwide solution adopted to guarantee a real-time response to the market gas requests, a high degree of elasticity in the management of production and transport structures, and the maintenance of "strategic" reserves. Seasonal and cyclical withdrawal and injection of gas induce an alike seasonal and cyclical oscillation (subsidence/rebound) of the ground surface, the so-called "earth breathing" phenomenon [12]. The magnitude and the extension of the phenomenon together with its time of occurrence, depend on a large number of factors including: the withdrawal/injection fluid volume and the induced pressure variation, the presence of surrounding aquifers, the depth, shape and volume of the reservoir formation (and its eventual aquifer) together with its petrophysical parameters, and the geomechanical properties of the reservoir and its surrounding formations [13]. The induced ground oscillation is assimilable to a sinusoidal "signal" which can be discriminated from the background noise; the latter includes the effects of ground water production, thermic and meteoric phenomena, sediment compression due to urban settlements, naturally occurring geological processes, among others. Therefore, the quality and the accuracy of the ground monitoring information can be extremely high if suitable acquisition/interpretation techniques are adopted.

The paper presents a case study concerning the ground movements induced by an UGS field operated in a depleted gas reservoir in the Po Plain area (Northern Italy). A reliable ground movement-monitoring plan via both InSAR and GNSS technologies, and a detailed geological research, at both reservoir and regional scale, provided the key information for the set up and calibration of the geological, fluid-flow and geomechanical numerical models using the Schlumberger Petrel™ software suite. The adopted multidisciplinary simulation approach deeply investigated and reproduced the main geological and structural features of the system together with the fluid-flow and poro-mechanical process related to UGS activities; furthermore, the achieved deep knowledge of the system and the involved processes guarantees the safety of local infrastructures and surrounding environment.

#### **2. Materials and Methods**

The aim of the surface monitoring plan was the determination of the planar and vertical components of movement and of the uplift/subsidence area ascribed to UGS operations via two complementary technologies: the satellite SAR interferometry and the GNSS techniques. SAR interferometry investigated an area larger than the field area providing relative measurements of movements whereas GNSS monitoring provided punctual and absolute measurements.

The SAR interferometry is widely used to detect and monitor slow terrain movements with millimeter accuracy along the satellite line of sight (LOS), covering long-term period with frequent updates and large areas with tens of thousands of measurements per square kilometer. The SAR image phase is dependent both on the distance between the radar antenna and the ground target, and on the radiometric characteristics of the transmission medium and of the target. From a stack of interferometric phases (i.e., phase difference between two SAR acquisitions), displacement information of the observed scene is extracted [14] for a set of sparse points (Persistent Scatterers—PSs). The PSs are physical targets characterized by temporally stable backscattering properties and they are commonly found in poorly vegetated and urbanized areas [15–17]. The present research was developed according to the Persistent Scatterer Pairs (PSP-IFSAR) approach; it is a multi-interferogram InSAR technique for deriving dense and reliable information in a set of sparse points (PSs). The adopted PSP-IFSAR method, which integrates both ascending and descending acquisition geometries, measures the same ground deformation along two-satellite lines of sight, and, consequently, the actual vertical and East–West planar components of the movements can be determined. The North–South component cannot be derived because the SAR sensors are insensitive to the North–South displacements.

The GNSS (formerly "GPS") is a space geodetic technique able to determine 3D coordinates of a monitoring permanent site with sub-centimeter accuracy and its over-time kinematic evolution (i.e., its velocity of movement along the planar and vertical components). This geolocation technique is based on more than one constellation of satellites (NAVSTAR for the GPS), numerous ground permanent stations receiving satellite signals and an International GNSS Service (IGS), which ensures high-quality GNSS data products (e.g., GNSS satellite ephemerides, Earth rotation parameters, velocities, etc.). For the present research, the GNSS data were processed using BERNESE software (developed by the Astronomical Institute, University of Berne, Switzerland) [18] through two different approaches (Table 1). The Network Approach adopts about nine reference stations of known coordinates, in addition to the estimating network, in order to set the terrestrial reference frame; the process uses double-differences of the GNSS observables (then a network of stations acquired simultaneously is required). The Precise Point Positioning Approach uses one-way GNSS observables; it does not require a network of stations acquiring simultaneously, but it strictly requires: GNSS data collected by the station to be estimated, precise GNSS orbits, related Earth rotation parameters, corrections for the clocks bias and drift, and specific mapping function for the atmospheric delay modeling (neutral part of the atmosphere).

**Table 1.** Main feature of the Global Navigation Satellite System (GNSS) data processing approaches. ERP: Earth rotation parameters. IGS: International GNSS Service. CLK: clock correction.


The acquired and interpreted ground monitoring data were used during the back analysis of the geomechanical modeling, which represents the final step of a complete geological, fluid-flow and geomechanical numerical modeling workflow (Figure 1). The 3D numerical modeling approach was based on a coherent and full dataset, including: 3D seismic survey, primary production and storage data, well-logs, petrophysical and fluid properties; in situ Image log and MDT stress tests, deformation and strength parameters from lab tests, logs and technical literature.

Starting from the detailed reservoir static model, describing all the main stratigraphic, structural, lithological and petrophysical characteristics of the reservoir, an extended (regional-scale) geological model was set up (around 450 10<sup>3</sup> cells). The latter effectively reproduces the main stratigraphic and

structural features of the investigation domain for the ground movement analysis, which includes the reservoir itself and its surrounding formations, up to the surface (Figure 2). A multiphase flow numerical model (FDM) was then set up based on the geological model and on all the available petrophysical and fluid properties together with well completion data. The reservoir consists in sand layers characterized by 18–23% porosity and 100–300 mD (milliDarcy) permeability ranging values; the permeability anisotropy at the reservoir scale reflects the occurrence of pelitic interlayers within the reservoir. The model was then initialized and calibrated to reproduce the pressure distribution in the reservoir during more than 60 years of production and storage history, considering data from 80+ wells. The pressure distribution in the reservoir is homogeneous as the reservoir has no internal compartmentalization; the cyclical pressure variations follow the cyclical gas storage volume variations that were compared with the ground movement measures during the InSAR and GNSS analyses. The dynamic model is periodically updated at the end of each injection/withdrawal phase as new pressure data are acquired from the monitoring wells.

**Figure 2.** 3D view of the regional geological model used for the land displacement analysis during geomechanical simulations, and of the reservoir model used for the flow simulations.

Finally, the ground displacement analysis was addressed via a stress–strain finite element (FEM) model, set up based on the extended geological study. The model was characterized by mechanical properties derived from lab and log data (Table 2). During the initialization phase, the initial pore pressure distribution and the original stress state were defined for all the investigated volumes according to the available data (i.e., initial static pressure, well logs and in situ tests). Dirichlet boundary conditions were imposed (i.e., zero displacement on the lateral and bottom boundaries, zero normal stress on the surface). Rock mechanics engineering analyses were developed according to the elastic-purely plastic constitutive law and adopting the Mohr–Coulomb failure criteria. As a matter of fact, different studies focused on ground subsidence induced by UGS for analogous case studies have shown the reliability of the elastic constitutive model in simulating the measured surface movement data [19–21]. It should be taken into consideration that the UGS-related pressure variations affect the formation both cyclically (following the withdrawal/injection phases) and in a very small time frame (5/6 moths): it clearly differs from a standard monotonic pressure decrease over some decades that typically characterizes the primary production, and the formations react accordingly. On the base of the preliminary sensitivity analysis results, the model of the case under study reacted within the elastic domain even at the end of primary production when the system experienced the maximum pressure variation. An adequate time stepping was set up accordingly for the UGS system analysis: the mechanical simulations were performed at the end of the primary production, and at the end of each production/injection period for all the InSAR monitoring time frames. Furthermore, a negligible variation of the petrophysical parameters due to pressure change was expected, in agreement with other UGS in analogous clastic formations studies [7,11,21]. Consequently, an uncoupled fluid-flow and stress–strain approach was adopted: the time and space pressure evolution represents the forcing function applied to the geomechanical model, whereas the petrophysical parameter values adopted in the fluid-flow model remain unchanged. The historical pressure evolution and the vertical ground displacement data were finally integrated in the numerical model via a back analysis approach: the adopted pseudo-elastic parameters of the reservoir and of the cap rock were fine-tuned to achieve a suitable match between simulated and measured movements. Previous experiences of the authors and the technical literature (please consider the abovementioned references) highlighted that the UGS-related ground movements for a clastic reservoir at about 1000-1500 m depth, placed in the specific geological contest under analysis (paragraphs 3.1 and 3.2), are mainly influenced by the stress–strain behavior of the gas bearing formations and, secondly, by the stress–strain behavior of the cap rock.


**Table 2.** Mechanical parameters used for characterizing the geomechanical model. E: Young's Modulus; ν: Poisson's coefficient; ρ: density; β: Biot's coefficient; c: cohesion; ϕ: friction.

#### **3. Results**

#### *3.1. Regional Geological Setting*

A deep knowledge of the geological characteristics of the examined formations is fundamental for the construction of accurate and representative 3D geological models and for the interpretation of ground monitoring data. Regarding the case study, a detailed geological analysis was developed at both reservoir and regional scales and it provided insight of the internal stratigraphic–structural architecture and of the heterogeneity of the investigated volume.

The studied UGS field is placed in the eastern Po Plain (North of Italy) (Figure 3). The Po Plain-Northern Adriatic region was extensively studied through exploration well data and 2D-3D seismic surveys acquired over the last 40 years (e.g., [22–25]). During Jurassic–Triassic time the region experienced east-west directed extensional phases which led to the formation of fault-bounded basins (e.g., [23,26]). In response to the subsequent Eurasia–Adria plate convergence, the rift-related framework was overprinted by Cenozoic compressional structures that are currently buried below the Po Plain sedimentary infill. The present-day architecture of the Po Plain is mainly the result of the Late Messinian–Pleistocene tectonic evolution; the latter led to the development of the Po Plain-Adriatic basin as a shared, complex foredeep separating the converging Southern Alpine (to the north) and Northern Apennine (to the south) chains, which locally collided in the Po Plain subsurface (e.g., [27,28]).

The buried Northern Apennines is an arc-shaped fold-and-thrust belt subdivided into three adjoining arcs named Monferrato, Emilia and Ferrara-Romagna arcs (from the west to the east, respectively) [29]. Its subsurface architecture consists of synclines and ramp anticlines associated with blind thrusts formed in response to a series of successive tectonic pulses that have led to the progressive North-Northeastward migration of the outermost thrust fronts (Figure 3). During the latest Miocene–Pliocene event, the eastern Po Plain has undergone a strong tectonic activity; the Late Pliocene–Pleistocene tectonic events led to the complete development of the outer fold and thrust system of Ferrara and to the uplift and tilting of the Bologna area [24]. The area of the UGS field is characterized by NE-verging, arcuate blind thrusts which involved the Mesozoic units, their pre-Mesozoic basement and the younger clastic infill of the Plio-Pleistocene foredeep (Figure 3). In the field area there is no evidence of the presence of thrusts dislocating the Pleistocene sediments or with reactivations during the Late Pliocene–Pleistocene times: none of the thrusts that created significant dislocation of the pre-Pliocene sequences propagated across the base of the Early Pleistocene deposits (e.g., [25,30,31]).

**Figure 3.** (**A**) Schematic structural map of the eastern Po Plain showing the geographic area of the underground storage of natural gas (UGS) under study (**B**). Interpreted seismic section showing the structural–stratigraphic architecture of the buried sequence across the area of interest (modified figure from [30]).

The Pliocene foredeep infill in the UGS field area includes the thick syntectonic, sand turbidites of the Porto Corsini Fm. (Early-Middle Pliocene) and Porto Garibaldi Fm. (Middle-Late Pliocene); meanwhile, the clay sequence of the Santerno Fm. (Pliocene–Early Pleistocene) was deposited above the

submarine paleo-highs and part of the foreland and foreland ramp. Finally, the late Pliocene–Pleistocene sedimentary infill recorded the gradual transition from deep-marine to continental environment (i.e., Asti Group), through the development of the basin-scale Po Plain prograding complex as a result of the east-ward progradation of the Po River delta. The overall Plio-Quaternary clastic sequence reaches a thickness of about 7-8 km in the deepest depocentral areas of the basin [22,24,32].

#### *3.2. UGS Gas Field Description*

The investigated gas field is located approximately 20 km NE of the Bologna city (Figure 3) at an average depth of 1300 m ssl. The gas primary production started in 1956 and, after two decades, the reservoir was converted into a gas storage system. Currently the system is managed at a maximum injection pressure close to the initial pressure of the reservoir with injection and production rates in the order of 3 to 25 MSm3/d and 5 to 40 MSm3/d, respectively; 50+ wells were drilled on the trap culmination for the storage operations and pressure monitoring.

The reservoir consists of a mixed (stratigraphic–structural) trap that combines an asymmetric anticline controlled by an NE-verging, blind thrust and the deposition of sandy levels onlapping on the structural high with local pinching-out features. The pinch-out closure over the structural high defines a non-deposition area in the SE part of the field. Minor, extensional faults, sub-parallel to the trap axis, affect the axial area. Neither the stratigraphic nor the structural features produced reservoir compartmentalization.

The reservoir sequence belongs to the Porto Garibaldi Fm. (Middle-Late Pliocene) and consists of an alternance of sand layers, meter-to-decameter thick, and clay layers, 1 to 10 m of average thicknesses. The caprock consists of clay and silty-clay of the upper member of the Santerno Fm.; the latter covers the entire reservoir area with an average thickness of 130 m ensuring the complete seal of the gas trap. Downward, the Porto Garibaldi Fm. is separated from the Porto Corsini Fm. (Early-Middle Pliocene) by the Santerno Fm. lower member.

In 2011, a 3D seismic survey was acquired and processed covering the entire UGS area (Figure 3). The approximate extension of the acquisition was 11.5 km × 7 km and the investigation depth almost 12 km. The results were integrated and interpreted together with data from the existing wells. The cross-section derived from the 3D seismic volume shown in Figure 4 highlights the geometry of the trap and of the main thrust bounding the reservoir to the North. The three main marker surfaces here represented correspond to the Top reservoir (Late Pliocene), Top Santerno (Calabrian, i.e., Early Pleistocene) and Top Prograding (Middle Pleistocene). The seismic data show that the thrust that bounds the trap ends within the Santerno Fm. without reaching the Top Santerno horizon and thus not affecting the Quaternary sequence as well as none of the faults inside the seismic volume. According to the literature, it confirms the absence of evidence of recent reactivation for the Neogenic faults present in the investigated area.

**Figure 4.** Cross section of the UGS area derived from the 3D seismic volume indicating the main stratigraphic horizons, the geometry of the trap and of the main thrust fault.

#### *3.3. Ground Surface Monitoring Results*

#### 3.3.1. SAR interferometry results

The huge archive of InSAR data (C-band) (Table 3) was analyzed via the PSP-IFSAR technique. Radarsat-1/2 data have provided a very long series of ground displacement measurements starting from 2003 with a temporal resolution of 24 days in both observation geometries. Starting from 2014, Sentinel-1 data have provided more detailed information on the last 5 years with a new acquisition every 6 days in both observation geometries. The satellite data cover an area larger than the field area, as shown in Figure 3.


**Table 3.** Main characteristics of the datasets. RSAT1: Radarsat-1; RSAT2: Radarsat-2; S1: Sentinel-1.

Based on ascending and descending results, for each sensor, the vertical and East–West planar movements were derived. The PSP measurements of the ground movement at each acquisition date allowed the evaluation of its mean velocity (as linear trend between the first and the last acquisition of each interferometric stack) and its mean amplitude (in case of sinusoidal behavior).

The mean velocity data show the overall stability of the field area with a slight subsidence trend; local zones with slight uplift or subsidence are evident in the vertical velocity map, with a notable subsidence phenomenon in the South–West area towards the city of Bologna (Figure 5). The mean vertical velocity in the area above the field ranges between −1.0 and −2.0 mm/year. The sub-millimetric

difference between the Radarsat and Sentinel results can be attributed to the different time spans they investigated (2003–2019 vs. 2015–2019, respectively).

**Figure 5.** Vertical (**a**) and horizontal (**b**) velocity maps calculated from Radarsat analysis and the field boundary projection on the surface (white line). The red triangle represents the reference point for Interferometric Synthetic Aperture Radar (InSAR) analysis. Positive values of the vertical and horizontal velocities mean upward and eastward movements, respectively.

Figure 6 shows the Radarsat maps of the mean vertical and horizontal amplitude of the ground displacement. The well-defined area of maximum vertical amplitude is located above the center of the UGS field (Figure 6a) and it coincides with the area of minimum horizontal amplitude (Figure 6b). Instead, the areas of the maximum horizontal amplitude identify two isolated peaks in the East–West direction around the center of the reservoir; zero amplitude values are detected along the North–South direction because of the limit of the InSAR technique.

**Figure 6.** Vertical (**a**) and horizontal (**b**) displacement amplitude maps calculated from Radarsat analysis and the field boundary projection on the surface (white line). (**a**) Location of the GNSS station and location of the points used for the displacement time series analysis.

Six representative points were selected in the InSAR-monitored area within and outside the field boundary (Figure 6a): their historical displacement along the vertical and the East–West directions were compared with the curve of the storage gas cumulative volumes of the whole field. The temporal evolution of the two ground displacement components for the P1, P2, P3 internal points shows a clear sinusoidal signal with a strong correlation with the trend of the gas cumulative volume curve in terms of amplitude and periodicity (Figures 7 and 8). The vertical displacement always shows the maximum and minimum peaks at the end of each injection and withdrawal period, respectively, with a delay of about 30 days, which could be attributed primarily to fluid-flow phenomena (i.e., time and space

propagation of the pressure sink in the reservoir and surrounding aquifer) (Section 2. Materials and Methods).

**Figure 7.** Comparison between the vertical displacement of the internal points from Radarsat (**a**) and Sentinel (**b**) analysis and the curve of the field gas cumulative volumes (orange line).

**Figure 8.** Comparison between the horizontal displacement of the internal points from Radarsat (**a**) and Sentinel (**b**) analysis and the curve of the field gas cumulative volumes (orange line).

The horizontal displacement shows alike clear sinusoidal signals: during the withdrawal periods P1-P2 points show positive trends (i.e., eastward direction) and P3 point shows negative trend (i.e., westward direction); vice versa, during the injection periods, P1-P2 points show negative trends and P3 point shows a positive one. The divergent responses can be attributed to the different locations of the points in respect to the area of maximum vertical displacement (i.e., reservoir center). The time series of the points outside the field area (P4, P5, P6) show displacement trends, oscillation characteristics and periodicity, both in vertical and horizontal directions, unrelated with the storage activity (Figures 9 and 10).

**Figure 9.** Time series of the vertical displacement of the external points from Radarsat analysis. The orange line represents the curve of the field gas cumulative volumes.

**Figure 10.** Time series of the horizontal displacement of the external points from Sentinel analysis. The orange line represents the curve of the field gas cumulative volumes.

Note that the Radarsat series are more sensitive to the long-term displacement trends whereas Sentinel series give a higher detail on the last monitoring years because of the different monitored time spans. The time series of the vertical displacement confirm a long-term trend of gentle surface downwarping in the field area.

#### 3.3.2. GNSS results

The permanent GNSS station was installed at the end of 2008 in the area of the storage plant located at the south-eastern edge of the gas field (Figures 3 and 6). The station was equipped with an ASHTECH UZ-12 receiver since the installation date in 2013, when it was replaced by a LEICA GR10 receiver. The GNSS antenna is a TOPCON, TPSCR4 model, with TPSH radome, dual frequency GPS (L1/L2).

The daily RINEX data of the GNSS permanent station were analyzed via the TEQC software (http://facility.unavco.org/software/teqc), an international standard for the pre-processing phase of GPS data and for the evaluation of their quality. The data quality was stated by using the TEQC parameters MP1 (multipath on L1) and MP2 (multipath on L2). Based on the values of the IGS reference stations, the values of MP1 (<0.5 m) and MP2 (<0.75 m) indicate a good quality station (Figure 11).

**Figure 11.** TEQC parameters: MP1 and MP2. The increase in the values of MP1 and MP2 for a short period in June 2016 was related to activities on the building where the GNSS is placed.

The results of the two processing approaches, the Network Approach and the Precise Point Positioning, are reported in Figure 12 in terms of coordinates time series: their consistency attests the robustness of the results. The areal and vertical trends of displacement indicate a marked NE-ward movement and a gentle subsidence of the monitored site, respectively.

**Figure 12.** Coordinates time series along the North, East and vertical component. Positive values mean northward, eastward and upward movements, respectively.

In order to discriminate the effects of the UGS activities along the three components of displacement, the signals ascribed to other known phenomena were removed. The typical seasonal signals of the GNSS time series are essentially due to load variations caused by the redistribution of fluid masses on the Earth's surface, (e.g., [33]). It includes changes in air mass, which results in changes in surface air pressure, in ocean level due to terrestrial and solar tides, wind and pressure atmospheric, and, above all, in soil moisture (surface hydrological load). Once removed the long-term velocity from the GNSS time series, the residual time series were obtained. The latter were compared with the displacements due to changes in the surface hydrological load (the so-called "surface hydrological mass loading") by using the data of the EOST Loading Service (http://loading.u-strasbg.fr/) and the data estimated by MERRA2 [34]. The seasonal vertical movement recorded by the GNSS station is largely due to the variation in the hydrological load (the hydrological model well followed the analytical seasonal model) and does not show an appreciable displacement-UGS correlation. On the contrary, the hydrological model does not explain the seasonal shifts observed in the planar components, attributed to the UGS activity (Figure 13). The two horizontal displacement components are seasonally dependent: during withdrawal periods the North–South component shows a positive trend (i.e., northward direction), and the East–West component shows a negative trend (i.e., westward direction), in agreement with the movement toward the peak of the maximum vertical displacement highlighted by InSAR data. On the contrary, during injection periods, the North–South component shows a negative trend, and the East–West component a positive trend. Furthermore, the main direction of planar movement is almost North–South (about N198◦).

**Figure 13.** Residuals time series of the GNSS station after removing the linear trend (from Network Approach): light blue line: residual displacement; yellow line: the curve of the field gas cumulative volumes.

#### *3.4. Geomechanical Simulation Results*

In the back-analysis phase of the mechanical model, the calibration of the pseudo-elastic parameters of the reservoir and the cap rock allowed a suitable match between the simulated vertical ground movements and the measured values discussed in the previous paragraph. Radarsat data were considered because of their longer period of acquisition. Figure 14 shows the comparison between simulated and measured values, in terms of relative variations within a single withdrawal/injection cycle, for the three points placed at surfaces within the boundary of the field (location in Figure 6) during a monitored time frame. For P2 point, a descending trend recognized by InSAR analysis, and not ascribable to UGS operations, was superimposed to the simulation results. The analysis of the InSAR data allowed the identification also of the areal extension of the ground surface cyclically affected by the UGS activities: Figure 15 shows an example of comparison between measured and simulated delta displacement maps at surface level, for a historical withdrawal period during UGS activity.

**Figure 14.** Comparison between measured and simulated vertical movements for the internal reference points.

**Figure 15.** Maps of measured (**a**) and simulated (**b**) land surface delta displacement within a historical withdrawal period (March-November). The geographic location of the points P1, P2, P3 is shown in Figure 6.

#### **4. Discussion**

#### *4.1. Ground Monitoring Data and Comparison with UGS*

Since 1950s, the eastern Po Plain area, where the studied field is located, has been experiencing a widespread and strong subsidence mainly due to anthropic causes, with a minor contribution (of a few mm/year) from ongoing natural processes (e.g., tectonics, sediment consolidation) [8,35,36]. In particular, the subsidence evolution in this region mostly depended on the groundwater pumping [8]. This phenomenon was particularly relevant in the Bologna area, where a subsidence rate in the order of several cm/year was historically measured; subsequently, the impressive subsidence reduction monitored during recent decades was a direct result of the reduction in groundwater pumping activities [36].

The huge spatial and temporal extension of the InSAR monitoring allowed the investigation, not only of the UGS field area but also of the surroundings; it ensured the consistency of the measured displacements with respect to regional trends over a large time span. The analysis and interpretation of the InSAR data have showed that:


The analysis and interpretation of the GNSS data have showed that:


In conclusion, the integration between GNSS and InSAR data provides a full and coherent panorama of the ground movement in the monitored area.

Note that the present analysis investigated the effects of the storage activities and did not deeply analyze other variables such as seasonal temperature variations, groundwater withdrawal or large-scale natural subsidence/uplifting.

#### *4.2. Discussion about the Geomechanical Simulation Results*

The pseudo-elastic parameters adopted for the mechanical characterization of the reservoir and the cap rock were derived from laboratory tests (triaxial compression tests under isotropically consolidated undrained conditions), from the interpretation of well logs (sonic and density logs) and corroborated by data from the technical literature. The obtained Young's modulus values differ by almost one order of magnitude: from few GPa (values from the lab data interpretation) to 8–10 GPa (values from the log data interpretation). As documented by analogous case studies reported in the technical literature [11,19], the mechanical behavior of clastic formations at a depth of around 1000–1500 m ssl could exhibit an important (non-linear) influence of the strain on the formation stiffness, especially in case of soft rocks [37,38]. Experimental values of the deformation parameters (i.e., Young's modulus) increase as the imposed strain levels decrease from small strain (i.e., the values obtained from the stress–strain measurements in a rock mechanical test) to very small strain (i.e., values obtained from the interpretation of acoustic acquisition) [39,40]. Furthermore, if acoustic acquisitions are available at wellbore level, the scale effect can also be taken into consideration.

The discrepancy in the experimental values could generate uncertainty in defining the real deformation behavior of the formations; it is trivial to say that high variability of input parameters leads to high uncertainty in model prediction.

For the investigated case study, the real deformation behavior of the formation was identified via the back-analysis of the mechanical model. The formation stiffness that resulted from the calibration process is around three times higher than the values from lab tests, in agreement with the lithologies under analysis and with the in-situ stress condition during UGS operations [11,19]. As a theoretical explanation, during primary production, the deformation behavior of the normally consolidated formations is controlled by the loading static elastic modulus EI (in the case presented in this paper, EI are the values obtained from the lab tests). In this phase, the stress–strain path follows the critical state line (CSL). Due to the formation compaction caused by storage activities, the formations could become slightly over-consolidated and their elastic deformation response during the storage phase is usually ruled by the unloading/reloading static elastic modulus EII, around 2-4 times higher than EI.

The good match obtained between measured and simulated ground displacements assuming a constant unloading/reloading elastic modulus over almost two decades of UGS operation denotes a sound and stable system response in time.

Once calibrated, the geomechanical model represents an effective tool for forecasting future scenarios.

#### **5. Conclusions**

The paper proposes a fully-integrated and multi-physics approach for the analysis of the ground movement related to a gas storage system (eastern Po Plain area). The gas field is hosted in a clastic Pliocene sequence involved in a mixed trap. The trap combines the onlap of sand levels and a structural high associated with an NE-verging blind thrust that ends within the Pliocene sequence without propagating up to the Quaternary clastic sequence. The InSAR data excludes surface evidence of an active growth of the buried anticline during the entire monitored period.

The monitoring of the ground movements was achieved via both the SAR interferometry (PSP-IFSAR) and the GNSS continuous technologies. The two complementary technologies provided quantitative information about vertical and planar (North–South and East–West) components of the surface movements (with millimetric accuracy) and the subsidence/uplift areal extension of the UGS influence. The InSAR data measured the East–West planar and vertical components; the robustness of the adopted methodology is attested by the use of two different sensors (Radarsat-1/2 and Sentinel-1) to get independent and independently processed measures, and by the consistency of the results. On the other hand, the GNSS station measured the punctual displacement component associated with the UGS activity along the North–South direction, missed by the InSAR acquisition. The information from GNSS and InSAR data provided a full and coherent panorama of the ground movement in the monitored area: in particular, the identification of a regional trend was used for interpreting the effects of UGS operation.

The reliable ground movement monitoring plan and the detailed geological analysis, both at reservoir and regional scales, are key information for the geological, fluid-flow and geomechanical numerical simulation process. The multidisciplinary simulation approach deeply investigated and reproduced the main geological and structural features of the system together with the fluid-flow and poro-mechanical processes induced by UGS activities. In particular, almost two decades of displacement data effectively supported the calibration process of the geomechanical model, which showed a sound and stable mechanical response of the system in terms of induced subsidence/uplift for a given variation of pressure.

Furthermore, the two calibrated dynamic and geomechanical models represent the most reliable tool in forecasting the performance and in verifying the safety of the UGS system under different management conditions.

**Author Contributions:** All authors collaborated for the preparation of the original draft and on revising the manuscript; A.M., C.B., G.C., G.T., V.R. proposed the main structure of this study; F.T., A.M., G.C. dealt with InSAR data analysis, C.F., A.M. dealt with GNSS data analysis, C.B., G.T. dealt with seismic data analysis, G.C., V.R. dealt with geomechanical analysis. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to express their gratitude to two anonymous reviewers for their comments and suggestions that considerably improved the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Technical Note* **Estimating Vertical Land Motion from Remote Sensing and In-Situ Observations in the Dubrovnik Area (Croatia): A Multi-Method Case Study**

#### **Marijan Grgi´c 1,\*, Josip Bender <sup>2</sup> and Tomislav Baši´c <sup>1</sup>**


Received: 29 September 2020; Accepted: 26 October 2020; Published: 29 October 2020

**Abstract:** Different space-borne geodetic observation methods combined with in-situ measurements enable resolving the single-point vertical land motion (VLM) and/or the VLM of an area. Continuous Global Navigation Satellite System (GNSS) measurements can solely provide very precise VLM trends at specific sites. VLM area monitoring can be performed by Interferometric Synthetic Aperture Radar (InSAR) technology in combination with the GNSS in-situ data. In coastal zones, an effective VLM estimation at tide gauge sites can additionally be derived by comparing the relative sea-level trends computed from tide gauge measurements that are related to the land to which the tide gauges are attached, and absolute trends derived from the radar satellite altimeter data that are independent of the VLM. This study presents the conjoint analysis of VLM of the Dubrovnik area (Croatia) derived from the European Space Agency's Sentinel-1 InSAR data available from 2014 onwards, continuous GNSS observations at Dubrovnik site obtained from 2000, and differences of the sea-level change obtained from all available satellite altimeter missions for the Dubrovnik area and tide gauge measurements in Dubrovnik from 1992 onwards. The computed VLM estimates for the overlapping period of three observation methods, i.e., from GNSS observations, sea-level differences, and Sentinel-1 InSAR data, are −1.93 ± 0.38 mm/yr, −2.04 ± 0.22 mm/yr, and −2.24 ± 0.46 mm/yr, respectively.

**Keywords:** Dubrovnik; GNSS; InSAR; satellite altimetry; Sentinel-1; tide gauges; vertical land motion

#### **1. Introduction**

Displacements of the Earth's surface are generally caused by spatio-temporal varying geophysical processes. These processes include slow and steady Earth activities such as Earth's core motion, glacial rebound, and plate tectonics, which generate secular movements of the land, periodic activities such as Earth rotation and tides, which lead to periodic responses of the Earth surface, and unpredictable occurrences such as earthquakes often induce major land motion [1]. In addition, the Earth's surface is affected by human activities directly or indirectly e.g., through natural resource exploitation, urbanization and constructional works, and environmental degradation along with global warming [2].

Monitoring of the vertical land motion (VLM) has been a well-addressed geodetic task performed from the 1970s using different geodetic methods. The first VLM estimations were computed from the leveling and/or tide gauge data, which commonly resulted in VLM trends at a limited number of sites (e.g., [3–5]). In the late 1980s and during the 1990s, a point VLM estimating was extended for the computations based on the Global navigation satellite system (GNSS) data that were obtained continuously or repeatedly, usually within the geodynamic campaigns at the sites of the interest (e.g., [6–9]). The GNSS method was often supplemented with the other technologies such as Satellite Laser Ranging (SLR) or Very Long Baseline Interferometry (VLBI) (e.g., [10,11]) as well as DORIS observations (e.g., [12]). Additionally, the 1990s have brought the first studies on the VLM calculations based on the comparison of the relative sea-level change rates derived from tide gauges, which are related to the adjacent land so they integrate both sea-level and VLM change, and absolute sea-level rates computed from satellite altimetry, which are referred to in the global geodetic reference frame and are not influenced by the VLM (e.g., [13,14]). At the same time, the first studies based on the Interferometric synthetic aperture radar (InSAR) technology have shown the potential of the technology for area land motion monitoring (e.g., [15,16]). With further advances of all the satellite-borne and in-situ technologies in the 2000s, the Earth surface motion monitoring became an easier and more common geodetic task.

Today the VLM monitoring mostly relies on the integration of globally available and reliable InSAR and sparse, but very accurate, GNSS observations, which are often compared against the other geodetic observation methods (e.g., [17–23]). InSAR technology gives high spatial resolution information on land motion over the large area with relatively high temporal frequency and up to subcentimeter-accuracy [19]. The surface displacements derived from InSAR are given only as one-dimensional variables along the radar line-of-sight (LOS), which can be converted into vertical displacements. The VLMs derived from InSAR in this study were evaluated and compared against VLMs derived from (1) GNSS, and (2) differences of the sea-level change rates derived from tide gauge measurements and satellite altimetry. The VLMs were computed for the periods for which the data existed, i.e., from 1992 for computations from sea-level differences, from 2000 for the GNSS estimates and from 2014 for computations using InSAR technologies. Different trends were afterwards compared and analyzed for the observation overlapping period from 2014 to 2020.

#### **2. Study Area and Previous Research**

A multi-method VLM analysis was performed on the area of the Dubrovnik seismic zone, which was identified by [24] using the combination of the deterministic and probabilistic methods as the seismic zone of the highest expected seismic peak in Croatia. The selected geographical area encompasses a coastal area that serves as the transition zone between the Adriatic Platform at the southeast and the Dinarides at the northwest (Figure 1). As shown in Figure 1 the transition zone along the Dubrovnik and Adriatic faults is very seismically active with more than 50 earthquakes of magnitude from 4.0 to 6.0 recorded over the last hundred years [25,26]. Furthermore, the historical data show that eight earthquakes of intensity IX or X ◦ MCS were registered over the 15–17th centuries [26]. The one from 6 April 1667, was one of the most significant natural hazards in the broader area [26].

Dubrovnik area is known for the natural and cultural heritage, which attracts tourists, intensive construction projects, and other economic activities (see [27]). Being of great value, the area was studied by many scientists. Ref. [28] derived VLM for Dubrovnik from the comparison of sea-level trends computed from tide gauge observations at Dubrovnik site and from Topex/Poseidon satellite altimeter mission. The study has shown rates of land subsidence at −0.7 ± 0.8 mm/yr defined for the period 1993–2001. In 2007, ref. [29] computed VLM rate of −5.8 ± 2.0 mm/yr for the period 1993–2001. The large differences between the trends computed by two studies could be due to the additional satellite altimeter missions ERS-1/2 used by [29]. The latter study has also reported on the VLM rate of −5.31 ± 0.52 mm/yr computed from Global Positioning System (GPS) observations.

In 2010, ref. [30] reported on the VLM trend of −0.7 ± 0.5 mm/yr using the GNSS data. Two years later, ref. [31] reported on the extended study on the VLMs in the Mediterranean, which has shown VLM rates in Dubrovnik of −0.50 ± 0.38 mm/yr, −0.32 ± 0.18 mm/yr, −0.92 ± 0.31 mm/yr, and −0.60 ± 0.17 mm/yr from classical comparison of tide gauge trends and altimetry, advanced approach of the latter, GNSS observations [32], and glacial isostatic adjustment (GIA) model SELEN [33], respectively. The same procedures and VLM rates were later used for research studies of the sea-level studies of the same area (see e.g., [34,35]).

**Figure 1.** The study area shown in a topographic map (ETOPO1 topography/bathymetry model [36]) related to a broader geographical area showing seismic activities (epicenters of all registered earthquakes of magnitude >4.0 from 1990, shown by '+'; data downloaded from [25]) and known plate tectonics (border of the Adriatic tectonic plate shown in red; according to [37]); pink circle represents the tide gauge station whereas yellow square shows collocated GNSS station.

#### **3. Data**

Three different data sets available globally were downloaded and used for the independent and conjoint analysis of VLM.

#### *3.1. InSAR and GNSS*

GNSS and SAR techniques present complementary characteristics regarding their resolution, both spatial and temporal, and their sensitivity to components of land motion. GNSS data recorded by continuously operating stations have high temporal and low spatial resolutions, while InSAR has a lower temporal resolution and a large number of measurements [38]. As radar is only capable of measuring path length differences in LOS direction, InSAR results need to be interpreted with care. A three-dimensional displacement vector is projected to one slant-range component [39]. The LOS is most sensitive to VLM, due to the usually steep incidence angle, but the near-polar orbital plane of SAR satellites means that it also contains any motion in the east-west and, to a much lesser extent, north-south directions [40]. While GNSS measures absolute movements in the geocentric frame, the Persistent Scatterer Interferometry (PSI) is a differential method, which means that the estimated velocities are relative to a reference point. For that reason, PSI is applied to study local movements in limited areas [38]. To provide absolute information, the motion of the reference point must be known with sufficiently high accuracy [41].

The SAR dataset used in this study consists of 75 Sentinel-1 images (VV polarization) in an ascending orbit (track 73) acquired over a period of six years (2014–2020). Data have been processed using the open-source European Space Agency's (ESA) Sentinel Application Platform (SNAP) and Stanford Method for Persistent Scattererrs (StaMPS) software packages, following well-established and described processing scheme (e.g., [42]). Data preparation for PSI analysis, coregistration and interferogram formation was carried out in SNAP. The topographic phase was removed using SRTM (Shuttle Radar Topography Mission) 1 arc-second digital elevation model (DEM). The images were cropped on the area of the interest to limit the processing time and computing resources. The single master PSI approach (see Figure 2) was utilized because area of the interest consists mostly of an urban area, while the non-urban area is covered in low vegetation and rocks, thus a sufficient number of PS pixels is expected [43]. This approach is based on the identification of a set of targets in the area of interest exhibiting stable radar response over time, named permanent scatterers (PS), which can be used to estimate and remove the atmospheric disturbances affecting the radar images [44]. Although phase stability is what constitutes the PS pixel, initial PS selection was done using the amplitude dispersion threshold, because of the statistical relation between amplitude and phase stability [45]. The deformation phase was separated from the atmospheric phase and noise by filtering in time and space, following the assumption that deformation is correlated in time, the atmosphere is correlated in space but not in time, and noise is uncorrelated in space and time [46]. StaMPS analysis resulted in an average PS density of 550 PS/km2.

**Figure 2.** Spatio-temporal relation-baseline plot; master image (red dot) and slave images (black dots).

The Dubrovnik-2 GNSS station (see subsection below), which is part of EUREF Permanent GNSS Network, is used to correct the precise high-spatial-resolution PS velocities and reference the motions with respect to an absolute reference frame. With only one GNSS station available, a calibration offset (*TREF*) was determined by computing the difference between the weighted mean displacement rate of the PS grid-cell where the GNSS antenna is grounded (*VLMPS*\_*GNSS*) and the GNSS velocity (*VLMGNSS*) projected to LOS direction [47]:

$$T\_{REF} = V L M\_{PS\\_GNSS} - V L M\_{GNSS}.\tag{1}$$

The computed InSAR velocities that were related to LOS direction, *LOSdisplacements*, were ultimately recomputed to the vertical displacements, *VLMInSAR*, using the following formula (e.g., [48–51]):

$$VLM\_{InSAR} = \frac{LOS\_{displacements}}{cosθ} \tag{2}$$

where *θ* is the InSAR incidence angle expressed relative to a flat horizontal terrain. For the InSAR data obtained in this study the angle was 36◦.

#### GNSS Data

GNSS observations give direct insight on three-dimensional land movements. GNSS solutions used within this study were retrieved from SONEL [52], which provides daily data computed by Nevada Geodetic Laboratory 2019 (NGL19) [53]. Two stations in Dubrovnik (see Figure 1), Dubrovnik-1 and Dubrovnik-2, obtained measurements for the period 2000–2012 and 2012–2019, respectively. Due to the different locations of the sites, the GNSS measurements were analyzed independently by estimating the linear trend for each site/period. Data obtained at the Dubrovnik-1 station were corrected for the offset of 15.2 mm for the period from 2007 caused by the GNSS instrument replacing at the site. Data gaps shorter than three months were filled using the average interpolator.

#### *3.2. Sea-Level Data*

Monitoring the sea-level usually relies on two technologies: in-situ tide gauge measurements and radar altimeter observations from satellites. Tide gauges capture relative sea-level change, i.e., sea-level change influenced by geophysical and the other characteristics of the site, including its geodynamics and geotectonics that produce vertical motion. On the contrary, satellite altimetry provides the absolute sea-level, which is referred to as the geocentric Earth model or the surface (ellipsoid). Therefore, the differences in the sea-level change trends computed from tide gauges and satellite altimetry can be attributed to VLM. More on this topic is discussed in [35,47,54], etc.

For this study, monthly tide gauge measurements at Dubrovnik station (see Figure 1) were downloaded for the altimeter period (1992-onwards) from Permanent Service for Mean Sea-Level (PSMSL) [55] and obtained from the Hydrographic Institute of the Republic of Croatia for additional two years of the observations. The tide gauge site in Dubrovnik is well collocated with the GNSS site by leveling methods, showing no variation over time regarding the vertical difference. Hence, the VLM difference at the tide gauge and at the GNSS site was considered not to have existed in this study. The tide gauge measurements were corrected only by the Dynamic Atmospheric Correction (DAC) from MOG2D model that accounts for high frequencies of the barotropic forcing and low frequencies of the inverse barometer forcing caused by wind and pressure [56].

All available satellite altimeter measurements were retrieved for the study area for the same observation period from Radar Altimeter Database System (RADS) [57]. These include data captured by ERS-1, ERS-2, Envisat, TOPEX/Poseidon, Jason-1, Jason-2, Jason-3, Cryosat-2, Saral, and Sentinel-3. The reference surface of the satellite altimeter data was defined as the mean sea-level DTU15 [58], which enables the data to be easily transformed to the other reference surface such as the geocentric ellipsoid or the global geoid model. Table 1 lists the main conditions the data have had to satisfy along with the main corrections applied to the data before processing. The data processing was conducted based on the methods explained in detail in [35]. That included computing the monthly solutions in resolution 0.01◦ × 0.01◦ for the altimeter period in the Dubrovnik area by using Inverse Distance to a Power interpolator. The sea-level change from altimetry at the tide gauge site in Dubrovnik was thereafter computed by bicubic interpolation from monthly grids.

Finally, a linear trend of one-dimensional crustal displacements (or VLM), i.e., vertical land velocity, was computed from the differences of the monthly sea-level estimates captured by tide gauges and satellite altimetry. A simple equation was used:

$$VLM = SL\_{TG} - SL\_{SA\prime} \tag{3}$$

where *SLTG* is a linear trend defined from tide gauge measurements, whereas *SLSA* is a linear trend computed from satellite altimetry. Such computations were performed for three different periods following the GNSS and InSAR data availability.

**Table 1.** Data quality constraints and corrections/models applied onto the radar satellite altimeter data in pre-processing.


#### **4. Integrated VLM Computation Procedures**

Figure 3 presents an integrated approach used in this study to determine the VLM in the Dubrovnik area. Three independent assessments were performed from sea-level measurements, GNSS data, and InSAR observations. A simple linear estimator was used to determine VLM trends. All the processing of the data was conducted based on the procedures described above.

Thorough analyzes of all the input data were conducted to meet the standards of the procedures and to enable further comparisons. Although defined in different reference frames and from different sources, because the comparison is based on trend analysis, the results are comparable.

**Figure 3.** An integrated approach on accessing the vertical land motion (VLM) in the Dubrovnik area based on different input measurements.

#### **5. Results and Discussion**

The processed and filtered data used in this study along with the trends derived from that data by linear interpolator are shown in Figure 4. Subplots A–D of the figure show the data availability and their variation over time. Subplots A and B refer to the computations of vertical displacements at the tide gauge site based on the sea-level data, i.e., comparison of the observations captured at tide gauge and the interpolated sea-level heights for the same point from the monthly grid solutions that were computed based on the available satellite altimeter data using Inverse Distance to a Power interpolation method (see e.g., [35]). Good tide gauge data availability (with a few data gaps fulfilled using average interpolator) and full data availability of satellite altimeter data enabled the sea-level comparison over the whole altimeter period, starting from 1992. The two data sets, as expected, generally show the same variation pattern with relatively good agreement. The continuous disagreement of two data sets defined as the trend difference hence reveals the vertical land displacement velocity of the ground the tide gauge is attached to. The difference is well-pronounced, giving the trend of −0.61 ± 0.45 mm/yr for a whole observed period. The sea-level measurement differences further point to land subsidence at the rates of −1.21 ± 0.41 mm/yr for the period starting from 2000 when GNSS became available, and, finally, to the rate of −1.93 ± 0.38 mm/yr for the period from 2014 when Sentinel-1 started capturing data.

Following the sea-level trend differences, the daily GNSS solutions at Dubrovnik stations resulted in confirmation of the land subsidence at the trend of −1.46 ± 0.24 mm/yr and −2.04 ± 0.22 mm/yr for the period 2000–2012 and 2011–2020, respectively (Subplot C of the Figure 4). The data obtained for the first time period of the observation have a significant data outage starting from December 2005 until July 2007 due to the instrument malfunction. The new instrument was installed in 2007 with an offset of 15.2 ± 0.2 mm from the previous one at the same station, so the later-obtained data were accounted for the offset. As the data outage at the GNSS station lasted for almost two years, the measurements were not fulfilled using any interpolation method. However, a newely-installed instrument has shown a similar variation pattern so the trend computations for the whole observation period of the instrument, 2000–2012, are firm. In 2012, the measurements at the first GNSS site were discontinued with the other instrument being installed in 2011 at the same site, a few meters away. The observation overlap from the two instruments ensured flawless transition and continuity of the measurements. To enable the direct comparison of the trends computed from the three methods, the VLM was also estimated for the InSAR data period, 2014–2020, solely. Those computations have shown a similar trend to the one for the period from 2011, i.e., −2.04 ± 0.22 mm/yr.

Lastly, the land subsidence was detected from the InSAR measurements shown in the Subplot D of the Figure 4 with the rates of VLM at −2.24 ± 0.46 mm/yr. The presented trend was averaged over the available InSAR estimates in the tide gauge sourroundings of 100 m. In addition to the point estimates of the VLM, the InSAR observations provided the VLM surface for the whole study area for the period from 2014 to 2020. The estimated land motion trends are shown in Figure 5. The trends show the land subsidence of the whole area up to −7.12 mm/yr with the average of −1.56 mm/yr for the whole area and no uplift noticed in the whole study area. The standard error of the measurements is centered at the 0.26 mm/yr for the whole area and ranges up to 0.50 mm/yr. It can be stressed that the computed subsidence of the coastal area shown in the figure directly amplifies the general effect and impacts of the sea-level rise in the study area.

**Figure 4.** Processed observations and computed VLM trends; subplot **A** presents original filtered and processed altimeter and tide gauge data; subplot **B** shows the differences between the sea-level data presented in subplot **A** and the trends for three defined periods-from 1992 (altimeter era), from 2000 (Global Navigation Satellite System (GNSS) observation period in the study area), and from 2014 (Interferometric Synthetic Aperture Radar (InSAR) observation period); subplot **C** presents GNSS data obtained in Dubrovnik at two GNSS stations for the periods 2000–2012 and 2011–2020 along with the computed VLM trends for two observation periods (the same trend is computed for InSAR period, from 2014); subplot **D** presents an averaged Sentinel-1 InSAR measurements over the period 2014–2020 at tide gauge location with the VLM trend line.

Overall, the results of this study confirm the subsiding of the study area, which is also addressed by the previous studies listed above. However, compared to those studies, this study has significantly longer observation periods for each of the observation methods, so the differences were expected. Large differences of the computed VLM estimates derived from sea-level data comparing to e.g., [29] that claimed VLM rate of −5.8 ± 2.0 mm/yr for the period 1993–2001 can be due to the utilization of the re-tracked data in this study with improved orbital parameters and better data calibration as well as better availability in the coastal areas, and using the improved methodology that has been evolved in the last decade (e.g., [59]). The differences of the GNSS vertical displacement estimates given in this

study and in the previous studies, e.g., in [32], which were at the rate of −0.92 ± 0.31 mm/yr, are due to the longer observation period in this study, better data availability (original data available) and the improved data processing methods computed by [52].

**Figure 5.** VLM trends computed from Sentinel-1 InSAR data for the period 2014–2020; each cell of approximately 100 × 100 m is represented by the weighted average of the VLM within the cell (3–15 estimates); the positions of the GNSS and tide gauge station are shown by the cross with yellow and cyan circle, respectively; the black box marks the old town of Dubrovnik; the orthophoto background map was downloaded from [60].

It should be emphasized that an area of the special architectural and historic interest within the study area, the old town of Dubrovnik on peninsula shown in Figure 5 (marked by the black box), shows firm trend of land subsidence at the average rate of −1.96 ± 0.44 mm/yr with no significant anomalies. However, some larger land subsidence at rates higher than 3 mm/yr are evidence for some smaller parts of the broader study area, possibly due to the more intensive underground water change, erosions, or slow landslides. Such areas could be marked as the potential focus areas in future studies.

#### **6. Conclusions**

Different technologies provide insight into vertical ground motion with limitations that are most often related to the time variability or to the spatial distribution of the data. Comparing at the single points in Dubrovnik, i.e., at the tide gauge site and the collocated GNSS site with the relation well established and accounted for using the repeated leveling measurements over time, for the period from 2014 to 2020, the computed VLM velocities are −2.04 ± 0.22 mm/yr, −1.93 ± 0.38 mm/yr, and −2.24 ± 0.46 mm/yr, from GNSS observations, sea-level differences, and Sentinel-1 InSAR data, respectively. A greater insight onto the surface VLM trends is given by the latter technology, showing the average trend of −1.56 mm/yr for the whole study area, and pointing toward some specific areas with the VLM motion trend exceeding −5 mm/yr.

GNSS data used in this study gave continuous high-precision and high-frequency insight into VLM outperforming the rest of the measuring techniques regarding the precision and measurement frequency. However, the VLM estimated at the GNSS station is point-limited. Further, the VLM estimates of the lower frequency and lower precision than of the GNSS measurements were computed from the trend difference of two sea-level measuring techniques. Such measurements could be employed for the areas without GNSS data available (and with tide gauge data available) or where longer time periods need to be considered. Finally, the InSAR technology provided the surface vertical land velocities over the study area with good reliability and coverage making it a great tool for land monitoring. Additionaly, as the VLM contributes to the impact of the sea-level rise, which is evident for the study area, the InSAR technology could ensure reliable sea-level rise risk assessment at the coast.

**Author Contributions:** conceptualization, M.G. and T.B.; methodology, M.G. and J.B.; software, J.B.; validation, M.G.; formal analysis, M.G. and T.B.; investigation, M.G. and J.B.; resources, J.B.; data curation, J.B.; writing—original draft preparation, M.G.; visualization, M.G. and J.B.; supervision, T.B.; project administration, M.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Technical Note* **Land Surface Subsidence Due to Mining-Induced Tremors in the Upper Silesian Coal Basin (Poland)—Case Study**

#### **Paweł Sopata, Tomasz Stoch, Artur Wójcik and Dawid Mroche ´n**

Faculty of Mining Surveying and Environmental Engineering, AGH University of Science and Technology, 30-059 Cracow, Poland; Pawel.Sopata@agh.edu.pl (P.S.); awoj@agh.edu.pl (A.W.); dawmro@agh.edu.pl (D.M.) **\*** Correspondence: tomst@agh.edu.pl

Received: 30 October 2020; Accepted: 27 November 2020; Published: 30 November 2020

**Abstract:** Seismic phenomena threaten land-based buildings, structures, and infrastructure and can transform land topography. There are two basic types of seismic phenomena, namely, tectonic and anthropogenic, which differ mainly in epicenter depth, surface impact range, and magnitude (energy). This article shows how a land surface was changed by a series of seven rock mass tremors of magnitude ML = 2.3–2.6 in March–May 2017. Their immediate cause was the "momentary" acceleration of void clamping, which was activated by local and short-term seismic phenomena caused by human activity. The induced seismic events resulted from the geological structure of the rock mass, which in the specific region of examination was classified as being highly prone to mining tremors. The authors focused on describing vertical surface displacements in the Upper Silesian Coal Basin in the south of Poland. The surface deformations were identified using DInSAR technology, which allows quasi-continuous monitoring of large areas of land surface. The present research used freely available data from the Copernicus Program and seismic data from the European Plate Observing System.

**Keywords:** DInSAR; mining-induced tremors; land surface deformation; Upper Silesian Coal Basin

#### **1. Introduction**

The seismic phenomena that occur in nature are associated with geological structure, tectonic-plate movement, volcanism, and human activities. The main cause of earthquakes is sudden release of stress in the Earth's crust with displacements of rock mass layers. Seismic phenomena of anthropogenic origin are a decline of equilibrium between the external and internal forces that affect rock mass elements and are caused by human activities [1,2]. Such seismic phenomena are caused directly by voids in the rock mass due to the extraction of mineral resources, and they are common in the deep mining of gold [3], copper-ore mining [4], hard coal mining [5], natural gas extraction [6], and open-pit lignite mining [7], for example. With favorable geological structure of the overburden, disturbing the force equilibrium in a three-axis stress configuration may cause abrupt cracking and movements of the rock mass toward the post-excavation caverns. High-energy mining tremors with energies of at least 10<sup>5</sup> J [8] pose a real threat and can be felt as far as several kilometers from the epicenter [9]. Numerous scientific studies indicate a significant relationship between strong mining tremors and damage to residential buildings and building facilities [3,10–14].

The development of interferometric synthetic aperture radar (InSAR) technology has led to new possibilities for evaluating the surface effects of rapid seismic phenomena [15,16]. InSAR has also been used widely to monitor land displacement processes related to landslides [17,18], mining-induced subsidence [19–21], subsidence induced by aquifer system drainage [22], upward movement of post-mining areas [23–26], and land subsidence due to groundwater extraction [27,28]. Satellite radar images allow assessment of the extent of the surface deformation area, as well as the vertical

and horizontal displacements. Previous research conducted within the Upper Silesian Coal Basin (USCB) [29–31] and the Legnica-Głogów Copper District (LGCD) [32–34] in the south of Poland indicates a risk of tremor-induced troughs from strong mining tremors (ML ≥ 3). In such cases, the surface deformations have relatively high dynamics (3–4 days in the LGCD), the maximum subsidence can amount to several centimeters, and the troughs can extend as far as 2 km [32,33].

A difficulty to date in observing the surface effects of mining-induced tremors has been due to having to make point-based observations of the surface on either profiles (observation lines) or scattered points using GPS technology [35]. Another difficulty is the inability to predict the time and (in most cases) the exact location of the epicenter to facilitate detailed measurements before and after a seismic event. As such, traditional land surveying measurement methods have led to only isolated cases of identified land displacements due to mining tremors [9,36].

The present paper presents the effects of a series of seven rockmass tremors ofmagnitude ML = 2.3–2.6 on a land surface. These changes occurred between March and May of 2017.

#### **2. Materials and Methods**

#### *2.1. Area of Interest*

The area of interest is the mining area of a hard coal mine operating in the northeastern part of the USCB in the south of Poland (Figure 1). In the USCB (approximately 5800 km2 within the boundaries of Poland), hard coal has been extracted since the second half of the 17th century [37]. Currently, the hard coal output is approximately 50 million Mg per year (in 2018). In the study area, coal was extracted from two seams (Figure 1) in the years 2001–2017. The first one was extracted along 6 longwalls at the depth of about 560 m in 2001–2007. The second seam was extracted at the depth of about 640 m along 3 longwalls in 2008–2014. After two years, the last longwall of seam no. 2 was extracted in 2016–2017. The phenomenon of land surface deformation discussed in this article is related to the extracted of this longwall. The USCB is characterized by high seismic activity, with more than 1000 tremors per year of magnitude ML of at least 1.5 [38] and several thousand minor tremors reported. The USCB is a highly urbanized region of Poland, where several cities (including Katowice, Sosnowiec, Gliwice, and Zabrze) form the Upper Silesian Conurbation inhabited by more than 2 million people.

**Figure 1.** Area of Interest: (**A**) location of the research area (**B**) mining area (**C**) mining operations in seam 1 and seam 2 in 2001–2017 and localization of geological cross section.

The changes in land surface morphology discussed herein occurred in a region that is prone to unfavorable dynamic phenomena covering high-energy tremors. They pose a real threat and can be felt by residents in the area, who report consequential damage to buildings and infrastructure [39]. This situation is impacted by both the rock mass geological structure (including fault zones; Figure 2) and post-excavation effects (coal deposits extracted from above seams cumulating the stresses from previous hard coal exploitations; Figure 1) and exploitations depth [39]. The main lithological reason is the presence of thick and stiff sandstone beds in deposit roofs [40,41] presented in Figure 2. In the carbon formations (between the hard coal seams), the thickness of the rigid sandstone layers reaches approximately 35 m [39]. In the immediate roof (above the hard coal deposits), the height of these layers (i.e., fine- and medium-grain sandstone) reaches more than 85 m locally. Tremor-generating qualities are also shown by high-strength sand slates in this area [40]. The simplified geological cross-section is shown in Figure 2, and its location in Figure 1.

**Figure 2.** Simplified geological cross section.

#### *2.2. SAR Data Acquisition and Processing*

The present land surface deformations were described and analyzed based on freely available data that came from the Sentinel-1 satellite mission of the Copernicus Program financed by the European Space Agency and the European Commission. An exemplary phase interferometric band showing land surface deformations related to mining operations in the USCB region over a 6-day period is shown in Figure A1. The seismic data (i.e., date, time, location, and local magnitude *ML*) were taken from the IS-EPOS platform [42], which is part of the European Plate Observation System (EPOS) project. The boundaries of the mining area were taken from the Central Geological Database of the Polish Geological Institute. Appendix A shows interferometric phase band of study area in the period from 30 March 2017 to 05 April 2017.

Between 17 and 27 March 2017 in the area of interest, there were seven reported tremors of the rock mass of magnitude ML = 2.3–2.6. Table 1 gives the detailed data concerning the analyzed seismic events. The first two tremors (*ML* = 2.3 and 2.5) were recorded within 24 h of each other on 17 and 18 March 2017, respectively, at around 17:00 UTC. On the 22–25 March, there were four tremors (*ML* = 2.5 or 2.6). The last one in the discussed series was reported before 21:00 UTC on 27 March 2017 (*ML* = 2.6). The horizontal distance between any two epicenters did not exceed 75 m (Figure 3).


**Table 1.** Characteristics of mining-induced tremors.

**Figure 3.** Surface distribution of mining tremors epicenters from 17 March 2017 to 27 March 2017 over longwall mining operations in seam 2.

The surface impact of the seismic events was evaluated using differential InSAR (DInSAR) technology. For this purpose, 15 radargrams from the Copernicus Open Access Hub covering the period between 6 March and 29 May 2017 were taken. The mission satellites Sentinel-1 A and B (each with a revisit time of 12 days) enabled the acquisition of radargrams with a 6-days interval (pass Ascending, relative orbit 102). SLC data acquired in Interferometric Wide swath mode and a ground resolution of 5 m × 20 m (range × azimuth) were used. Table 2 shows the radar images used for the investigation.



SNAP (SeNtinel's Application Platform) and SNAPHU (Statistical-Cost, Network-Flow Algorithm for Phase Unwrapping) software dedicated for processing the satellite images, were used to processing radargrams. In the first step, co-registration of radargrams was performed. It consisted in spatial matching of the corresponding pixels of two processed master and slave radargrams acquired in 6-day intervals. This process took into account the data on the precise orbits of satellites and the DEM obtained from the 1" Shuttle Radar Topography Mission. In the next step, the phase difference for individual pixels was calculated, creating an interferogram. At the same time, a coherence value was calculated for each pair of pixels. During the generation of interferograms, the part of the phase resulting from the curvature of the earth was removed. It was possible thanks to the precise calculated coordinates of the satellites and the used DEM. Additionally, based on the DEM, the component responsible for the terrain topography was removed from the determined value of the phase difference. In order to reduce the noise of the radar signal, the interferogram was filtered by the Goldstein method. As a result, a differential interferogram was obtained, in which the phase difference of individual pixels indicated changes in the height of the land surface over a period of 6 days. Information on the value of the phase difference, ranging from 0 to 2π, was developed with the use of SNAPHU software (the so-called phase unwrapping). The developed wavelength was calculated to vertical displacement *s* by Equation (1).

$$s = \frac{\text{unuvrapped phase} \ast \text{unuv length}}{-4\pi \ast \cos(\text{incident angle})} \text{ [mm]} \tag{1}$$

In the last step, the obtained image of the altitude changes of the terrain surface was geo-referenced (EPSG: 2177).

#### *2.3. Determination of Subsidence Trough*

Figure 4 shows the stages of data processing in the Geographic Information System (GIS) environment. The final subsidence trough was digitized as a sum up all 6-days subsidence once.

Step A presents 6-days subsidence based on radargrams. In step B, the raster file with subsidence (Figure 4A) was automatically converted into a vector image of subsidence isolines (Figure 4B). However, the isolines so generated do not represent fully the actual spatial image of the effects of the analyzed deformation, the reason being the residual signal noise (e.g., low coherence, plant vegetation, ground moisture, weather conditions), which could not be eliminated fully in the radargram processing. During the three-dimensional (3D) modeling of the surface deformation, attempts were made to eliminate the occurring irregularities of the isoline distribution of surface subsidence (Figure 4B; blue lines) and the subsidence profile (Figure 5) caused by the radar signal noise. This modeling was in two steps: first, by interpolating the subsidence profile using eighth-degree orthogonal polynomials (red line, Figure 5). At the intersection of the P-P' profile and the approximated subsidence, marked as green points (Figures 4C and 5), corrections for individual subsidence isolines [39] were determined as the absolute value of the difference of subsidence Θ, as follows in Equation (2):

$$\Theta = |\text{ s}\_{\text{RASTER}} \text{ - spOL} \,\text{y}\,\, | \,\tag{2}$$

where SRASTER is subsidence based on raster data, SPOLY is subsidence based on orthogonal polynomials approximation.

**Figure 4.** Steps of determining 6-day subsidence isolines on the basis of a raster image and P-profile: (**A**) vertical displacements based on interferometric data, (**B**) automatic vectorization of land surface subsidence isolines, (**C**) approximation of subsidence profiles along the P-profile and smoothing of subsidence isolines, (**D**) surface distribution of corrected subsidence isolines.

**Figure 5.** The approximation of the subsidence profile (SPOLY, red line) and the subsidence profile based on raster data (SRASTER, light-blue line).

In the second step manually entering the corrections to the surface distribution of isolines for each individual subsidence once by manual digitalization (Figure 4C; red lines). The presented procedure may be used for several field profiles. Interpolating the profiles allowed correct assignment of subsidence values for the corrected isolines. This provided a spatial image of the deformations (3D model) close to the actual conditions (Figure 4D).

To obtain an adjustment of an orthogonal polynomial, the mean-square criterion for minimizing the distances between the approximating function and the discrete data set was used. The main problem in such cases is determining the optimum degree of the approximation orthogonal polynomial. Hejmanowski and Kwinta [43] proposed using the following criterion from Equation (3):

$$p = \min\_{j=0,1,\ldots,n} (j) \; \land \; m\_{\mathfrak{s}}(j+1) - m\_{\mathfrak{s}}(j) \le M\_{\mathfrak{s}} \cdot m\_{\mathfrak{s}}(0) \tag{3}$$

where *j* is the optimum degree of the approximating orthogonal polynomial, and *ms*(*j*) is the mean error in approximating subsidence *s* with the orthogonal polynomial of degree *j*; *Ms* is a coefficient that describes the stochastic variability of the subsidence Equation (4).

Polish measurements suggest that *Ms* is between ±0.01 and ±0.04; here, *Ms* = 0.01 is assumed.

$$M\_s = \frac{\sigma\_s}{|\mathcal{S}\_{\max}|} \tag{4}$$

where σ*<sup>s</sup>* is a standard deviation of observed subsidence average course, and *sˆmax* is maximum subsidence of observed subsidence average course.

The analysis of fitting the subsidence profiles using the proposed criterion indicated that (I) the loss function decreases as the polynomial degree is increased; (II) after a certain degree of the orthogonal polynomial, the approximation error stabilizes and consequently the graphs of approximating functions have a similar course; and (III) in the case of surface subsidence, the optimum orthogonal polynomial degree for different sets of input data (vertical displacements observed in different locations) may not be the same. Based on the performed analyses [43,44], it can be assumed that the optimum orthogonal polynomial degree *p* for surface subsidence is between 6 and 8. Herein, *p* = 8 was assumed because of the iterative nature of the process of reaching the final model of vertical displacements (Figures 9 and 10).

#### **3. Results**

Based on the available radar data and their processing according to the presented methodology, a series of 6-days vertical displacements was obtained for 6 March to 29 May 2017. In the case of the last two interferograms for 17 May to 23 May 2017 and 23 May to 29 May 2017, insufficient coherence (due mainly to the intense plant vegetation) made it impossible to determine the vertical displacements to a sufficiently high reliability. Connected to this fact, some data were lost for the end of the planned study period, this being due directly to the limitations of the DInSAR method. Figure 6 shows the obtained 6-days land surface deformations between 6 March and 17 May 2017.

In the subsequent stage of data processing, the vertical displacements between 6 March and 17 May 2017 were summed up. Six-days subsidence isolines were also generated (starting on 12 March 2017) to obtain the incremental values. Consequently, this allowed field profiles to be developed, particularly ones crossing the point at which maximum surface subsidence occurred (Figure 7, yellow point). On this basis, it was possible to present the development of surface deformation over time. Figure 8 presents an example of the field profile (in the SW-NE direction).

Based on the summed-up raster data (Figure 9A) and abovementioned methodology of data processing (Section 2.3) the total vertical displacements were generated (Figure 9B). Final subsidence trough with 50-mm isolines is shown in Figure 9D. According to the prepared methodology, an example of field profiles is presented in Figure 10. Combining the information from interpolation (two-dimensional) with the digitized surface distribution of subsidence enabled generating a spatial model (3D) of the surface deformation caused by the series of rock mass tremors (Figure 11).

**Figure 6.** Days subsidence maps in the period 6 March to 17 May 2017 with registered epicenters of mining-induced tremors.

**Figure 7.** The final subsidence trough in the period from 6 March 2017 to 17 May 2017 with profile lines (black lines). The SW-NE red line represents the main profile line.

**Figure 8.** Development of subsidence trough in the period from 6 March 2017 to 17 May 2017.

**Figure 9.** Results of final land surface subsidence modeling based on methodology described in Section 2.3. (**A**) summed-up vertical displacements from 6 March 2017 to 17 May 2017, (**B**) automatic vectorization of land surface subsidence isolines, (**C**) approximation of subsidence and smoothing of subsidence isolines, (**D**) surface distribution of corrected subsidence isolines.

**Figure 10.** Orthogonal polynomial interpolation of incremental subsidence.

**Figure 11.** 3D model of final subsidence trough.

#### **4. Discussion**

The observed phenomenon of surface transformation was a consequence of deformations in the rock mass. It was characterized by dynamic conditions varying in time, related to the changing speed of incremental surface subsidence (Figure 12). The authors' experience [45–47] shows that the impact of drainage of aquifers in the conditions of Polish coal mining is small and does not exceed 3–4 mm in the analyzed period. For this reason, dehydration subsidence was not considered in this article. In this case, the impulse that initiated the deformation process was a series of seven tremors that occurred within 10 days of each other between 17 and 27 March 2017. The dynamics of the variations in the land surface morphology were characterized by analyzing the subsidence profiles in time for a point in the center of the resulting subsiding trough (Figure 7, yellow point). The initial speed of land subsidence caused by the moving rock mass material toward the cavity was approx. −3 mm/d. A few days after the last tremor (on 27 March 2017), the speed started to increase nonlinearly. In the middle of April, the subsidence speed reached its maximum value of approximately −14 mm/d. In the subsequent period, a slowing decreasing rate was reported. In the middle of May, it was between −5 and −6 mm/d, indicating slow fading of the surface displacements.

This phenomenon can be explained by the delay related to the time required for deformations to pass through the rock mass on the way from the cavern to the land surface [48,49]. This can be compared to the propagation of electromagnetic waves in a damping medium. Changing speeds during the process of rock mass deformation indicate similarity to the exponential models of natural processes in time, in which a damping factor is present (e.g., mechanical systems with motion resistance; oscillating systems with energy loss; charging and discharging electrical systems).

The relatively low quality of radar data in the following period (17–29 May 2017) made it impossible to determine exactly the date on which the land surface subsidence ceased. This was caused by the high intensity of plant vegetation after 17 May 2017, which reduced the legibility of the acquired data to a level below the threshold for interpreting the processed DInSAR images. This movement could have continued for between two and four weeks, increasing the subsidence of the trough by no more than between −5 and −10 cm.

**Figure 12.** Subsidence and speed of surface subsidence initiated by mining-induced tremors.

Strictly speaking, the land surface deformation presented herein does not show features that are characteristic of the effects of the ongoing mining operations [2,50]. This is indicated by the generated maps of vertical displacement (Figure 6), presenting continuous data about the subsidence that occurred between 6 March and 17 May 2017. The time-based subsidence (6-day) detected in the subsequent images covered almost the same surface area all the time (800 m × 700 m) and are strongly correlated with the epicenters of the tremors that preceded the examined phenomenon. In the case of progressive exploitation of the hard coal longwall mining, the increasing size of the voids formed in the rock mass would cause subsidence on the surface that covers more and more area. The indicated temporary subsidence would then cover an increasing area on the surface. Subsidence development would follow the direction of the longwall face. In the conditions of the USCB, increments in subsidence can reach values higher than −20 mm/d [51]. Additionally, the expected land subsidence amounts to about 70–80% of the thickness of the extracted coal seam, and 95% of this subsidence is revealed in the period from 3 to 9 months [52]. In the case of the analysis of subsidence accumulation dynamics, a significant similarity can be shown between the development of post-tremor and mining-induced subsiding trough. This is related to the physical and mechanical parameters of the rock medium in which deformation occurs. In both the first case (post-tremor, Figure 12) and the second case (mining-induced), three basic phases of deformation development can be distinguished: (1) The initial phase (occurring in early April) involves the first reaction of the rock mass on the cause of subsidence, in which the subsidence speed began to rise; (2) the main phase (ending in the third decade of April), in which the subsidence speed stabilized at a specific and usually relatively high level (the maximum speed occurs in this phase); (3) the ending phase (start of May and later), where the impact dynamics abates and as a result the subsidence speed is reduced to zero.

All the observed phases of deformation development, particularly the few weeks with symptoms of deformation showing on the surface, are contrary to the earlier experiences resulting, for example, from observations performed in previous years in the LGCD [32,33]. Approximately 80% of the total subsidence values can be detected within 48 h from the occurrence of the rock mass tremor. The final subsidence *smax* is between −8 and −12 cm (speed of subsidence: *s˙* = 20–35 mm/d). The shapes of the formed subsiding troughs are close to an ellipse with dimensions of approximately 1.5 km × 2.5 km. Research conducted so far in the USCB [29,31] refers to the detection of subsidence caused by seismic phenomena and the distances from the points of maximum subsidence to the tremor epicenters, but it does not cover the aspect of dynamics (time) of detecting the surface impacts.

Sroka [53–55] presented the dynamics of the deformation process in German and Polish mining operations. That work shows that significant damage of buildings and structures in the mining areas is caused, among other factors, by the high and variable speeds of deformation of the rock mass and land surface. Factors that determine a threat to buildings and structures under the influence of dynamics of the deformation process include the speed of land surface subsidence *s˙*. Table 3 gives the upper limit *s˙*<sup>B</sup> of the boundary speed of surface subsidence along with the corresponding categories of building strength and mining areas in Poland.


**Table 3.** Boundary speed of surface subsidence s˙B in relation to building strength categories and category of mining areas.

The incremental subsidence speed was analyzed based on the similarity of the dynamics of mining-induced and post-tremor subsidence processes. The dynamics of this process (Figure 12) indicate that the process of gradual increase in subsidence complies with the exponential model of the dynamics of other natural processes in which damping effects can be found. Such a model exists in the description of the land surface subsidence process in time, as a result of mining operations. Equation (5) present the formula for the course of subsidence in time at a given point [56]:

$$s\_l = s \begin{pmatrix} 1 - c^{-ct} \end{pmatrix} \tag{5}$$

where *s* is the final subsidence as a result of mining operations, *st* is the subsidence at time *t*, and *c* is a parameter describing the damping effect (delay in reaching the final subsidence at a given point).

The strength of buildings to the dynamic impact of mining exploitation is characterized by the subsidence speed (s˙B). The maximum subsidence speed (Figure 12) where the subsiding trough forms is s˙max = −13.8 mm/d. This is a value that poses a damage hazard for infrastructure that is qualified as resistant to the 3rd category of mining impact. The strength category of buildings and structures is used in the context of mining area protection. Building structures are assigned to a category of resistance against mining impacts [57,58] described by the values of the so-called deformation factors, e.g., changes in the area tilt T and horizontal deformation ε (relative change of the distance between observation points) [59,60]. This indicates that the building construction can transfer the deformation impact caused by mining operations. In the five-stage scale of strength categories, buildings classified to 3rd category are resistant to significant land surface deformations (−6 mm/d < s˙B ≤ −12 mm/d [55], 5 mm/m < TB ≤10 mm/m, 3 mm/m < |εB| ≤ 6 mm/m [59]). Therefore, it can be concluded that buildings (with strength category 0–3) exposed to a dynamic impact whose magnitude is comparable to that described here, depending on their strength, could be damaged. However, because of the lack of such examples, it is not possible to verify this hypothesis, but nevertheless attention should be paid to that type of threat.

#### **5. Conclusions**

Rock mass tremors and their related effects on urbanized areas are definitely unfavorable phenomena and can threaten buildings, infrastructure, and people. In the case of tremors induced in areas where rock mass voids exist (natural caverns, closed workings, fault zones), apart from the surface tremors caused by the course of a seismic wave, significant morphological changes of this surface should also be expected. These changes can be observed for example as continuous deformations, such as subsiding troughs. The general availability of seismic data (e.g., GEOFON program from GFZ German Research Centre for Geosciences, IS-EPOS, European-Mediterranean Seismological Centre) and satellite imaging (Sentinel mission from European Space Agency) and also the possibilities for data processing in software licensed under the GNU General Public License allow examination of large-scale area deformations. They also allow the assessment of the level of risk to the area covered by such deformations. An excellent example of such use of data is the case analyzed herein. Assessing the impact of rock mass tremors on the land surface required developing a bespoke method for processing the satellite data, leading to the obtained model of a post-tremor subsiding trough. As a result, the characteristics and dynamics of the surface deformation process caused by a series of seven tremors of the rock mass with magnitudes reaching *ML* = 2.6 were determined in the studied region. The tremors occurred between 17 and 27 March 2017. Moreover, based on experience in protecting mining areas, the threats to building structures were evaluated according to the classification used in mining areas. That evaluation was done to indicate potential threats for the developed area on the surface. In the discussed case, building structures with resistance lower than the limit values of subsidence speed (*s˙B* = −12 mm/d; Table 3) can be exposed to the risk of damage.

Because of the significant scale of occurrence of seismic activities caused by humans, for example, related to the excavation of mineral deposits, it is important to examine the effects of such phenomena

regarding public safety. To date, this problem has not been identified sufficiently despite the large number of seismic events that occur, such as those in the USCB. This creates favorable conditions for performing further studies in this field.

**Author Contributions:** Conceptualization, P.S., T.S., A.W. and D.M.; methodology, P.S., T.S., A.W. and D.M.; software, A.W. and D.M.; validation, T.S. and A.W.; formal analysis, P.S. and T.S.; investigation, D.M.; resources, A.W. and D.M.; data curation, A.W. and D.M.; writing—original draft preparation, T.S. and D.M.; writing—review and editing, P.S. and A.W.; visualization, A.W. and D.M.; supervision, P.S. and T.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the research subvention of AGH University of Science and Technology, grant number 16.16.150.545.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Figure A1.** Interferometric phase band of study area in the period from 30 March 2017 to 05 April 2017.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Workflow Based on SNAP–StaMPS Open-Source Tools and GNSS Data for PSI-Based Ground Deformation Using Dual-Orbit Sentinel-1 Data: Accuracy Assessment with Error Propagation Analysis**

**Francesco Mancini, Francesca Grassi and Nicola Cenni**



**Citation:** Mancini, F.; Grassi, F.; Cenni, N. A Workflow Based on SNAP–StaMPS Open-Source Tools and GNSS Data for PSI-Based Ground Deformation Using Dual-Orbit Sentinel-1 Data: Accuracy Assessment with Error Propagation Analysis. *Remote Sens.* **2021**, *13*, 753. https://doi.org/10.3390/rs13040753

Academic Editor: João Catalão Fernandes

Received: 11 January 2021 Accepted: 13 February 2021 Published: 18 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Abstract:** This paper discusses a full interferometry processing chain based on dual-orbit Sentinel-1A and Sentinel-1B (S1) synthetic aperture radar data and a combination of open-source routines from the Sentinel Application Platform (SNAP), Stanford Method for Persistent Scatterers (StaMPS), and additional routines introduced by the authors. These are used to provide vertical and East-West horizontal velocity maps over a study area in the south-western sector of the Po Plain (Italy) where land subsidence is recognized. The processing of long time series of displacements from a cluster of continuous global navigation satellite system stations is used to provide a global reference frame for line-of-sight–projected velocities and to validate velocity maps after the decomposition analysis. We thus introduce the main theoretical aspects related to error propagation analysis for the proposed methodology and provide the level of uncertainty of the validation analysis at relevant points. The combined SNAP–StaMPS workflow is shown to be a reliable tool for S1 data processing. Based on the validation procedure, the workflow allows decomposed velocity maps to be obtained with an accuracy of 2 mm/yr with expected uncertainty levels lower than 2 mm/yr. Slant-oriented and decomposed velocity maps provide new insights into the ground deformation phenomena that affect the study area arising from a combination of natural and anthropogenic sources.

**Keywords:** Interferometry; PSI; SNAP-StaMPS; Sentinel-1; Ground deformation

#### **1. Introduction**

Satellite radar interferometry is recognized as an effective technique in different applications focused on deformation phenomena that occur on the Earth surface. The successful application of this methodology is strictly related to the ability to depict the displacement of ground targets at a very high level of accuracy (1–2 mm/yr) using synthetic aperture radar (SAR) images and time series analysis of displacements oriented along the Line Of Sight (LOS) directions. Since the first applications in the early 1990s, with the first generation of SAR satellites (ERS-1 and Radarsat-1), the temporal and spatial resolution increased, and a number of research fields have taken advantage of satellite radar interferometry. Among these, applications to natural hazard assessment, ground deformation, and investigations into polar cap dynamics, slope failures, and instability processes are relevant to the present study [1–7].

Differential SAR interferometry (InSAR) was initially used to measure deformation of the land surface through interpretation of interferograms, combined with digital elevation models (DEM) to remove the topographic contributions. DEM were obtained with the SAR interferometry by using one or two additional SAR images and it was called three, or four, pass interferometry, respectively. Today, approaches based on the stacking of interferograms are adopted more frequently, and the methodologies can be differentiated

into persistent scatterer interferometry (PSI), small baseline subsets (SBAS), and methods based on the search of distributed scatterers rather than the persistent ones (Homogeneous Distributed Scatterer Interferometry, HDSI). These PSI [8–11], SBAS [12,13], and HDSI [14] methods require different strategies for the generation of the interferogram stacks and the statistical approach to find the velocity of stable and highly coherent scatterers from the time series of displacement between satellite passes. Regardless of the methodology used, one of the main advantages with respect to geodetic terrestrial surveys is that it is possible to find the displacement of a multitude of points with significant improvement in spatial density, and better delineation of the displacement phenomena at the required spatial scale. These methodologies are based on differential approaches, and the displacement derived is referred to a reference position. Thus, the absolute velocity of the scatterers is required whenever the relative displacement field must be transposed into an absolute reference frame. Major limitations can be associated with spatial and temporal decorrelation phenomena, image gaps, the impossibility to detect fast motions due to the used wavelengths (typically the radar bands C and X), and the ability to detect a minor part of the actual displacements due to the sampling along the slant direction. In particular, also depending on the acquisition geometry, the methodologies are more sensitive towards displacement that occurs along the vertical and East-West directions. To obtain the final velocity field in a geodetic reference frame, displacements along the vertical and horizontal directions are required and a vector decomposition analysis must be applied. However, during the decomposition from ascending and descending slant velocities to different directions the uncertainties propagate. See [15] for details about the application of the error propagation law to geodesic and surveying science.

A number of SAR images provided by ended (e.g., European Remote Sensing [ERS]1– 2, Environmental Satellite [Envisat], Radarsat 1, ALOS PALSAR), and active (e.g., Radarsat 2, TerraSAR-X, ALOS-2, COSMO SkyMed) satellite radar missions have been used for deformation analyses through the InSAR approaches, with some remarkable results obtained. In the framework of the Copernicus Earth Observation programme, Sentinel-1A and Sentinel-1B (S1A, S1B) satellites equipped with radar sensors were launched in 2014 and 2016, respectively, as part of the Sentinel constellation of satellites. Under the umbrella of Copernicus, the European Space Agency (ESA) started to provide worldwide free-access SAR images at average spatial resolution and with very short revisiting times, which opened new perspectives for continuous ground-surface monitoring [16]. The time series of the S1 SAR images is now long enough to process reliable displacement maps with more reliable detection of seasonal and/or long-term trends. Moreover, the extension of available time series makes the validation procedure based on comparisons with external data more reliable. The studies focused on ground deformations phenomena from the processing of S1 SAR image could be based on single or dual-orbit S1 data. These studies depicted ground deformation maps of large areas [17–22] with no further decomposition analysis and validation of average vertical and East-West oriented displacements by comparison with external data complemented by uncertainty analyses.

In addition, users now have the opportunity to process SAR images using trusted and freely accessible toolboxes. Over the last years, the ESA has distributed the Sentinel Application Platform (SNAP) software with incorporated utilities for interferogram generation and stacking, with this complemented by the Statistical Cost Network Flow Algorithm for Phase Unwrapping (SNAPHU) package, for phase unwrapping [23–25]. The displacement histories of persistent scatterers can be derived from the stacks of interferograms using the free Stanford Method for Persistent Scatterers (StaMPS) tool, which only requires interferograms from the SNAP since it uses an unwrapping algorithm derived from SNAPHU [26]. A few studies have explored the combined use of SNAP and StaMPS for displacement analysis by the PSI method from ascending and descending orbits [18,19,27] but a validation procedure complemented by uncertainty analysis has rarely been carried out. Manunta et al. [20] used accurate continuous global navigation satellite system (CGNSS) positions to refer the displacements for the whole Italian territory to a global reference frame, as

obtained from descending S1 data and the SBAS method. Delgado Blasco et al. [19] developed a workflow based on the SNAP and StaMPS tools to process dual-orbit S1 data with successive vector decomposition, to obtain the actual vertical motion component. They used the velocities available from CGNSS sites as reference, with no further validation steps. To date, the scientific literature does not provide evidence on the accuracy assessment of ground deformation maps based on dual-orbit S1 data and the SNAP–StaMPS workflow with respect to a defined global reference frame. However, a similar approach was described by [28,29].

The present study introduces the strategies adopted in the processing of S1 dual-orbit data using the SNAP and PSI–StaMPS open-source routines with constraint of single orbit products to a global reference frame, decomposition analysis and accuracy assessment at validation sites complemented by an error propagation analysis. To align PSI slant velocities to a reference frame and validate results on selected sites, we processed GNSS observation from continuous stations distributed over the study area. In particular, we discuss the performances of the SNAP–StaMPS workflow for the investigation of ground deformation over metropolitan areas located in the southeastern border of the Po Plain (Italy). This represents a unique case study where significant ground deformation phenomena have been reported over a large area due to the combination of natural (e.g., quaternary sediment compaction of the Plio-Quaternary deposits and deep tectonics) and anthropogenic (e.g., mainly pumping of groundwater for industrial and drinking purposes or deep gas field exploitation) factors [30,31]. Subsidence phenomena have been previously investigated through analyses based on PSI from data provided by both single and dual ERS and Envisat satellite orbits [29,32,33]. However, no updates have been provided using data from dual orbits modern satellite radar sensors and CGNSS data. To cover this gap, we processed the full range of the available ascending and descending S1 data, along with the observations provided by the CGNSS permanent stations installed in the area. We retrieved the vertical and East-West oriented velocity fields after vector decomposition and validate results following the above discussed approach. We obtained a ground deformation map which updates the present knowledge about phenomena occurring in the investigated area and discussed a few case studies where vertical and horizontal displacements can be linked to different factors.

#### **2. Study Area**

This paper investigates an area of the eastern sector of the alluvial Po plain (Figure 1, inset). This area is characterized by the historical and present lowering of the ground soil due to a combination of anthropogenic and long-term geological processes. The human-induced contribution is mostly connected to the demand for the groundwater from a well-developed multi-aquifer system that started during the second half of the 20th century, in parallel with the increasing industrial activities [34,35]. A further minor human contribution has come from the extensive agricultural and zootechnical techniques and from gas exploitation [36], the latter strongly limited to the extents of productive areas. Indeed, the loading and compaction of the Holocene sediments is the main source of subsidence from natural processes in the alluvial plain, in addition to minor contributes due to deep active tectonics. As reported by [37], the modern rate of subsidence is at least an order of magnitude greater than the historical rate. The major cities included within the study area are Bologna and Ferrara, settled on areas that are characterized by very different geological settings and superficial deformation processes [37] (Figure 1).

**Figure 1.** Map of the study area, including the main urbanized areas and seismogenic faults (dotted lines, SHARE project, [38]. Eleven continuous global navigation satellite system (CGNSS) sites used in the present work as constraint to a reference frame and validation purposes are located. Labels, composed of four alphanumeric codes, are printed above the CGNSS locations.

#### **3. Interferometric Data Processing and Analysis of GNSS Observations**

In this section, main theoretical foundations for the SAR interferometry and the workflow adopted to process the S1 SAR images and observations from the CGNSS sites are discussed.

#### *3.1. Interferogram Generation*

To produce the stack of interferograms, the whole dataset of ascending and descending S1 SAR images were processed with respect to a single master image. The selection of the master image was based on minimization of the possible geometrical and temporal decorrelation effects [7,8,39]. The selection of the master image might thus have a significant role in the effectiveness of the whole procedure. The joint effects of sources of decorrelation have to be considered, with estimation of the total correlation coefficient [40],

$$\rho\_{\text{total}} = \rho\_{\text{temp}} \cdot \rho\_{\text{spat}} \cdot \rho\_{\text{Doppler}} \cdot \rho\_{\text{therm}} \approx \left(1 - f\left(\frac{T}{T\_{\text{C}}}\right)\right) \cdot \left(1 - f\left(\frac{B\_{\text{PEEP}}}{B\_{\text{PEEP,C}}}\right)\right) \cdot \left(1 - f\left(\frac{F\_{\text{DC}}}{F\_{\text{DC,C}}}\right)\right) \cdot \rho\_{\text{therm}} \tag{1}$$

where:

$$f(\mathbf{x}) = \begin{cases} \mathbf{x}, \text{ for } \mathbf{x} \le 1\\ 1, \text{ for } \mathbf{x} > 1 \end{cases}$$

and *ρtemp* is the temporal correlation coefficient, *ρspat* is the spatial correlation coefficient, *ρDoppler* is the Doppler correlation coefficient, *ρtherm* is the thermal correlation coefficient, *BPERP* is the perpendicular baseline, *T* is the temporal baseline, *FDC* the Doppler centre frequency difference between two images, and the term C is the critical one. *BPERP*,*<sup>C</sup>* and *FDC*,*<sup>C</sup>* terms do not represent limiting factors in S1 data processing being the S1 mission designed to minimize these sources of decorrelation. Then, assuming *ρtherm* is constant, the master is chosen where the value of ∑*<sup>N</sup> <sup>i</sup>*=<sup>1</sup> *ρtotal* is maximized, where *N* is the total number of images.

The interferometric processing of data acquired using the Terrain Observation with Progressive Scan (TOPS) mode follows the approach presented by [41,42]. The interfero-

metric operations are performed at a burst level, and thus the selection of the same swaths and bursts for the master and slave images is required. Then, the coregistration of the SAR images is performed by exploiting the precise orbit information (provided after two weeks from the acquisition time) and an external digital elevation model (DEM) of suitable accuracy and spatial resolution. After coregistration of the image pairs, the merging of adjacent bursts in azimuth direction (i.e., the deburst step) and the selection of the study area can be performed. The image pairs are then used in the generation of the stack of interferograms. Finally, the contribution of the topography to the interferometric phase must be removed, with the support of the external DEM.

#### *3.2. PSInSAR*

It is worth noting that the developed workflow for the selection of persistent scatterers and the computation of their displacement history follows the theory of the persistent scatterer InSAR technique. The basic method to identify and select persistent scatterers considers a constant velocity model for targets motion [9,39]. Once these pixels are selected, their phase history is analyzed, and only the pixels with a history similar to the functional model are included as final persistent scatterer candidates [9]. This approach can limit the application of the technique when the actual deformation model differs from that stated in the initial hypothesis. According to the StaMPS method, initial selection of the persistent scatterer candidates is performed based on amplitude analysis [11]. First, by exploiting the statistical relationship between the amplitude and phase stabilities, a high value of the amplitude dispersion as defined by [9] is used as a threshold value:

$$D\_A = \frac{\sigma\_A}{\mu\_A} \tag{2}$$

where *σ<sup>A</sup>* is the standard deviation, and *μ<sup>A</sup>* is the mean of the amplitude of the pixel; a value of 0.4 was suggested by [10] for such parameter. Then the phase stability of each pixel is analyzed by estimating the phase noise contribution [11] through an iterative process and checking that it does not obscure the signal and, in particular, the displacement term. Indeed, the wrapped phase of the x-th pixel and the i-th interferogram is given by the following equation:

$$\psi\_{\mathbf{x},i} = \mathcal{W}\{\phi\_{D,\mathbf{x},i} + \phi\_{A,\mathbf{x},i} + \Delta\phi\_{S,\mathbf{x},i} + \Delta\phi\_{\theta,\mathbf{x},i} + \phi\_{N,\mathbf{x},i}\}\tag{3}$$

where *W*{·} is the wrapping operator, *φD*,*x*,*<sup>i</sup>* is the phase change due to displacement along the line of sight (LOS) of the pixel, *φA*,*x*,*<sup>i</sup>* is the phase contribution due to atmospheric refraction, Δ*φS*,*x*,*<sup>i</sup>* is the residual phase that depends on the satellite orbit inaccuracies, Δ*φθ*,*x*,*<sup>i</sup>* is the residual phase due to look angle error, and *φN*,*x*,*<sup>i</sup>* is the phase noise term [11].

The final selection of eligible persistent scatterers is performed through a statistical approach. In this step, pixels that do not have stable behavior along the period spanned by the data and/or are affected by the response from neighbouring pixels can be removed. Indeed, even if scatterers remain dominant in a specific pixel, as its characteristics can change over time, so can its response to radar signals. Such a pixel would not be considered a persistent one, and it must be discarded. Moreover, pixels adjacent to a persistent scatterer are considered to belong to the same scatterer, and the physical location corresponds to the pixel with the highest *γ<sup>x</sup>* (see Equation (20) in [11]). The spatially uncorrelated part of the signal is primarily due to the spatially uncorrelated part of the look angle error (Δ*φ<sup>u</sup> θ* ) and the contribution of the master (*φ*ˆ*<sup>x</sup> <sup>m</sup>*,*u*). This part is estimated and subtracted from the unwrapped phase. This step is mandatory, because the spatially uncorrelated contribution can affect the wrapped phase, and thus the success of the whole unwrapping process. The wrapped phase can be written as:

$$\mathcal{W}\{\psi\_{\mathbf{x},i} - \Delta\phi^{\mu}\_{\theta} - \hat{\phi}\_{\mathbf{x}}^{\ \ m,\mu}\}\tag{4}$$

Then, the unwrapping, performed by the method discussed in [26], yields the unwrapped phase:

$$
\hat{\phi}\_{\mathbf{x},i} = \phi\_{\mathbf{D},\mathbf{x},i} + \phi\_{A,\mathbf{x},i} + \Delta\phi\_{S,\mathbf{x},i} + \Delta\phi^{\mathbf{c}}\_{\theta,\mathbf{x},i} + \phi\_{N,\mathbf{x},i} - \hat{\phi}\_{\mathbf{x}}{}^{m,\mu} + 2k\_{\mathbf{x},j}\pi \tag{5}
$$

where Δ*φ<sup>c</sup> <sup>θ</sup>*,*x*,*<sup>i</sup>* is the spatially correlated part of Δ*φθ*,*x*,*i*, and *kx*,*<sup>i</sup>* is the remaining unknown integer ambiguity. The spatially correlated phase terms are estimated using temporal and spatial filtering and subtracted in order to retrieve the displacement term. Finally, the atmospheric phase can be removed.

#### *3.3. Atmospheric Filtering*

For the estimation of the tropospheric phase contribution, we used the linear phasebased model included in the Toolbox for Reducing Atmospheric InSAR Noise (TRAIN), a module developed by the University of Leeds to reduce the atmospheric noise [43,44]:

$$
\Delta\phi\_{tropo} = K\_{\Delta\phi}h + \Delta\phi\_0 \tag{6}
$$

where Δ*φtropo* is the interferometric tropospheric phase contribution, *h* is the topography, *K*Δ*<sup>φ</sup>* is the coefficient that relates the phase to the topography, and Δ*φ*<sup>0</sup> is a constant phase shift that is applied to the extent of the interferogram. The atmospheric filtering is performed in a final step of the PSI processing to remove residual and local artifact.

#### *3.4. Vector Decomposition with Uncertainties*

The LOS-oriented displacements, detected from the interferometric processing of the ascending and descending orbits, can be decomposed into the actual displacement along the North, East, and vertical directions using the procedure introduced in [45]. In general, for a satellite orbit with heading *α<sup>H</sup>* the mean velocity along the LOS can be expressed following the theory of [46], where instead of the displacement, the mean velocity of the persistent scatterer is considered, without losing the validity of the discussion, as shown in the equation:

$$v\_{LOS} = v\_{Up}\cos(\theta\_{inc}) - \sin(\theta\_{inc})\left(v\_N\cos\left(a\_H - \frac{3\pi}{2}\right) + v\_E\sin\left(a\_H - \frac{3\pi}{2}\right)\right) \tag{7}$$

where *θinc* is the incidence angle, *vUP* is the mean velocity component in the vertical direction, *vN* is the mean velocity component in the North direction, and *vE* is the mean velocity component in the East direction. Equation (7) can be adopted to project absolute velocities provided by the processing of CGNSS time series into ascending and descending slant range directions. In the present work, this procedure is used to constrain and align slant-oriented velocities to a global reference frame. This is performed prior of further decomposition analyses. In addition, from Equation (7), the propagation of uncertainties that affect the CGNSS site velocities along the North, East and vertical (up) directions can be performed to obtain the level of uncertainty of the LOS-oriented velocities:

$$\sigma\_{V\_{\rm LOS}}^2 = \cos^2 \theta\_{\rm INC} \cdot \sigma\_{V\_{\rm IDP}}^2 + \sin^2 \theta\_{\rm INC} \left[ \cos^2 \left( a\_\pi - \frac{3}{2} \pi \right) \cdot \sigma\_{V\_{\rm IN}}^2 + \sin^2 \left( a\_\pi - \frac{3}{2} \pi \right) \cdot \sigma\_{V\_{\rm E}}^2 \right] \tag{8}$$

Following the theory proposed by [46], the sensitivity to displacement oriented along the North-South direction is very low. The S1 orbit is designed with a mean incidence angle of approximately 37◦ and a heading angle of 190◦. Given the acquisition geometry, the sampling of the actual displacement oriented along the North-South direction by the slant looking shows a sensitivity of −0.10. Consequently, the vector decomposition procedure aims to compute the displacement along the Earth-West (E) and vertical (Up) components from the ascending and descending passes only. Using the matrix notation, the expression given in the following equation holds:

$$
\begin{bmatrix} \upsilon\_{UP} \\ \upsilon\_{EAST} \end{bmatrix} = A^{-1} \cdot \begin{bmatrix} \upsilon\_{LOS}^{ASC} \\ \upsilon\_{LOS}^{BESS} \end{bmatrix} \tag{9}
$$

where:

$$A = \begin{bmatrix} \cos(\theta\_{ASC}) & \sin(\theta\_{ASC})\cos(\mathfrak{a}\_{H,ASC}) \\ \cos(\theta\_{DES}) & -\sin(\theta\_{DES})\cos(\mathfrak{a}\_{H,DES}) \end{bmatrix}.$$

Note that *θASC* has a negative value because of its counterclockwise counting. From Equation (9), the variance–covariance matrix of decomposed velocities can be expressed through the uncertainty error propagation law as:

$$\mathbb{C}\_{V\_{decup}} = A^{-1} \mathbb{C}\_{V\_{LOS}} \left( A^{-1} \right)^{T} \tag{10}$$

Equation (10) allows the computation of uncertainty level to be assigned at velocities along the Earth-West (E) and vertical (Up) directions after the decomposition of slantoriented velocities.

#### *3.5. Accuracy Assessment with Error Propagation Analysis*

To provide an accuracy assessment of PSI-derived velocities a comparison with CGNSS time series can be used. This comparison could be performed at the level of time series projected along a common direction or by using decomposed velocities. The comparison of average annual velocities from regression analysis applied to PSI and CGNSS time series could be poorly sensitive to unmodelled effects. However, its significance increases when comparing linear trends of slow ground deformation phenomena.

The approach followed in the present paper uses the error propagation law. As reported in Section 3.4, a level of uncertainty can be assigned to the decomposed PSI velocities and used to establish the statistical significance of the comparison with CGNSS actual velocities. This approach requires: (a) the assessment of uncertainties introduced by the alignment of PSI-derived slant velocities to the reference frame and (b) the error propagation analysis to obtain uncertainties of decomposed PSI velocities in the vicinity of a CGNSS sites used in the error assessment procedure. The decomposition procedure requires a preliminary sampling and averaging of PSI-based velocities over regular grid for ascending and descending products. Then, during the step (b), uncertainties linked to averaging procedure of slant velocities and decomposition procedure (the latter using Equation (10)) propagate. The comparison between PSI and CGNSS velocities could be therefore complemented by the uncertainty level.

#### *3.6. GNSS Data Processing Strategy*

The absolute reference frame for the differential displacement provided by the PSI methodology was defined using a network of CGNSS sites that are unevenly distributed within the study area. Moreover, CGNSS time series can provide fundamental information about the tectonic and geodynamic behaviors of the region and help in the understanding of the superficial kinematic phenomena. Unfortunately, CGNSS sites belonging to international networks are usually a limited number over an area of interest and stations from other positioning services could be of interest in the validation processes. Figure 1 shows the distribution of the CGNSS stations used in this work to provide a constraint to the global reference frame and reference absolute velocities at sites used in the assessment of PSI-derived velocity. These CGNSS sites are managed by different public and private agencies for different purposes, and for this reason, they are not evenly distributed, and the observation time spans might differ. Some sites were developed for scientific investigations by public research authorities, while others make up part of the national services for real-time positioning. However, [47,48] demonstrated that CGNSS sites designed and

realized for technical operations do not introduce noise into the data recorded, with respect to sites used for scientific applications.

In the present study, we have considered only the CGNSS sites located in the study area with observations acquired between 1 March 2015, and the end of 2019 were processed using the GAMIT-GLOBK software package (Herring et al. 2018a, b). In particular, the processing strategy excluded the sites with data from <2.5 years and efficiency <50% (i.e., number of days with regular observations with respect to the whole period of observation). During the processing, a loose constraint approach (i.e., 100 m) was applied to a-priori coordinate of each site, while tight constraints were used for the precise ephemerides and the Earth Orientation Parameters. The FES2004 ocean-loading model [49] and an atmospheric propagation delay based on the global mapping function [50] were used in addition to the absolute antenna phase center model provided by the International GNSS Service for satellites and ground stations. At the end of the first step of the processing, the daily solutions were loosely constrained using the GLOBK package [51] to the European Terrestrial Reference Frame (ETRF2014; [52]) using a seven-parameter Helmert transformation, where the differences between the positions of the 14 International GNSS Service stations surrounding the study area were minimized [53]. The daily time series of the North, East and vertical components estimated was modelled as:

$$y\_j(t\_i) = D\_k + v\_k t\_i + \sum\_{k=1}^{N} d\_{\text{jk}} H(t\_i - T\_k) + \sum\_{k=1}^{M} \left( A\_{1k} \cos\left(\frac{2\pi t\_i}{P\_k}\right) + A\_{2k} \sin\left(\frac{2\pi t\_i}{P\_k}\right) \right) + \varepsilon\_k(t\_i) \tag{11}$$

where *j* represents the North (*j* = 1), East (*j* = 2) and vertical (*j* = 3) components. Assuming a linear trend of displacement in the time series within the analyzed time span, the initial position (*D*) and velocity (*v*) obtained from an ordinary least squares linear regression analysis were used to model the long-term trend. The N discontinuities due to instrumental changes and seismic events that occur near the sites are modelled with the *djk* terms, *εk*(*ti*) is time-dependent noise, and *A*1*<sup>k</sup>* and *A*2*<sup>k</sup>* are the amplitudes of the *M* (*M* ≤ 5) seasonal signals with period *Pk*. During the first phase of the analysis, only the initial position (*Dk*), velocity (*v*), discontinuities (*dki*), and noise (*εk*(*ti*)) were estimated by a weighted least-squares method. To obtain a residual time series, a model of the motion was computed using the parameters estimated and, successively, this motion was removed from the time series. Residuals were used to compute the standard errors of the linear regression analysis and, successively, analyzed using a non-linear least-squares technique to estimate the spectra, following the Lomb–Scargle approach [54,55]. The spectrum of each component was analyzed to estimate the period, *P*, of the five (statistically meaningful) principal signals in the interval between 1 month and half of the observation time span. Periodicities were therefore used in Equation (11) to estimate discontinuities, the intercept, and velocities, with amplitude *Ak* = *A*2 <sup>1</sup>*<sup>k</sup>* + *<sup>A</sup>*<sup>2</sup> <sup>2</sup>*<sup>k</sup>* and phase *<sup>φ</sup>* <sup>=</sup> *arctg <sup>A</sup>*2*<sup>k</sup> A*1*<sup>k</sup>* for each component.

#### **4. Available Data and Processing Workflow**

The interferometric processing discussed here was performed using free and opensource tools and software. Indeed, the complete workflow is based on the sequential use of the software packages of SNAP for interferograms generation, and StaMPS to process persistent scatterer time series, with the following reduction of the atmospheric effects by the model discussed in Section 3.1 and included in TRAIN.

#### *4.1. Sentinel-1A and -1B Radar Dataset*

The satellite radar data used in this study was delivered under the Copernicus program and acquired by the S1A and S1B platforms placed in near-polar Sun-synchronous orbits, with inclination of 98.18◦ at an altitude of 693 km. Each satellite has a repeat cycle of 12 days, while for joint use of S1A and S1B, this is reduced to 6 days. Each satellite has a C-band radar with a frequency range of 4–8 GHz and a wavelength range of 37.5–75.0 mm. The acquisition is carried out with the TOPS acquisition mode, in which the antenna beam

is steered in terms of both range (to acquire data from three sub-swaths, as in ScanSAR mode) and azimuth (from backward to forward for each sub-swath). As a consequence, the targets are illuminated by the entire azimuth antenna pattern, thus strongly reducing the scalloping effect and leading to constant signal-to-noise ratio and ambiguities along the azimuth direction. The main drawbacks are the reduced spatial resolution along the azimuth direction caused by the fast azimuth beam steering and the need for extremely precise coregistration. Single Look Complex (SLC) images acquired by the Interferometric Wide (IW) acquisition have been downloaded from the Sentinel Scientific Data Hub (scihub.copernicus.eu) of the Copernicus Open Access Hub.

The IW SLC images have a nominal resolution of 20 m × 5 m (azimuth × range, respectively) and a ground swath width of 250 km. Moreover, they are composed of three sub-swaths (IW1, IW2, IW3) and by three images in single polarization and six images for dual polarization. Each of the sub-swaths has a different mean incidence angle (IW1, 32.9◦; IW2, 38.3◦; IW3, 43.1◦) and consists of a series of bursts in the azimuth which partially overlap in both directions to guarantee the continuity of the image spatial coverage.

In the present study, the ascending and descending Sentinel-1A and Sentinel-1B IW SLC products reported in Table 1 were selected for the successive processing. The temporal and perpendicular baselines between the master and the slave images are reported in Figure 2.



**Figure 2.** The temporal and perpendicular baselines for the ascending (**a**) and descending (**b**) S1 datasets.

#### *4.2. Software*

The Sentinel Application Platform is an architecture that is provided by the ESA and incorporates all the necessary tools to ensure visualization and processing of the Sentinel mission data. The Sentinel-1 Toolbox includes the functionalities for the processing of both ESA (e.g., Sentinel-1, ERS-1, ERS-2, Envisat) and third party (e.g., Cosmo SkyMed, Radarsat-2, TerraSAR-X, ALOS PALSAR) SAR missions. In particular, the SNAP TOPSAR capabilities allow InSAR processing for integration with the StaMPS package, to provide full PSI time series analysis [27]. StaMPS is a software package that was developed by Stanford University (Stanford, USA), University of Iceland (Reykjavík, IS), Delft University (Delft, NL), and University of Leeds (Leeds, GB) that implements the PSInSAR method to extract ground displacements from time series of synthetic aperture radar acquisitions. StaMPS works even in the case of non-steady deformation [56]. The TRAIN tool is available in StaMPS for reduction of the atmospheric noise. In addition to the phase-based linear correction used in this work, TRAIN includes several other models: phase-based power-law correction; spectrometer correction (Moderate/Medium Resolution Imaging Spectrometer [MODIS/MERIS]); weather model correction (European Reanalysis-Interim [ERA-I]; Modern-Era Retrospective Analysis for Research and Applications [MERRA/MERRA-2]; Generic Atmospheric Correction Online Service [GACOS]; Weather Research and Forecasting model). See [43,44] for further details.

#### *4.3. The Processing Workflow*

The SNAP–StaMPS workflow used in the present study is summarized in Figure 3 and further commented on in this section, with details of the single steps performed during the processing.

The SNAP InSAR processing strategy used in the present study includes all the necessary steps for the preparation of products needed for persistent scatterer processing in StaMPS, which are, in particular:


The processing parameters for PSI analysis by StaMPS are set using the available mt\_prep\_snap script. Parameters like the amplitude dispersion (*Da*) and the number of overlapping pixels in the azimuth (*na*) and range (*nr*) are left as default (*Da* = 0.4; *nr* = 50; *na* = 200), while the number of patches is selected depending on the computational characteristics of the computer. As an example, for the descending orbit, five patches were produced, both for azimuth and range. The amplitude dispersion parameter was left as default because, due to the prevalence of vegetated land cover in the study area, no improvements were found increasing its value as total number of candidate PS.

**Figure 3.** Combined SNAP–StaMPS workflow adopted in the present study for the S1 data processing.

The stack of interferograms is processed using MATLAB procedures implemented in StaMPS, complemented by linear-phase-based tropospheric corrections using TRAIN [43,44]. In more detail, this includes:


At the end of StaMPS processing, the TRAIN linear phase-based correction is applied. Further details about these steps can be found in the reference StaMPS/MTI Manual [56,57].

At this stage, the user can select a reference area/site and introduce an average velocity for it (usually known from the CGNSS data). The whole persistent scatterer dataset will be referred to this reference velocity. However, by default, StaMPS uses the whole extent as a reference area. The mean of the interferometric unwrapped phases from all the persistent scatterers is computed, and its value is set to zero. In this way, the mean stability behavior across the entire scene is assumed. This step is shown in Figure 3 as the SAR calibration and introduces a constraint to a global reference frame. Whenever data have to be visualised in a GIS environment, a vector shapefiles format can be exported, which contains the coordinates of each persistent scatterer and an attribute table including the time series of slant-oriented displacements for relevant epochs and the average velocity from a least squares regression analysis. We added a MATLAB procedure for vector decomposition and computing of the vertical and East-West oriented displacement. In the decomposition procedure, the variation of the incidence angle all over the study area was considered and the correct value for each resolution cell was used.

#### **5. Results**

The interferometric processing of S1 SLC images acquired from ascending and descending orbits generated a geocoded LOS velocity map over the study area. To filter out scatterers that showed low temporal coherence, we assumed as coherent pixels those characterized by Da < 0.4. As indicated, the displacements are initially computed under the hypothesis of mean stability behavior across the entire scene. Then, to align the ascending and descending LOS velocities to a common reference frame, we used the local deformation trend provided by the processing of a single CGNSS station included in the study area during an overlapping period. In particular, the BOLG CGNSS station was chosen as the reference by a criterion based on the time length of the series, number of days with regular observations with respect to the whole period of observation and quality of the monument used in the antenna installation. BOLG is part of the EUREF Permanent GNSS Network. It guarantees a reliable constraint to the global reference frame (see Figure 1). The velocities processed at the BOLG site using Equation (11) and referred to an observation time span that overlaps with the interferometric dataset are: *VUP* = −1.0 ± 1.2 mm/yr, *VE* = 0.3 ± 0.5 mm/yr and *VN* = 4.6 ± 0.5 mm/yr. The residual horizontal components (North and East) are estimated by removing the Eurasian plate movements modelled with the parameters provided by [52]. The CGNSS velocities projected along the ascending and descending LOS directions are processed with uncertainties according to Equations (7) and (8), respectively. The following values are obtained: *VLOS*,*ASC* = −1.4 ± 0.9 mm/yr,

*VLOS*, *DESC* = −1.1 ± 1.0 mm/yr. Figure 4 illustrates the LOS-oriented velocity maps within the study area.

**Figure 4.** Line-of-site (LOS)-oriented and geocoded velocity maps (mm/yr) obtained from the ascending (**a**) and descending (**b**) orbits within the study area. In the legend, red correspond to a range increase (movement away from the satellite). Areas enclosed by rectangles identify zones characterized by significant deformation phenomena: (1) Bologna metropolitan area; (2) municipality of Minerbio; and (3) municipality of Budrio. The four alphanumeric characters code of the CGNSS sites indicate the locations of the sites within the study area. BOLG represents the reference CGNSS site adopted for the calibration for the LOS-oriented velocities to the European Terrestrial Reference Frame, ETRF14.

The uncertainties introduced by the alignment procedure arise from the linear combination of errors characterizing the average velocities from CGNSS and PSI time series adopted for alignment purposes. The uncertainty included in the average velocity of radar targets located in the vicinity of BOLG station could be quantified by a statistical analysis applied to the group of selected PSI time series used. In particular, a sample of 36 slant velocities in both orbits has been used and averaged to find reliable average values to be aligned with slant projected CGNSS velocities. Figure 5 provides a sample of ascending PSI time series for targets located in the vicinity of BOLG reference site. The uncertainty of average slant velocities could be placed at 0.1 mm/yr.

As shown by time series of Figure 5, the scattering of points displacements is within a couple of mm among the epochs. Successively, uncertainties introduced by the alignment procedure are transferred to single slant velocities within the study area and will be further propagated in the final error propagation analysis.

The LOS-oriented velocity maps in Figure 4 show deformation phenomena to different geographic extents. It can be noted, three zones (1, 2, and 3 in Figure 4) where the ground deformation is higher than the surrounding areas. As pointed out by a number of studies carried out in this area [32,47,48,58], the displacement pattern is mainly related to lowering of the ground soil due to mainly anthropogenic causes. Previous studies in the area that were based on satellite radar interferometry did not perform decomposition analysis from ascending and descending slant displacements. However, a contribution from horizontally oriented displacement cannot be excluded a priori. This assumption fails in the case of phenomena characterized by the horizontal displacement of the same order as the vertical one, as detected by previous studies in the Po plain area [47,48,53,58]. These investigations reveal horizontal and vertical displacement rates from GNSS data of comparable magnitude. Thus, to achieve full knowledge of the superficial velocities field as seen by the SAR acquisition geometry, there is the need for vector decomposition analysis from ascending and descending velocity maps acquired in the slant geometry complemented by the error propagation analysis (see paragraph Section 3.4).

**Figure 5.** A sample of LOS-oriented time series provided by the PSI analysis in the vicinity of BOLG station. Each sub-figure (1–6) represents the time series of a single PS.

Before the decomposition procedure can be performed, the ascending and descending dataset are sampled over a common grid with a cell size comparable with the spatial resolution of the Sentinel-1 data. A grid of 20 m × 20 m in the ground range is selected and Equation (9) is used to compute the velocities along the vertical (up) and East-West directions for cells with at least 1 persistent scatterer included in both orbits. After the decomposition procedure with a 20-m grid size, vertical and horizontal velocities for 111,306 targets were obtained. The use of larger grid size in the decomposition analysis produces a decrease in the cells associated with decomposed velocities, and a smoothing effect on the resulting decomposed velocities. For instance, a sampling grid size of 30 m × 30 m corresponds to a reduction of the final cells by 10%. The results obtained from this decomposition procedure are shown in Figure 6, where the displacement oriented in the vertical (up) and East-West directions are shown.

In particular, Figure 6a shows the significant vertical deformation patterns over a wide area in the western sector of the metropolitan area of Bologna (area 1), where vertical rates of up to 20 mm/yr can be seen over an extent of ~200 km2. The horizontal velocity pattern shows increasing eastern movement in the same area where the greatest subsidence velocities are seen. These data suggest that the subsidence processes also involve the horizontal displacement.

Velocity maps visible within the areas 2 and 3 differ and an insight will be provided with the aim to show the ability by the maps to inform on ground deformations even by using decomposed velocities. Figure 7 shows the vertical and horizontal velocity maps within the area 2, located in the municipality of Minerbio (district of Bologna).

**Figure 6.** Velocity maps (mm/yr) oriented in the vertical (up) (**a**) and East-West (**b**) directions. Rectangles represent areas of interest, as in Figure 4. Negative rates in (a) indicate subsiding phenomena. Negative rates in (**b**) indicate West-oriented horizontal displacement, and positive rates indicate East-oriented horizontal displacement.

**Figure 7.** Vertical and horizontal velocity maps (mm/yr) within area 2. Negative vertical rates in (**a**) indicate subsiding phenomena; in (**b**) negative rates indicate West-oriented horizontal displacement, positive rates indicate East-oriented horizontal displacement. Module of horizontal displacements are better described using row instead of points. The white dot denotes the main site of industrial activity.

> As visible in Figure 7a, the area exhibits subsidence rates up to 10 mm/yr. Moreover, a prevalent East-West component is well depicted in Figure 7b over a delimited area. In this area, an industrial activity (gas storage in depleted deposits) takes place. However, the authors do not have any other evidence for this relationship and detailed study of the causes of such superficial deformation field is far from the aims of this study. Similarly, Figure 8a,b, shows vertical and horizontal velocity maps processed in the municipalities of Budrio (area 3).

**Figure 8.** Vertical and horizontal velocity maps (mm/yr) within area 3. Negative vertical rates in (**a**) indicate subsiding phenomena; in (**b**) negative rates indicate West-oriented horizontal displacement, positive rates indicate East-oriented horizontal displacement. Module of horizontal displacements are better described by the use of row instead of points. White dot denotes the main site of industrial activity.

In the remaining part of the investigated areas the horizontal velocity pattern shown on Figure 6b indicates widespread and eastward movement (i.e., positive values) in the range of a few mm/yr, especially in the northern area, where the city of Ferrara is located. The eastern sector of this area was stricken by the seismic events in May–June 2012. This area below the Po plain sedimentary cover is characterized by a buried arcuate thrust system and growth folds related to the Apennines chain [59], the principal lineaments of which are also shown in Figure 6b. The moderate uplift in this area (Figure 6a) might be connected to the tectonic movements of the buried Apennines chain, as suggested by several studies for the eastern movements shown in Figure 6b [60]. Finally, minor contributions to the eastward displacement might derive from actual displacements along the North direction, which would contribute to a lesser extent in the slant direction and in successive decomposition analysis.

To provide quality assessment for the vertical and horizontal velocity maps after the decomposition procedure, the rates produced by the PSI-based analysis were compared with GNSS measurements relevant to the period of the InSAR analysis. The CGNSS observations and time series were specifically processed here for this purpose. All measurements are complemented by the uncertainty levels (provided as standard deviations) and, based on the procedure discussed in Sections 3.4 and 3.5, the error propagation analysis is performed to provide errors in the final comparison between GNSS and decompose PSI velocities. The results are given on Table 2.

**Table 2.** Results of the comparisons between the velocities at the GNSS sites and the average persistent scatterer interferometry velocities for a minimum set of five persistent scatterers located around the CGNSS site at a distance ≤500 m.


#PS, number of averaged persistent scatterers; SD, Sampling distance around the CGNSS sites to achieve the minimum of five persistent scatterers needed to compute the averages.

In this comparison, the PSI-based velocities were computed as the average values from a minimum of five persistent scatterers located in the vicinity of the GNSS site. A circle area with a 50-m radius was used to sample and average the persistent scatterer velocities. To reach the minimum number of required persistent scatterers, the search radius was increased 50 m until the required number of persistent scatterers was achieved. Table 2 shows the number of averaged persistent scatterers and the resulting sampling distance used around the CGNSS sites. These data reported in Table 2 show differences between the GNSS and PSI-based velocities of <±2 mm/yr within the reference period, with larger values detected for the East component for the GNSS site of FERR, and in the vertical (up) component for FERA. FERR and FERA are in the urbanized area of Ferrara, at a distance of about 40 km from the reference site. This distance might have a role. Table 2 shows levels of uncertainties <2 mm/yr in the majority of comparisons.

#### **6. Discussion**

The processing of the full ascending and descending S1 radar data and the use of contemporary GNSS observations provided the complete workflow based on PSI analysis and combined SNAP–StaMPS open-source software for ground deformation monitoring. The uncertainty levels introduced by the calibration procedure were also estimated, and the comparison with the velocities provided by a limited number of GNSS sites provided an assessment of the reliability for depicting the displacement field over an area affected by subsidence phenomena. In particular, the contemporaneity of the ascending and descending SAR data and the GNSS observations used in the present study offered us the opportunity to perform a careful calibration procedure with respect to a global reference frame, and for decomposition analysis with successive assessment of the ground deformation phenomena. A similar approach can be found in [28,29], where the PSI velocities from dual orbits were aligned and then the SAR velocities were corrected using the GNSS values. In [61], a S1 descending orbit, dataset was also processed, and the results were combined with in-situ measurements, such as groundwater and GNSS measurements, for more complete interpretation of the data. However, these previous studies based on S1 data lack contemporaneity between the datasets, and none of them explored the vertical and horizontal displacements fields with calibrations in an absolute reference frame, followed by assessment of the accuracy of the uncertainties.

Moreover, the processing workflow of the S1A and S1B data using open-source tools and then the vector decomposition procedure offers some points of discussion. As shown in Figure 5, the SNAP–StaMPS workflow provided a final persistent scatterer time series that can depict linear trends and periodical signatures along the slant direction with a points dispersion of few mm, with no further filtering procedures. After the persistent scatterer time series was processed, we selected a CGNSS site with reliable time series to align the persistent scatterer velocities from single orbits in an absolute reference frame (ETRF2014).

From the analysis of the reference site 3D GNSS time series of BOLG site, we found a propagated uncertainty level on the LOS-projected velocities at the millimeter level, which further propagates during the decomposition in the vertical (up) and East-West velocities. Such propagated errors are reported in Table 2. The comparison between CGNSS and PSI-derived velocities allowed to define differences of 2 mm/yr in worst cases with uncertainties at mm/yr. In this computation, the incidence angle must be referred to the area where the reference site is placed. The persistent scatterer time series appeared less reliable at the beginning of the reference period, when only the S1A data were available, and a void of data was found in the first half of the year 2016, due to the lack of S1 radar images in the repository used.

The use of GNSS time series processed for even longer time periods in comparison with PSI-derived series is suitable for accuracy assessments. This is due to long-term fluctuations which could be not completely detected by the period of the satellite passages. For this reason, the use of CGNSS velocities referred to different time spans could therefore

affect the calibration procedure. The use of GNSS velocities available at the relevant time is therefore advised. Nevertheless, a comparison of GNSS and PSI techniques based on average annual velocities could be poorly sensitive to noises and unmodelled effects. To remove part of these noises, a comparison at the level of time series of displacements should be used. In Figure 9 the LOS-oriented ascending time series for a single target (black dots, see the PS time series number 6 in Figure 5) is compared with the CGNSS positions projected along the LOS. CGNSS are reported for the full time series (grey dots) and for a subset of points sampled at the time of satellite passages (red dots). An epoch-byepoch comparison could be assessed by the root mean square error (RMSE). The RMSE for comparison of series represented in Figure 9 is 2.9 mm. This is mainly due to periodic fluctuations in the CGNSS time series that are not clearly visible in the response from radar target.

**Figure 9.** Comparison between BOLG and a single PS time series. Both series are oriented along the ascending LOS and aligned with respect to a common first epoch. The time in the x-axis is referred to the period under investigation and is provided in the form day/month/year.

The linear phase-based model adopted in the present study to estimate the atmospheric phase contribution is shown to be suitable for mapping deformation in the study area. In particular, the atmospheric corrections processed for the interferograms affected the hilly area that was part of the central Apennines. After the decomposition, this area shows null vertical deformation rates, with only delimited areas characterized by very small positive rates (Figure 6, southernmost sector). This behavior is consistent with the theories related to the tectonics of the Apennine chain. Residual horizontal velocities visible on Figure 6 in correspondence with the hilly area are mostly due to slope instabilities. Removal of the atmospheric phase with more complex models included in TRAIN could be the objective of further studies.

The workflow is revealed to be useful in the detection of deformation phenomena at the millimeter level, like those affecting the Po Plain sedimentary basin. The spatial and temporal variability of the ground deformation driven by anthropogenic causes, such as gas and groundwater exploitation, might be significant, and this dual approach based on GNSS and InSAR offers a solution. A study of the spatial distribution of the phenomena is also necessary, especially when the techniques adopted have an important difference in terms of the spatial distribution: the single point for the GNSS measurements and a relatively large number of persistent scatterers on the ground for the SAR technique. For these reasons, we processed the same observation time span (2015–2019) for the GNSS and SAR observations, and the comparison between the results obtained shows good agreement between the vertical and horizontal (East-West) velocity values (Table 2). Moreover, the use of dual-orbit radar data for subsidence studies might be a requirement for the S1 satellites that are designed to work with a wider look angle with respect to past sensors, and that are less sensitive to actual vertical motion.

The workflow based on PSI analysis and combined SNAP–StaMPS open-source software has some drawbacks. The processing was time demanding using computers equipped with Intel Core i7 and 24Gb of RAM memory. The processing of ascending and descending dataset in the investigated areas took a couple of months with the majority of the processing time (roughly the 90%) used by SNAP. However, computing can be fastened using high-performance or cloud computing. The open code of StaMPS represents a positive point, and the users are able to introduce personalized routines, as we did in the calibration and vector decomposition steps.

The ground deformation map produced in the present study updates previous investigations commissioned by the Regional Agency for Environmental Prevention, on behalf of the Emilia-Romagna Region. Such studies were initially based on traditional topographic surveying [33]. From 2005, single-orbit SAR interferometry based on radar images provided by ERS, Envisat and Radarsat satellites was used to depict the subsidence phenomena over the area, thus assuming a negligible horizontal motion. These previous studies provided updates up to the year 2016 [62] and a reliable benchmark for the present investigations, even though no studies based on dual-orbit and vector decomposition are available.

The maps of Figures 6–8 depict the current subsidence phenomena in the metropolitan area of Bologna, in addition to ground deformation phenomena that affect the nearby settled areas. This might represent a recent update on the modern subsidence rates, and a valuable tool for stakeholders involved in finding countermeasures to land subsidence after the continuous reduction in the amount of groundwater pumped by supply stations located within the study area [63].

#### **7. Conclusions**

In this study, we presented a PSI-based workflow to process dual-orbit S1 radar data with open-source tools complemented by the use of GNSS observations as constraints for the global reference frame and final accuracy assessment of the vertical and East-West oriented velocity maps.

The workflow allowed the investigation of ground deformation due to subsidence phenomena over large extents of Emilia Romagna (Italy). The combined use of the SNAP and StaMPS processing tools offers an opportunity for users who are interested in processing freely available S1 radar images with calibration of velocity maps and use of algorithms included in TRAIN for the atmospheric phase removal. We have added a procedure to process contemporary CGNSS site velocities to refer differential LOS-projected velocities provided by the InSAR approach to the modern ETRF2014, and an algorithm for decomposition analysis at preferred spatial resolution, with successive accuracy assessment carried out at 10 CGNSS sites. The validation of the velocity maps through the comparison of the decomposed SAR and GNSS annual rates provided differences at the millimeter level with larger values at 2 mm/year, thus showing substantial agreement between the PSI-based and GNSS measurements. Moreover, the computed difference values are compatible with the uncertainty level provided by the error propagation analysis. Although a similar procedure was already presented in other studies, they lacked the contemporaneity between the SAR and GNSS datasets, and none of them explored the vertical and horizontal displacement fields, with alignment in an absolute reference frame complemented by an accuracy assessment with error propagation analysis. Even though the accuracy assessment provided a satisfactory outcome, more effort can be focused on investigation of the algorithms other than the linear phase-based model for the atmospheric phase correction. However, the study area is characterized by a prevalently flat topography, and these data can be considered as an improvement in the understanding of ground deformation processes that affect the area as a response to underground resource exploitation.

**Author Contributions:** Data curation, F.G. and N.C.; Formal analysis, F.M., F.G. and N.C.; Investigation, F.M., F.G. and N.C.; Project administration, F.M.; Software, F.G.; Supervision, F.M.; Writing—original draft, F.G. and N.C.; Writing—review & editing, F.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** Authors are grateful to Angelo Galeandro (Eng.PhD), who supported the development of the decomposition algorithm. The Sentinel-1A and 1B data used in the present study were accessed through the Sentinel Scientific Data Hub (scihub.copernicus.eu) of the Copernicus Open Access Hub. The authors would like to thank the Department of Physics and Astronomy of the University of Bologna, and the Department of Physical Sciences, Earth and Environment of the University of Siena, which have kindly made the GNSS daily observations available. We also thank all public and private institutions that manage and distribute the GNSS data used in this study, and in particular: EUREF Permanent GNSS Network, Italian Space Agency (ASI), Leica Geosystem for the ItaPos network, FOGER (Fondazione dei Geometri e Geometri Laureati dell'Emilia Romagna), GeoTop srl for the NetGeo network.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Monitoring of Land Subsidence in the Po River Delta (Northern Italy) Using Geodetic Networks**

**Nicola Cenni, Simone Fiaschi and Massimo Fabris**


**Abstract:** The Po River Delta (PRD, Northern Italy) has been historically affected by land subsidence due to natural processes and human activities, with strong impacts on the stability of the natural ecosystems and significant socio-economic consequences. This paper is aimed to highlight the spatial and temporal evolution of the land subsidence in the PRD area analyzing the geodetic observations acquired in the last decade. The analysis performed using a moving window approach on Continuous Global Navigation Satellite System (CGNSS) time-series indicates that the velocities, in the order of 6 mm/year, are not affected by significant changes in the analyzed period. Furthermore, the use of non-permanent sites belonging to a new GNSS network (measured in 2016 and 2018) integrated with InSAR data (from 2014 to 2017) allowed us to improve the spatial coverage of data points in the PRD area. The results suggest that the land subsidence velocities in the easternmost part of the area of interest are characterized by values greater than the ones located in the western sectors. In particular, the sites located on the sandy beach ridge in the western sector of the study area are characterized by values greater than −5 mm/year, while rates of about −10 mm/year or lower have been observed at the eastern sites located in the Po river mouths. The morphological analysis indicates that the land subsidence observed in the PRD area is mainly due to the compaction of the shallow layers characterized by organic-rich clay and fresh-water peat.

**Keywords:** land subsidence; Po River Delta; integrated monitoring; time-series analysis

#### **1. Introduction**

River deltas, which host large population and extensive economic activities, are among the territories most vulnerable to land subsidence [1,2]: the effects of this phenomenon are linked to environmental degradation, morphological changes of coastlines and emerged surfaces, damage to buildings, and interruption of services [3–5].

Land subsidence afflicts many areas of the world, in particular the ones located along transitional environments, such as coastal areas, deltas, wetlands, and lagoons, which are becoming increasingly vulnerable to flooding, storm surges, salinization, and permanent inundation [6–9]. In these areas, subsidence can be usually considered as a consequence of a complex combination of natural and anthropogenic factors: the compaction of Holocene sediments, tectonic movements, sinkholes formation, volcanism, thawing permafrost, and the Glacial Isostatic Adjustment (GIA), are generally considered as the main natural sources of land subsidence [10–12]; aquifer-system compaction associated with groundwater/oil/natural gas depletion and storage, drainage of organic soils, underground mining, hydro-compaction and stress given by new constructions, are the principal drivers of the anthropogenic land subsidence [13–19]. Moreover, the effects of climate change can dramatically increase the subsidence-related problems due to the rising of sea levels: the 2012 Intergovernmental Panel on Climate Change (IPCC, www.ipcc.ch, (accessed on 9 April 2021)) report, in fact, highlight an increasing occurrence of coastal and fluvial flooding,

**Citation:** Cenni, N.; Fiaschi, S.; Fabris, M. Monitoring of Land Subsidence in the Po River Delta (Northern Italy) Using Geodetic Networks. *Remote Sens.* **2021**, *13*, 1488. https://doi.org/10.3390/rs13081488

Academic Editor: Lorenzo Solari

Received: 5 March 2021 Accepted: 9 April 2021 Published: 13 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

extreme weather events and sea-level rise as a consequence of climate change during the XXI century. Such events can contribute to changes in the environmental conditions of river deltas with a major impact on both the ecosystems and the human activities in these areas.

The monitoring of these complex areas can be carried out through ground-based and remote-sensing high precision techniques such as repeated precise geometric leveling [20], Global Navigation Satellite System (GNSS) [21], and Interferometric Synthetic Aperture Radar (InSAR) [3].

Historically, the repeated precise leveling for decades has been the primary source of data for the monitoring of land subsidence in plain and deltaic areas. However, the high costs and very long execution times of this technique limit its application in present monitoring activities, despite the high accuracy achievable, in the order of 2 mm/km [20,22].

The GNSS first and, more recently, the InSAR techniques have provided more advanced ways to monitor the evolution in time of ground movements at reduced costs and operational time [23,24]. The InSAR technique has been widely applied by the scientific community to the monitoring of land subsidence in different areas, for example: [8] the Vietnamese Mekong Delta; [9] the Nile delta; [25] the Mexicali Valley (Mexico); [26] the entire central Mexico Region; [27] Miami Beach and Norfolk (USA). Several authors (e.g., [3,28–30]) investigated rates and extension of land subsidence in both the Po plain and PRD area using InSAR data. Other authors (e.g., [31–38]) used GNSS observations acquired using Non-Permanent Sites (NPS) and continuous stations to estimate the land subsidence rate in river deltas, showing that despite its benefit, these methods are strongly limited to the number of points that it is possible to measure.

The combination and integration of InSAR and GNSS observations is potentially the best approach to be adopted to monitor the spatial distribution and the temporal variability of the land subsidence in deltaic areas [39]. The data and information provided by both methods are necessary to protect river deltas and low-lying communities from the risks related to the increasing spread of land subsidence [2].

In this study, we present the new PO DELta NETwork (PODELNET), a low-cost network of NPS and Continuous GNSS stations located in the PRD area (Figure 1). The network was developed to improve the 3D information of land subsidence in the PRD area, which has been monitored until 2016 with only two CGNSS sites located inside the study area (TGPO—Taglio di Po and PTO1—Porto Tolle) and three in the surrounding zones (CGIA—Chioggia, CODI—Codigoro, and GARI—Porto Garibaldi) (Figure 1). Furthermore, the increased number of available GNSS points allows a more robust validation of both the already available and the future InSAR data and a more robust integration of the results obtained with both space- and ground-based approaches with great benefits in terms of 3D information and spatial coverage. This study illustrates and discusses the GNSS kinematic pattern obtained using the data acquired during two PODELNET measurements (2016 and 2018) and the comparison with the available InSAR Sentinel-1A/B results (2014–2017) already published by [3]. The deformation pattern derived from GNSS and InSAR data was also superimposed to the geomorphological map of the study area to highlight the possible correlation between the underground characteristics of the PRD area and the pattern of the observed vertical movements.

**Figure 1.** Geomorphological map of the Po River Delta (PRD) area (scheme integrated and modified from "Geomorphological map of Po Plain at 1:250,000 Scale (1999)"). Blue squares and red diamonds indicate the position of the permanent sites respectively located in the PRD and surrounding area. The yellow circles show the location of the Non-Permanent Sites (NPS) belonging to the PODELNET network; in the figure are also shown the positions of the seven Po River terminal branches.

#### **2. The Po River Delta**

The Po is the largest river in Italy. It opens in the northern sector of the Adriatic Sea with a large delta of about 400 km2 that extends seawards for about 25 km [28] (Figure 1). In the delta, the main river (Po di Venezia) is divided into seven branches: Po di Volano, Po di Goro, Po di Gnocca or di Donzella, Po di Tolle, Po di Pila, Po di Maistra and Po di Levante (Figure 1).

The formation of the modern delta is the result of natural processes and human interactions, such as the filling of the wetlands area and engineering endeavors [40]. The eastern part of the delta is mainly characterized by reclaimed territories, as shown in Figure 1. At present day, the reclamation in the PRD area is strongly reduced and part of the territory (about 200 km2) is characterized by a large system of shallow water bodies. Most of the reclaimed territories are actually below the mean sea level and are poorly supplied by sediments because all its river branches have major artificial levees [41]. The complex-wide sandy beach ridges elongating from south to north can be considered as the natural western border of the reclaimed territories (Figure 1). The compaction of the Plio-Quaternary alluvial deposits in the PRD area is an important driver of the natural

ground settlement in these territories: in particular, [30] suggests that the land subsidence rates in the delta area are significantly correlated with the age of the highly compressible Holocene deposits that compose the shallowest 30–40 m of the sedimentary sequence. Other authors [42,43] obtained similar results by analyzing the correlation between the age of the Holocene sediments and the compaction rate observed in the Southern Louisiana Mississippi delta. Their results show land subsidence rates of about 4–5 mm/year for deposit of ages and thicknesses similar to those present in the PRD.

An additional contribution to the natural land subsidence can be related to tectonic movements: some authors (e.g. [11,44]) suggest a contribution of about 1 mm/year due to a crustal flexing induced by the south-westward subduction of the Adriatic plate under the Apennine belt. This subduction model seems to be not compatible with the evidence provided by seismic surveys, which do not show an evident and well-developed slab [45], and with a maximum depth of more than 90 km of the earthquakes occurring in the Northern Apennines. An alternative model suggests that the present deformation pattern observed in the Apennine belt and the Po Plain is driven by the northward motion of the Adriatic plate, instead of the rollback of the Adriatic subducted margin [46].

The sedimentary basin of the Po Plain is also characterized by a complex system of non-confined, semi-confined, and confined aquifers [47]. The extensive exploitation of these aquifers can be considered as one of the principal causes of the anthropogenic land subsidence observed in both the plain and the delta areas. Because of the important historical and natural patrimony located in the Po Plain and the large concentration of industrial facilities as well as the intensive farming in this area, it is necessary to systematically monitor the occurrence and development of the land subsidence phenomenon.

The first leveling performed in the PRD was at the end of the XIX century (1877–1903) and repeated at the beginning of the 1950s [48]. These surveys were carried out before the Italian economic growth: analyzing these data, the vertical kinematic pattern measured by [49] shows in the eastern sector of the Po Plain values of land subsidence lower than 10 mm/year, while the western part is characterized by even less intense subsidence (up to 2/3 mm/year) and by upward motion (uplift). In detail, the leveling surveys carried out in that period provided average land subsidence rates of about 5–7 mm/year, the highest in the Po Plain; because the leveling measurements were carried out when the contribution of the industrial and agricultural activities of the XX century were negligible, the reported values can be considered representative of the natural contribution to the land subsidence and/or the compaction process due to the land reclamations.

The land subsidence rates increased during the Italian economic growth, after the end of the Second World War: in that period, the multi-aquifers system of the Po Plain was extensively exploited for anthropic uses (industrial and agricultural). The methane-water withdrawals from onshore and offshore reservoirs have also contributed to the increase of the land subsidence values. [28,30,50–52] used precise leveling surveys performed in the 1950–1957 in the PRD and Po Plain to measure land subsidence rates up to 300 mm/year. Subsequently, leveling data obtained from 1962 to 1974, highlighted a progressive reduction of the rates with values of about 30 mm/year, in agreement with the progressive reduction of the anthropogenic withdrawals [50,53]. During the late 1970s and at the beginning of the 1980s, the construction of new public aqueducts exploiting surface waters significantly reduced the aquifer overdraft and yielded a general progressive head recovery with a significant decrease of the land subsidence rate [54]. Recent studies [3,31–33,52,55] indicate a decrease of the land subsidence rates with values less than 10 mm/year, probably as a consequence of the successful policies applied 40 years ago to reduce the anthropogenic deformations in the PRD.

#### **3. The PODELNET Network**

The PODELNET network was developed in 2016 in agreement with the IGMI (Istituto Geografico Militare Italiano) and the Veneto Region (Unità di Progetto per il Sistema Informativo Territoriale e la Cartografia and Unità Organizzativa Genio Civile di Rovigo): it is made of 46 non-permanent sites monitored with the GNSS technique. This low-cost and low-impact network is the first application of the densification at 7 km of the IGM95 network in the Veneto Region and, for this reason, the points are named according to the IGMI procedures.

The PODELNET is used also as a reference frame for the activities carried out by the public authorities in the area, mainly the Veneto Region and the Po Delta reclamation consortium.

The PODELNET sites were chosen to take into account both the existing points belonging to the IGMI leveling network and other points belonging to local Authorities and Institutions: all points are located on stable and consolidated foundations; 16 new sites were also identified to obtain a uniform distribution of the network and avoid gaps. Among the 46 non-permanent GNSS sites of the PODELNET network (Figure 1), 24 were connected with the nearest leveling benchmark, providing the orthometric elevation of the points referred to the last IGMI geometric leveling measurements, made in 2005. The others 22 NPS are the benchmarks of the IGMI leveling network. Each point of the PODELNET is connected with its closest point by using a baseline of between 5 km and 7 km.

The first measurement was performed in June/July 2016, and the second survey was carried out two years later, but in the same months to reduce the potential influence of the seasonal signals in the measured velocities. The amplitudes of these seasonal signals are usually not negligible and can introduce biases in the velocity estimation [56,57]. The two surveys were also designed to reduce the impact of the measurement conditions on the estimation of 46 NPS positions. In particular, the two campaigns were carried out surveying the same baselines, when possible, with the same sampling rates, minimum time stationing, and instruments. The minimum observation time at each site was of 3 h, with a sampling rate of 15 s.

#### **4. Materials, Methods, and Processing**

#### *4.1. GNSS Data and Analysis*

The GNSS observations acquired in 2016 and 2018 were analyzed using the Gamit software (release 10.7) [58]. The Earth Orientation Parameter (EOP) and precise ephemerides provided by the International GNSS Service (IGS) were included in the processing with tight constraints. The IERS/IGS 2003 models were adopted to reconstruct the temporal evolution of diurnal, semidiurnal, and terdiurnal solid Earth tides; the pole tide corrections were applied according to the same IERS standards [58]. The FES2004 model [59] was applied to model the effects of ocean loading. The atmospheric propagation delay was then implemented using the "global mapping function" [60]. We assigned also loose constraints to the coordinates of each station included in the processed network.

The daily loose solutions were translated into a local reference frame through a tridimensional movement, where the three translation parameters were estimated with the Globk-Glorg package software [61], using the coordinates and velocities of three external CGNSS sites located near the PRD area (Figure 1): CGIA, CODI, and GARI. The daily solutions of each survey were also combined in a unique solution using the Globk-Glorg package software [61] to provide a unique value of position for each PODELNET site.

The data of the two CGNSS sites located inside the PRD (TGPO and PTO1, Figure 1) were also added to the processing of the PODELNET. The velocities of these two sites, obtained using only the observations acquired during the surveys of the network, were compared with those obtained using the entire time-series by means of the IGS sites as reference stations, to evaluate the uncertainties due to the use of a local reference system.

The coordinates and velocities of the five CGNSS sites were measured applying a procedure similar to that described in [21,33,62]. The present kinematic patterns in the Italian peninsula and surrounding regions described in these works were upgraded using the entire time-series available for each site from 1 January 2001 to 31 December 2019, aligned to the ITRF2014 reference frame, and in particular to the ETRF2014 realization, where the horizontal Eurasian plate motion was removed [63]. As suggested by [21,33,62],

the time pattern of the north, east and vertical components of daily position can be modeled with the following equation:

$$y\_j(t\_i) = A\_j + V\_j t\_i + \sum\_{l=1}^{N} D\_{jl} H(t\_i - T\_l) + \sum\_{k=1}^{M} \left[ B\_{1jk} \cos\left(\frac{2\pi t\_i}{P\_{jk}}\right) + B\_{2jk} \sin\left(\frac{2\pi t\_i}{P\_{jk}}\right) \right] + \varepsilon\_j(t\_i) \tag{1}$$

where *j* = 1, 2, 3 for the north, east, and vertical components; *Aj* and *Vj* are the intercept and trend/velocity of the best fitting straight-line model, respectively; the *Dj* terms are the *N* instrumental or seismic discontinuities eventually occurred at the Tl epochs; *H* is the Heaviside step function; *B*1*jk* and *B*2*jk* are the amplitude of the five (*M* = 5) more significant seasonal terms with period *Pjkj*; *εj*(*t*) represents the time-dependent noise.

The periods of the five more significant seasonal signals were searched in the estimated power spectrum analyzing the residual time-series by the Lomb–Scargle method [64,65]. The residual time-series were obtained removing the estimated linear trend model with only discontinuities (1). These parameters were computed in the first step of the procedure by a weighted least-squares approach (see [21,33,62] for more details).

The spectrum was analyzed finding the period P of the five statistically meaningful signals in the interval between 1 month and half of the observation period. This condition may represent a problem in the procedure because some climatic and geological/geophysical processes could introduce seasonal signals with relatively long periods (greater than half of the observation period) in the GNSS daily time series. Figure 2 shows the power spectrum of the three components related to the five CGNSS sites located in the PRD: it can be noted that some sites are characterized by seasonal signals with a period greater than half of the observation period (about 5–6 years, e.g., GARI). Hence, we have extended the search of the five statistically meaningful signals in the interval between 1 month and the entire observation period. The velocities obtained using equation 1 and associated uncertainties using the method suggested by [66] are reported in Table 1, and the values of the surrounding sites (CGIA, CODI, and GARI) are used to translate in the same reference frame the measured PODELNET positions.

The values were obtained under the hypothesis that the velocity is constant in the entire period of analysis; these periods are greater than the one between the two PODELNET surveys. With this hypothesis, the velocities obtained using the entire period are equal to the rate between the two PODELNET campaigns. This assumption cannot always be satisfied when the GNSS observations are used to monitor land subsidence because the anthropogenic contributions can cause quick changes in the lowering of the ground surface and, therefore, in the land subsidence velocity. These changes can have a limited temporal extension, for example, less than the recommended 2.5 years, which is the minimum period suggested to obtain reliable GNSS velocities [67]: for these reasons, the values shown in Table 1 could not represent the real velocities occurred between the two PODELNET surveys; moreover, the use of these data in the referencing procedure could introduce not negligible biases in the estimated land subsidence velocities for the PODELNET points. We attempted to overcome this problem by using the Velocity Moving Window Approach (VMWA), as suggested by other authors [24,68].

The daily time-series were analyzed with a moving window of 755 days (the time span between the two PODELNET campaigns) to obtain the "instantaneous" velocity (IV): however, the seasonal signals shown in Figure 2 can introduce noise and biases in the IVs, because they can represent local phenomena not representative of a common signal in the PRD area; for this reason, we analyzed the daily time-series where the meaningful significant (M) seasonal signals obtained with equation 1 were removed. Figure 3 shows the IV time series obtained by a window of 755 days, rejecting the velocities obtained with observations shorter than 300 days. The obtained IV corresponding to the PODELNET surveys is also shown in Table 2.

**Figure 2.** Normalized power spectrum of the five Continuous Global Navigation Satellite System (CGNSS) sites located in the PRD and surrounding areas. The spectrums were estimated with the Lomb–Scargle approach [64,65].

In Table 2 are shown (Campaigns columns) the velocities obtained with only the CGNSS observations acquired during the PODELNET surveys. The differences between the values obtained with VMWA and the ones obtained using the observations acquired during the PODELNET measurements are less than 2 mm/year.

**Table 1.** Velocities in mm/year of CGNSS stations located in the PRD and surrounding areas, adopted in the processing of the network as references (CGIA, CODI, GARI) and benchmarks (PTO1, TGPO). In the first and second columns are reported the CGNSS station codes and the date (T), in decimal year, of the first observation used in the processing, respectively. V are the velocities in the ETRF2014 reference frame [63] of the North (VN), East (VE), and Vertical (VV) components. The uncertainties associated with the velocities were estimated using the method suggested by [66]. The Weighted Root Means Square values (σ) in mm/year of the analyzed time-series are shown in the last three columns.


**Figure 3.** IV time-series of the CGNSS sites located in the PRD area obtained using the Velocity Moving Window Approach (VMWA) with a window of 755 days (21 June 2016–18 July 2018). The windows containing less than 300 days were rejected. The black bars represent the uncertainties associated with the IV. The horizontal red lines represent the velocities obtained using the entire time-series and reported in Table 1. The purple vertical line represents the data related to the last day of the second PODELNET survey (18 July 2018). The IV associated with this data are the values corresponding to the time span between the two measurements.

Table 3 and Figure 4 shows the velocities of the PODELNET NPS sites obtained by means of the differences between the 2016 and 2018 coordinates: the uncertainties associated with these values cannot be estimated using the methods proposed in the literature (e.g., [66,69]), where the time-correlated noise in the GNSS time-series is taken into account. We adopted a conservative uncertainty of about 4 mm/year, twice the maximum difference reported in Table 2. The velocities of the CGNSS sites located in the study area shown a kinematic pattern characterized by velocities of about 1–2 mm/year in the north direction (Table 1), as also found in other studies (e.g., [21,33,46]). The uncertainties

associated with the PODELNET velocities are greater than the horizontal kinematic values, which suggests that these rates cannot be used for a detailed study of the horizontal velocity field.

**Table 2.** Velocities of the three components (North VN, East VE, and Vertical VV), in mm/year, of the continuous GNSS stations estimated by means of the Velocity Moving Window Approach (VMWA columns) using a window of 755 days (from 21 June 2016 to 18 July 2018). The values of the CGNSS sites estimated using only the observations acquired during the two campaigns are shown in the Campaigns columns, in mm/year; in the last three columns (Differences) are reported the differences between the corresponding values. The uncertainties associated with the VMWA velocities were obtained by a weighted least square approach adopted in the estimation of the velocities.


**Table 3.** Velocities of the PODELNET sites were estimated by comparing the positions obtained analyzing the two surveys (2016–2018). The code of the PODELNET sites, defined by the IGMI, is reported in the first column; in the subsequent three columns are shown the geographical coordinates: Longitude and Latitude in centesimal degrees and the Height (H) in the WGS84 ellipsoid, in meters. The velocity components in mm/year along the North (VN), East (VE), and Vertical (VV) directions are also reported. The velocities along the LOS direction of Sentinel-1A/B satellites are shown in the VG column, in mm/year. The InSAR velocities, in mm/year, in the PODELNET sites are reported in VS column; DV are the differences, in mm/year, between the InSAR (VS) and GNSS (VG) velocities (DV = VS − VG). Np is the number of InSAR points used to obtain the VS rate, and R is the minimum radius, in meters, adopted to find the number of points. The values of PTO1 and TGPO permanent sites were obtained using only the observations acquired during the two PODELNET surveys and processed with the data of these campaigns.



**Table 3.** *Cont.*

#### *4.2. InSAR Data and Analysis*

The data and results used in this study are the ones already published in [3]. In that work, the authors analyzed three different InSAR data sets with the Small Baseline Subset (SBAS) [70] approach over the entire PRD area. In particular, they obtained mean velocities and displacement time-series from the C-band ERS-1/2 (1992–2000), ENVISAT (2004–2010), and Sentinel-1A/B (2014–2017) satellites. The SBAS approach takes into account the socalled distributed scatterers, which are objects on the ground with similar backscattering properties within a SAR image resolution cell (ground pixel). For each of the detected ground pixels, the time-series analysis is performed on the stack of available images by following five main steps: (a) each image is connected multiple times (high redundancy) to generate a well-connected interferogram network; (b) the interferograms are generated and

then co-registered to a selected reference image; (c) the phase of the pixels with coherence higher than a selected threshold is unwrapped using a Delaunay Minimum Cost Flow (MCF) method [71]; (d) the atmospheric phase components (Atmospheric Phase Screen) are estimated and removed by using low-pass and high-pass filters; (e) the velocities and displacements time-series are finally calculated and classified in an ArcGIS environment. For further details on the adopted processing workflow and processing parameters please refer to [3].

**Figure 4.** Vertical velocity field obtained analyzing the observations acquired during the two PODELNET surveys. Diamonds and squares show the position of the three CGNSS stations chosen as reference and benchmark sites, respectively. The background image is the geomorphological map shown in Figure 1 (**a**) and the Po Delta lobe generations map ((**b**), modified from Figure 13 in [72]); are also shown the positions of the seven Po River terminal branches.

#### **5. Results of the GNSS-InSAR Integrated Monitoring**

Among the available InSAR datasets, only the Sentinel-1A/B entirely overlaps with the GNSS observations obtained from the CGNSS stations and the PODELNET network. To compare the InSAR- and GNSS-derived velocities, it was necessary to align the two datasets into the same reference frame. Since the InSAR images do not entirely cover the study area and the CGIA and GARI sites are not covered with InSAR data, we have chosen the GNSS velocities of CODI to translate the InSAR velocities from LOS to the GNSS reference system. The velocities of PTO1 and TGPO stations were used to evaluate the accuracy of the translation process.

The first step of the GNSS-InSAR integrated monitoring procedure is to estimate the GNSS values taking into account only the period overlapping with the InSAR observations: we analyzed the daily time-series without discontinuities and principal seasonal signals with the MVWA method, using a window of 1128 days, which is the time span between the first and the last InSAR observations (17 November 2014–17 May 2017). In Figure 5 are reported the obtained IV time-series.

A limitation of the procedure to align the InSAR velocity field with a local GNSS reference frame is related to the estimation of the InSAR rate at the CGNSS site: a simple method is to find the InSAR point closest to the CODI site and compare the two velocities to measure the translation value. However, these points could not be exactly located at the CODI station and therefore the local movements of the InSAR points are potentially not representative of the CGNSS site, or noisy InSAR data could introduce biases in the alignment procedure. To avoid these potential sources of error, we adopted an averaging procedure, where the velocities of the InSAR points located at a distance less than a defined search radius are used to obtain the InSAR velocity. A critical point is to find the correct size of the search radius: values obtained with a small radius can be strongly influenced by outliers, while a large radius can include other processes not belonging to the CGNSS site.

**Figure 5.** IV time-series of the CGNSS sites located in the PRD area obtained using the VMWA with a window of 1128 days (the time span between the first and last InSAR images processed). The windows containing less than 500 days were rejected. We have analyzed the time series using the model described in equation 1 without discontinuities and significant seasonal signals and estimated by a weighted least square method. The black bars represent the uncertainties associated with the IV obtained with the weighted least square method used to compute the rates. The horizontal red lines represent the velocities obtained with the entire time series and are reported in Table 1. The purple vertical line represents the data related to the last day of the InSAR observations (17 May 2017). The IV associated with this data are the values corresponding to the time span between the InSAR images.

The InSAR observations acquired by the Sentinel-1A/B satellites are characterized by a native ground resolution of about 20 m × 5 m while the SBAS technique used to process the data takes into account the averaged contribution of many independent scatterers within a larger ground resolution cell of about 20 m × 20 m [3]. Thus, we tested the stability of the averaged values at different radiuses, starting from 20 m (Figure 6a): it can be noted that increasing the radius size, the InSAR velocity increases, even if the differences are less than the weighted root mean square values adopted as uncertainties associated with the average values (Figure 6a). As shown in Figure 6b, the number of InSAR points in the search domain increases significantly with the size of the radius, but the average distance between the points and the CGNSS site is less than the search radius (Figure 6c). The relatively high number of InSAR points located in the surroundings of the CGNSS site (Figure 6b) can be explained by the presence of the urban area where the CGNSS station is located.

**Figure 6.** Average InSAR velocity along the LOS direction at the CODI CGNSS site. (**a**) Average velocities with different search radiuses (from 20 m to 300 m, with a step of 10 m); the vertical bars are the calculated Weighted Root Mean Square (WRMS) values; (**b**) Number of InSAR points (Num. Point) used to obtain the average rate; (**c**) The mean distance (Mean Dist.) from the CGNSS site to the points included in the search area.

We found that the 30 m radius was the one that better represents the CODI site with the InSAR velocities considering that the differences between the different radiuses are less than the weighted root mean square values assumed as uncertainties associated with the average values (of around 0.6 mm/year) (Figure 6a).

After the correction of the InSAR kinematic pattern in the local GNSS reference system (at CODI site), we have compared the InSAR and the GNSS velocities of the two benchmark sites (PTO1 and TGPO) using the selected 30 m radius, which allowed us to find at least 6 InSAR points. The obtained InSAR velocities in correspondence of the benchmark sites are shown in Table 4 together with the GNSS velocities along the LOS direction obtained using the formulas suggested by [73]: it can be noted that even if the differences with the GNSS values are less than 1 mm/year, we have adopted a conservative threshold of 2 mm/year as the uncertainty associated with the InSAR velocities translated in the GNSS reference frame.

**Table 4.** Velocities of the three components (North VN, East VE, and Vertical VV), in mm/year, of continuous GNSS stations obtained using the Velocity Moving Window Approach (VMWA columns) with a window of 1128 days. The velocities of the site CODI (Figure 1) were used in the reference translation procedure to align the InSAR velocities in the GNSS reference frame. VL, in mm/year, is the GNSS velocity along the LOS direction. VS, in mm/year, is the InSAR rate obtained averaging all the points in a 30 m radius. The differences between the InSAR rate and the LOS GNSS velocities of the benchmarks sites are also reported, in mm/year, in the VD column.


We obtained the average InSAR velocities at each GNSS PODELNET station by selecting at least 6 InSAR points enclosed in a search circle centered on the NPS sites with a radius from 20 m to 300 m and with a step of 10 m. The NPS velocities along the LOS and the differences between these rates and the InSAR values are shown in Table 3: it can be noted that for 7 sites the InSAR velocities are not reported because the number of InSAR points enclosed in the largest search circle is less than 6, probably as a consequence of the high vegetation coverage of these areas.

The minimum number of 6 InSAR points included in the domain around the NPS GNSS sites is not always reached with small circles. For some of the sites, in particular for the ones located in the poorly urbanized territories of the PRD, the minimum search radius necessary to obtain an adequate number of points is greater than 150 m (Figure 7b).

**Figure 7.** Comparison between the InSAR and GNSS velocities; (**a**) The differences DV (VInSAR − VGNSS), in mm/year, between the interferometric and GNSS velocities; (**b**) The InSAR velocities were calculated averaging at least the rates of 6 points located at distances less than the radius. The figure shows the minimum radius needed to find at least 6 points in the searched area; (**c**) differences DV (VInSAR − VGNSS), in mm/year, between the rates compared to the radius used to obtain the InSAR velocities; (**d**) differences DV (VInSAR − VGNSS), in mm/year, compared to the number of points (Np) used to estimate the interferometric rates.

Furthermore, we have investigated if the differences in velocities found in respect to the NPS points are linked to the selected search radius or to the number of InSAR points used to calculate the InSAR velocities: Figure 7c,d shows that the differences are not related to the radius or the number of points included in the domain.

Figure 8a shows the InSAR kinematic pattern in the GNSS local reference frame.

**Figure 8.** (**a**) Velocity map along the LOS obtained analyzing the Sentinel-1A/B observations acquired from 2014 to 2017. The color of the circles represents the LOS PODELNET GNSS velocities (Table 3) overlapped to the InSAR kinematic pattern using the same scale; (**b**) Differences between InSAR and GNSS velocities (Table 3) at the PODELNET sites (VInSAR − VGNSS); are also shown the positions of the seven Po River terminal branches.

The differences between the InSAR and the NPS GNSS velocities along the LOS are shown in Figure 8b. The sites characterized by high positive differences (VInSAR – VGNSS), correspond to the points where the PODELNET data observed high rates of land subsidence. These sites shown in Figure 4a with yellow and green points are probably affected by noise or bias due to some problems that occurred during the surveys. These relatively high velocities can also be related to very localized land subsidence not captured by the InSAR technique. Only future GNSS measurements can provide the information necessary to understand which of the previous two hypotheses are correct.

In the Discussion section we take into consideration only the PODELNET sites where the differences between the InSAR and the GNSS values are less than 10 mm/year (Figure 8b).

#### **6. Discussion**

The vertical kinematic pattern shown in Figure 4 has highlighted several aspects of the spatial and temporal distribution of the land subsidence in the PRD area. The velocities of the analyzed CGNSS sites (Table 1), obtained using the observations acquired in the last decade, are significantly less than -10 mm/year: the values related to the vertical component of PTO1 and TGPO CGNSS stations, are in agreement with those obtained by [55]. The comparison between these values and the ones reported by other authors (e.g., [21,30,31,49,74]) using observations acquired between the past centuries and the first years of the 2000s, clearly indicates a decrease in the land subsidence rate in the PRD area.

The land subsidence rates of about 6 mm/year measured in the last decade (Table 1) agree with the values linked with the consolidation of the late Holocene sediments [30,41,42,52]. However, anthropogenic activities can also contribute to these rates providing variations with periods shorter than the ones of the available CGNSS observations. For this reason, we have investigated the possibility of having short changes of the CGNSS velocities during the entire observation period using the VMWA method. This approach provided also the velocities necessary to align the PODELNET results and the InSAR LOS rates in the same local reference frame, using only the daily positions of the CGNSS stations acquired in the overlapping period between the NPS GNSS surveys (Table 2) and the InSAR observations (Table 4): Figures 3 and 5 show the results obtained with two different windows of about 2 years and 3 years, respectively. The horizontal components (north and east) of the sites located in the PRD area (PTO1 and TGPO) and the surrounding territories (CGIA, CODI, and GARI) are not characterized by significant variations in respect to the average values obtained taking into account the entire time-series (Table 1): the differences are less than the uncertainties associated to the average values (Tables 1, 2 and 4). The east component of the CGIA is characterized by significant variations, especially in the 2017–2019 period of the VMWA time-series obtained with a window of 755 days (Table 2, Figure 3): this result can be explained with changes in the amplitude of the strong annual seasonal signal that affect the temporal evolution of the movements at this station (Figure 2); the variations of the amplitude of the annual signal can be related to the location where the site was installed, more sensitive to the climatic and meteorological conditions of the area.

The vertical components of the VMWA time series of the CGNSS sites show some significant variations (Figures 3 and 5), especially for the stations located outside the PRD area (CGIA, CODI, and GARI). The IV time series obtained with the shorter window (Figure 3) show more significant changes compared to the ones obtained with the 1128 days window (Figure 5). Climatic and/or meteorological processes and human activities are usually characterized by seasonal and/or non-periodic signals with variable amplitudes and relatively short periods (e.g., few months and/or years): these characteristics can be detected by the MVWA method providing variations in the IV time-series obtained with short windows; the analysis performed using large windows can attenuate or remove the variations related to these phenomena in the 'instantaneous' time-series. The variations observed in the PTO1 and TGPO sites seem to present these characteristics: the time-series obtained with short windows (755 days, Figure 3) show some variations in the vertical velocities, while the series obtained with large windows (1128 days, Figure 5) do not show significant changes. Similar results were obtained analyzing the time series of the CGIA site, whereas the other two stations, CODI and GARI, are characterized by significant variations also in the IV time series obtained with a large window. The signals in the GARI and CODI time-series can be linked to the locations of the sites, which can highlight climatic processes or local human activities: for example, the GARI CGNSS site is located on the quay of the Porto Garibaldi harbor and, for this reason, the changes in the vertical velocities could be connected to sea-level changes or local movements of the site related to the tide levels. The analysis of these particular displacements is one of the goals of our future researches: however, the difficulty to obtain adequate information in these sites is one of the major challenges to better understand these temporal variations.

Results of Table 1 (CGNSS sites) indicate also a small but dominant N-S movement, not detectable with the InSAR data. The preliminary horizontal PODELNET velocities reported in Table 3 could indicate that the N-S component is not always dominant. The E-W component in some sites can be possibly related to technical problems that occurred during the measurements or to local effects; therefore, in this case, the comparison with the InSAR velocities can highlight the sites affected by such problems. On the other hand, land subsidence processes can introduce not negligible local E-W movements that can provide a significant contribution to the measured InSAR LOS velocity: in this case, the information inferred from the observations acquired by a ground-based GNSS network can help to

detect the areas not covered by InSAR, and/or to remove possible biases in the InSAR vertical velocity field.

Figure 4a shows the preliminary vertical kinematic pattern obtained analyzing the two PODELNET surveys: the green and yellow sites, characterized by velocities greater than 15/20 mm/year, have not been considered because they are most likely outliers: particular attention must be taken for the future measurements to verify if these high velocities are linked to a particular land subsidence process or just to a technical problem occurred during the surveys. The comparison between the InSAR and GNSS data along the LOS direction shows differences lower than 10 mm/year (Figure 8b): it can be noted that the NPS sites located on the complex sandy beach ridge of the PRD area (Figures 1 and 4a) are characterized by positive velocities or low/moderate negative rates (up to −2 mm/year), whereas the ones on the land reclamation and sedimentary areas are characterized by more heterogeneous velocity patterns, with land subsidence rates lower than 5 mm/year (Figure 4a). However, the obtained kinematic pattern is in agreement with both the results provided by [55] using COSMO-SkyMed and ALOS-PALSAR InSAR data and those of [30] derived from ERS-1/2 InSAR data.

Figure 4b shows the comparison between the vertical kinematic pattern and the highstand PRD lobe map: here there is a good correlation between the PODELNET land subsidence rates and the "younger" territories (0–200 B.P. and 200–350 B.P.) (Figure 4b, dark yellow and pink areas). Land reclamation activities in these areas were carried out by the Republic of Venice starting from the XVII century, with the main intervention that forced the Po river to flow southward through an artificial delta mouth to prevent the sedimentary infilling of the Venetian Lagoon [72]. These results indicate that the land subsidence is mainly due to the compaction of the depositional elements, organic-rich clay, and fresh-water peat located in the shallowest layers of the stratigraphic succession [72]. A similar relationship between sediment age and the land subsidence velocity was found by [30]: they provided land subsidence rates of about 2 mm/year on the sandy beach ridges, and higher velocities (10–15 mm/year) in the eastern deltaic area, in agreement with the results obtained in this work. These authors assumed that only vertical movements contributed to the LOS displacements measured with InSAR: however, the preliminary horizontal velocity field reported in Table 3 indicates that in some areas this hypothesis may not be verified.

The area of the Po di Goro mouth, located in the southern sector of the PRD on 'younger' territories, seems to be characterized by low values of land subsidence, as observed in three different PODELNET sites (Figure 4a,b). A similar pattern can be observed in the Po di Tolle mouth area, where two PODELNET sites present low values of the land subsidence rates (Figure 4a,b). Similar rates are also visible in other two PODELNET sites, located in the area of Albarella Island (Figure 4a,b) and in agreement with the InSAR data: these NPS are located on a sandy beach ridge (Figures 1 and 4a) that correspond to the territories reclaimed during the reorganization of the drainage system occurred in the Medieval age (Figure 4b): this result confirms the previous conclusion about the NPS located on the sandy beach ridges in the western portion of the study area, i.e., the territories reclaimed before the XV–XVI century.

Figure 8a shows the PODELNET and LOS InSAR kinematic patterns: it can be noted that the points located on the sandy beach ridges (Figure 1) and in the western PRD area are characterized by LOS rates lower than the ones located on the easternmost "younger" [72] reclaimed territories, in agreement with the results provided by [30,55]; this confirms that the main contribution to the land subsidence is the compaction of the superficial layers caused by the recent reclamation processes.

#### **7. Conclusions**

The present land subsidence rates in the PRD area were monitored by integration of CGNSS, NPS GNSS, and InSAR data. The results obtained analyzing 5 CGNSS stations show vertical land motions of about −6 mm/year in the deltaic area and about

−3 mm/year in the surrounding territories, without significant variations in the last decade. The preliminary PODELNET NPS GNSS vertical kinematic pattern, obtained analyzing the data acquired through a network of 46 sites in the 2016 and 2018 surveys, shows a heterogeneous land subsidence pattern in the PRD area: the lower vertical land motion rates (values greater than −5 mm/year) are located on the sandy beach ridge that corresponds to the Bronze/Iron Age stranded beach ridges. On the contrary, the higher land subsidence velocities are measured in the easternmost PRD territories formed during the last three centuries. Some areas located in the south-eastern sector of the PRD are characterized by low subsidence velocities (Po di Goro and Po di Tolle mouths): this could indicate more efficient management of the land subsidence processes, but further studies are needed to confirm such hypothesis. Considering the fact that in general the InSAR technology does not perform well in high vegetated areas or areas with high temporal variability, the 46 PODELNET sites represent an important improvement for the monitoring of the land subsidence in the PRD and can represent a strategic monitoring tool for the management of the PRD area in the next future. A scheduled continuous campaign and/or the transformation of some sites in semi-permanent continuous GNSS sites can provide further useful information to better understand the deformations in these territories. An integrated monitoring system, combining GNSS and InSAR data, allows to overcome the limits of the two techniques: the processed InSAR images, validated using the available CGNSS stations, can be integrated with the points of the PODELNET network, increasing the spatial and the temporal resolution together with the cost-effectiveness of both approaches.

This integrated monitoring system find considerable applications in coastal management, especially for the safeguarding of natural ecosystems and human activities; the monitoring of the defense systems against flooding is crucial in areas largely below the mean sea level. Further developments of this study will be linked to the monitoring of the 350 km-long primary levees network of the eastern PRD using new high-resolution InSAR data and GNSS observations.

Finally, the method presented in this work, which includes the integration between different types of data for the evaluation and monitoring of land subsidence, can be applied worldwide in other coastal areas where significant vertical land motion is expected as a consequence of human activities and increasing rates of sea-level rise due to climate changes.

**Author Contributions:** Data curation, N.C., S.F., and M.F.; formal analysis, N.C.; investigation, N.C., M.F., and S.F.; project administration, N.C. and M.F.; software, N.C.; supervision, M.F.; writing original draft, N.C., S.F., and M.F.; writing—review and editing, N.C., S.F., and M.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Department of Civil, Environmental, and Architectural Engineering of the Padova University (Italy) with the Projects "Integration of Satellite and Ground-Based InSAR with Geomatic Methodologies for Detection and Monitoring of Structures Deformations" and "High-Resolution Geomatic Methodologies for Monitoring Subsidence and Coastal Changes in the Po Delta area".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank the following Institutions: EUREF Permanent GNSS Network, Italian Space Agency (ASI), FOGER (Fondazione dei Geometri e Geometri Laureati dell'Emilia Romagna), LEICA-HxGN SmartNet, RING-INGV, GeoTop srl, which have kindly made available GNSS recordings data; the "Laboratory of Geomatics and Surveying" staff of the of the University of Padova, the "Unità di Progetto per il Sistema Informativo Territoriale e la Cartografia" and "U.O. Genio Civile di Rovigo" (Veneto Region) for the contribution in the PODELNET survey activities. Several figures were generated with GMT software [75]. The authors would like also to thank the reviewers and editors for their constructive and helpful comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Remote Sensing* Editorial Office E-mail: remotesensing@mdpi.com www.mdpi.com/journal/remotesensing

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18