*2.2. Vertical Land Motion Rates*

SAR images from the Advanced Land Observing Satellite (ALOS) L-Band and Sentinel-1A/B C-Band satellites (search.asf.alaska.edu (accessed on 14 May 2021)) were gathered for the periods July 2007 to February 2011 and February 2017 to February 2021, respectively (see details in Table 1). The ALOS data were acquired along ascending orbit geometry, whereas ascending and descending orbits acquisitions were used for Sentinel-1. Interferograms for selected pairs of images were created and unwrapped using GAMMA, an InSAR processing software [40]. The analysis began with co-registering Single Look Complex (SLC) images to a reference image, which includes a standard matching algorithm using a digital elevation model (DEM), precise orbital parameters, and amplitude images [41]. For the Sentinel-1A/B datasets, the step above was followed by an enhanced spectral diversity (ESD) approach [42,43]. Using this dataset, a set of high-quality interferograms was generated, where only those interferograms with short perpendicular and temporal baselines were processed. A multi-looking operator of 20 and 4 pixels in range and azimuth was applied to obtain a ground resolution cell of approximately 46.0 m × 56.4 m (Table 1). The corresponding figures for ALOS are 6 and 3 pixels in range and azimuth, resulting in a ground resolution of 38.5 m × 23.6 m (Table 1).



To calculate and remove the effect of topographic phase and flat earth correction [44], a 1-arcsecond (~30 m) Shuttle Radar Topography Mission DEM [45] and precise satellite orbital information were used. To identify the elite (i.e., less noisy) pixels, only those pixels with average coherence larger than 0.65 were considered in the analysis, a value that has been used in coastal landscapes with conditions similar to the study area [12,46]. To retrieve the absolute (unwrapped) phase values, a minimum cost flow (MCF) algorithm adapted for sparsely distributed elite pixels was applied (see unwrapping in in Figure 3). Although the precise orbits are used, a few interferograms were still affected by a ramp-like signal, which was removed by fitting a second-order polynomial to their unwrapped phase [47]. Several wavelet-based filters were further applied to correct for effects of spatially uncorrelated topography error and topography correlated atmospheric delay [48] (see atmospheric and orbit correction in Figure 3). Subsequently, a re-weighted least square approach was iteratively applied to invert the corrected measurement of the unwrapped phase at each elite pixel and solve the surface deformation time series. The effect of residual atmospheric errors was further reduced by applying a high pass filter based on continuous wavelet transform to the time series of surface deformation at each elite pixel. Finally, an estimate of the long-term line-of-sight (LOS) deformation rate was obtained as the best-fitting linear slope (hereafter velocity) to the time series of surface deformation at each elite pixel. Together with the LOS velocities, an estimation of the residuals for each elite pixel is provided by plotting the standard deviation of the LOS velocities.

For the Sentinel-1A/B dataset, by combining the ascending and descending LOS velocities in Equation (1), the vertical velocities were retrieved in Cartesian coordinates over those areas where scatterers of both tracks overlap or are in close proximity [49].

$$\text{ddz} \approx \text{(DlosAsc} + \text{DlosDesc)} / 2 \cos \theta \tag{1}$$

where DlosAsc and DlosDesc are the LOS velocities along the ascending and descending orbits, θ is the incidence angle, and dz is the projection of the displacement along the vertical Cartesian axis.

Only the Sentinel-1A/B dataset, therefore, yielded true estimates of VLM. In the case of ALOS data, a qualitative comparison was carried out between its LOS velocities and those obtained after processing the Sentinel-1A/B dataset. It should be noted that in using LOS velocities, the horizontal and vertical velocities are not separated for the observed period. The qualitative comparison between the velocities was made after performing an inverse distance weighted (IDW) interpolation model for transforming the discrete points representing VLM and LOS velocities of individual scatterers into a continuous surface. IDW has been used to interpolate InSAR scattered data points where a relationship or influence over neighboring data is not proven [49]. After running the IDW using a 120-m neighborhood search radius weighted with a cubic exponential power, a 100-m resolution, continuous velocity surface was obtained. A Moran's I test was used to support the choice of distances above for the IDW interpolation. The steps and parameters to obtain the VLM

are summarized in Table 1 and Figure 3. Lastly, aggregations of both the VLM and coastline change transects at 1.5 and 5 km resolutions were performed to examine the effect of scale on our observations. The results are presented in the Supplementary Materials.

**Figure 3.** Flow diagram of interferometric processing of InSAR data (modified from Li et al. [50]).

Uncertainties in VLM rate estimates arise from the backscatter signal being affected by soil structure and moisture [51,52], resulting in apparent ground deformation values that may reach 0.5 radians [53]. Considering this value as the maximum uncertainty level, an annual uncertainty of ±0.2 cm/yr and ±0.7 cm/yr is expected for Sentinel-1A/B and ALOS, respectively. Adding to soil variability, the amount of radar wave penetration in areas with vegetation coverage depends on the structure and density of the canopy and the radar wavelengths; accordingly, wavelengths less than 10 cm (C-band, e.g., Sentinel-1A/B) mostly sense the upper portions of vegetation whilst wavelengths longer than 20 cm (Lband, e.g., ALOS) penetrate deeper into the canopy and can interact with the ground [54], yielding higher coherence values than C-band data [55]. Consequently, ALOS has proven to be more effective in reaching scatters below sparse vegetation [56,57].

Last, the deformation values obtained using SAR data were compared and validated with the velocity values from two permanent GNNS stations located in Barranquilla (11.0197◦/−74.8496◦) and Santa Marta (11.2253◦/−74.1870◦). Despite these stations being situated outside the study area, the satellite frames covered their locations, rendering them helpful to establish whether the trends found using these two approaches were alike (see Figure 2 for the location of GNSS stations).
