## Data and Methodology
### PSI Technique and Data Processing  
The approach proposed by [PERSON] et al. (2019) can increase the number of PS in natural areas by identification based on phase characteristics and finding points of low amplitude with a stable phase. This algorithm also uses the spatial correlation of the phases rather than phase history. The residual phase value of pixel \(x\) and the interferogram \(i\) contain the sum of five components:  
\[\phi_{x,i}=\phi_{def_{x,i}}+\phi_{\alpha_{x,i}}+\phi_{orb_{x,i}}+\phi_{e_{x,i}} +n_{x,i} \tag{1}\]  
where:  
\(\phi_{def}\)--phase change due to terrain deformation,  
\(\phi_{a}\)--phase change due to atmospheric contribution,  
\(\phi_{orb}\)--phase due to orbit inaccuracies,  
\(\phi_{e}\)--the residual topographic phase due to error in the DEM,  
\(n\)--phase noise. For the mean phase value in the radius \(L\) for pixel \(x\), it can write:  
\[\overline{\phi}_{x,i}=\overline{\phi}_{def_{x,i}}+\overline{\phi}_{\alpha_{x,i }}+\overline{\phi}_{orb_{x,i}}+\overline{n^{\prime}}_{x,i} \tag{2}\]  
The difference of Equations (1) and (2) give:  
\[\phi_{x,i}-\overline{\phi}_{x,i}=\phi_{e_{x,i}}+n_{x,i}-\overline{n^{\prime}} _{x,i} \tag{3}\]  
where \(\overline{n^{\prime}}=\overline{n}+(\overline{\phi}_{def}-\phi_{def})+( \overline{\phi}_{a}-\phi_{a})+(\overline{\phi}_{orb}-\phi_{orb})\). The residual topographic phase is proportional to the perpendicular component of the baseline \(B\)  
\[\phi_{e_{x,i}}=B\perp_{x,i}\cdot K_{e,i} \tag{4}\]  
where \(K\) is a proportional constant estimated in the sense of least squares. For each pixel, \(x\) is defined as a measure based on temporal coherence:  
\[\gamma_{x}=(1/N)\mid\sum_{i=1}^{N}\exp j(\phi_{x,i}-\overline{\phi}_{x,i}- \phi_{e_{x,i}})\mid \tag{5}\] \begin{table}
\begin{tabular}{c c c c c c c} \hline \hline
**Satellite** & **Track** & **Pass** & **Years** & **Master** & **Number** & **Spatial** & **Time** \\  & & & & **Image** & **of Scenes** & **Baseline (m)** & **Baseline** \\ \hline ERS-2 & 36 & Ascending & 1995–2004 & 22.04.1999 & 55 & \textless{}800 & \textless{}365 \\ Envisat & 186 & Descending & 2002–2008 & 11.01.2004 & 21 & \textless{}700 & \textless{}365 \\ Sentinel-1A & 73 & Descending & 2014–2019 & 09.10.2016 & 52 & \textless{}150 & \textless{}180 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Characteristic of the SAR data used in calculations. At the beginning of the algorithm, candidates for PS are identified on the basis of amplitude dispersion with a high threshold value of 0.4. In the neighborhood of each candidate, the mean phase values (Equation (2)) of the other candidates are calculated and estimate \(K\) and \(\gamma_{x}\). In subsequent iterations, candidates with low \(\gamma_{x}\) are removed and are recalculated the means based on the remaining candidates and \(\gamma_{x}\) for each candidate. The \(\gamma_{x}{}^{*}\) threshold that maximizes the number of PS is determined statistically. More details are contained in the [PERSON] et al. [9]. The ERS-2 and Envisat data were processed with the StaMPS scripts [15], in the DORIS and Matlab environments [47; 48]. Phase unwrapping was performed with the Snaphu ver. 1.4.2 application [49]. Wave phase correction in relation to ground surface was performed with the SRTM-C (3\") model [50]. To process Sentinel-1A SAR data, the GMTSAR software [51] was used to calculate the interferograms, whereas the calculation of displacements using the PSI technique was carried out with the StaMPS. The processed series of radar images provided mean values of displacement velocities for persistent scatterers in the direction towards the satellite (LOS velocity in mm/year). The reference system was adjusted to a set of points having similar displacement values. Without the use of SAR corner reflectors integrated with the GNSS station, it is not possible to determine a reliable absolute reference system.