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

: 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 speciﬁc 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.


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][4][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][7][8][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][18][19][20][21][22][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.

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]).  [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.

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

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/km 2 . 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 (T REF ) was determined by computing the difference between the weighted mean displacement rate of the PS grid-cell where the GNSS antenna is grounded (VLM PS_GNSS ) and the GNSS velocity (VLM GNSS ) projected to LOS direction [47]: The computed InSAR velocities that were related to LOS direction, LOS displacements , were ultimately recomputed to the vertical displacements, VLM InSAR , using the following formula (e.g., [48][49][50][51]): 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.

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: where SL TG is a linear trend defined from tide gauge measurements, whereas SL SA is a linear trend computed from satellite altimetry. Such computations were performed for three different periods following the GNSS and InSAR data availability.  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.

Integrated VLM Computation Procedures
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.

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. 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]. 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.

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.

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