## 3 Data Coverage and Method
### Uncertainty Estimation  
The cumulative phase \(\sum_{i=1}^{K}\varnothing_{i}\) is the sum of phases, including surface deformation \(\sum_{i=1}^{K}\varnothing_{def,i}\), topography errors \(\sum_{i=1}^{K}\varnothing_{topo,i}\), atmospheric delays \(\sum_{i=1}^{K}\varnothing_{atm,i}\), satellite orbit errors \(\sum_{i=1}^{K}\varnothing_{orb,i}\), and thermal noise \(\sum_{i=1}^{K}\varnothing_{noise,i}\) in the stacking process:  
\[ \begin{array}{l}\sum_{i=1}^{K}\varnothing_{i}=\sum_{i=1}^{K}\varnothing_{ def,i}+\sum_{i=1}^{K}\varnothing_{topo,i}+\sum_{i=1}^{K}\varnothing_{atm,i}+\sum_{i=1}^{K} \varnothing_{orb,i}+\sum_{i=1}^{K}\varnothing_{noise,i},\\ \cong-\frac{4\pi}{\lambda}\cdot\sum_{i=1}^{K}t_{i}-\frac{4\pi\Delta}{\lambda \text{K}\text{S}\text{m}}\sum_{i=1}^{K}\mathbb{B}_{\perp i}+\sum_{i=1}^{K} \varnothing_{atm,i}+\sum_{i=1}^{K}\varnothing_{orb,i}+\sum_{i=1}^{K}\varnothing_ {noise,i},\end{array}  \tag{1}\]  
where \(\varnothing_{i}\) is the unwrapped phase, \(K\) is the number of differential interferograms, \(t_{i}\) and \(B_{\perp i}\) are the time span and perpendicular baseline of one interferogram, \(\lambda\) is the wavelength of SAR images, \(R\) is the distance between satellite and ground and \(\Delta z\) is the residue DEM error. In order to obtain the precise surface deformation, the potential spurious phase, i.e., orbital errors, topographic errors, and atmospheric delays, should be taken into account and minimized. The orbital error can be corrected by a 2-D linear or quadratic model in unwrapped interferograms [15] The atmospheric delays relating to topography can be neglected in this study due to the flat terrain, and the variance of neutral atmospheric delay has been mitigated by our modified stacking method. The topographic error also can be neglected here because the cumulative perpendicular baseline is close to 0 m [26], which can cause an error of 0.6 mm, 1.5 mm, and 3.5 mm with a 5 m DEM error for JERS, ENVISAT and ALOS1  
Figure 3: Atmospheric delay and differential interferograms before and after pair wise logic calculation. (**a**) interferogram between image 20090209 and 20070204; (**b**) interferogram between image 20100212 and 20090209; (**c**) stacking interferogram after summing Figure 3a,b, aiming to reduce the atmosphere effects of image 20090209; (**d**) interferogram between image 20100212 and 2070204 after using direct interferometry; (**e**) the coherence histogram of Figure 3c,d. Black dashed lines in Figure 3a,b show areas with serious atmospheric delays and red rectangles in Figure 3c,d represent zones selected for coherence comparison. [PERSON] and blue solid lines in Figure 3e represent the coherence histogram of Figure 3c,d. (Date format: yyyymmdd). respectively. We used several reference points to calibrate the results of three different platforms, assuming these points are stable with respect to the main movement areas. So we obtained the mean velocity (\(V\)) after considering the source errors listed above. \[V=-\frac{\lambda}{4\pi}\sum_{i=1}^{K}\varpi_{i}/\sum_{i=1}^{K}t_{i}, \tag{2}\]  
We calculated the variances of residual error in each selected interferogram and these results (up to 99, 183 and 140 mm\({}^{2}\) for ENVISAT, JERS and ALOS1) indicated that the atmospheric delay is sensitive in this study (see Figure S3). So we adopted a 1-D exponential covariance function \(Cov_{mn}=\sigma^{2}e^{-d_{mn}/\alpha}\) to estimate the uncertainty of mean velocity [27], where \(\sigma^{2}\) is the variance, \(Cov_{nm}\) is the covariance, \(d_{mn}\) is the distance between pixels, and \(\alpha\) is a variable window size of each acquisition named e-folding wavelength. Figure 4 shows the 1-D covariance function of deformation velocity in three different sensors and their variance for JERS, ENVISAT, and ALOS1 are 3.0, 5.1, and 4.1 mm\({}^{2}\)/year\({}^{2}\) with the e-folding wavelengths as 2.0, 2.5, and 8.8 km respectively.