Remote Sensing Application for Landslide Detection, Monitoring along Eastern Lake Michigan (Miami Park, MI)

: We assessed the nature and spatial and temporal patterns of deformation over the Miami Park bluffs on the eastern margin of Lake Michigan and investigated the factors controlling its observed deformation. Our approach involved the following steps: (1) extracting bluff deformation rates (velocities along the line of sight of the satellite) using a stack of Sentinel-1A radar imagery in ascending acquisition geometry acquired between 2017 and 2021 and applying the Intermittent Small Baseline Subset (ISBAS) InSAR time series analysis method; (2) generating high-resolution (5 cm) elevation models and orthophotos from temporal unmanned aerial vehicle (UAV) surveys acquired in 2017, 2019, and 2021; and (3) comparing the temporal variations in mass wasting events to other relevant datasets including the ISBAS-based bluff deformation time series, lake level (LL) variations, and local glacial stratigraphy. We identiﬁed areas witnessing high line-of-sight (LOS) deformation rates (up to − 21 mm/year) along the bluff from the ISBAS analysis and seasonal deformation patterns associated with freeze-thaw cycles, suggesting a causal effect. The acceleration of slope failures detected from ﬁeld and UAV acquisitions correlated with high LLs and intensiﬁed onshore wave energy in 2020. The adopted methodology successfully predicts landslides caused by freezes and thaws of the slope face by identifying prolonged slow deformation preceding slope failures, but it does not predict the catastrophic landslides preceded by short-lived LOS deformation related to LL rise.


Introduction
The Laurentian Great Lakes, the largest group of freshwater lakes by area and the second largest by total volume, support the fishing, agriculture, and recreation industries for millions of citizens living near its shorelines. Slope failure along the Great Lake's coastlines arouses attention from government agencies and local homeowners [1][2][3][4] and is exacerbated by recent high lake levels (LLs) that intensify slope failures [5,6]. After more than a decade of below-average LLs, Lake Michigan began a rapid rise from a record low in January 2013 to near a record high in July 2020, reaching 6 cm short of the late 1980s LLs [7]. Michigan LLs have been on the rise since 2013 at a rate of 16 cm/year, amounting

Data and Methodology
We adopted a three-fold method. First, the ISBAS method was applied to extract ongoing slow line-of-sight (LOS) deformation over the study area (0.18 km 2 in Miami Park) using the Copernicus Sentinel-1 data (2017-2021), retrieved from Alaska Satellite Facility (ASF) (accessed on 5 June 2021), processed by ESA. The LOS deformation refers to the displacement measured along the direct path of the radar wave from the source to the

Data and Methodology
We adopted a three-fold method. First, the ISBAS method was applied to extract ongoing slow line-of-sight (LOS) deformation over the study area (0.18 km 2 in Miami Park) Remote Sens. 2022, 14, 3474 5 of 23 using the Copernicus Sentinel-1 data (2017-2021), retrieved from Alaska Satellite Facility (ASF) (accessed on 5 June 2021), processed by ESA. The LOS deformation refers to the displacement measured along the direct path of the radar wave from the source to the target. Second, the bluff failures were mapped using high-resolution elevation models and orthophotos collected by multi-temporal UAV surveys. Finally, a comparison was made between the spatio-temporal variations in InSAR-based bluff deformation rates and bluff failure events and other relevant datasets (e.g., LL variations, local lithology) to constrain the nature of the causal effects.

Lake Level
Monthly peak lake level (PLL) and minimum lake level (MLL) (period: October 2000 to January 2022) were obtained for the nearest monitoring station (Station ID: 9087031; Holland, MI; https://tidesandcurrents.noaa.gov/waterlevels.html?id=9087031 [accessed on 5 January 2022]), located some 35 km north of the study area ( Figure 1a). The LL data are maintained by the Center for Operational Oceanographic Products and Services and the National Oceanic and Atmospheric Administration (NOAA). The monthly data includes the day and time of the measured PLL and MLL and is referenced to the International Great Lakes Datum of 1985 (IGLD 85). The difference between IGLD 85 and NAVD 88 datums in the study area was calculated (0.14 m) and applied to estimate LLs relative to the NAVD 88 datum. The conversion enabled comparisons between LLs and the digital elevation models (DEM)s derived from lidar data [57] (spatial resolution: 1 m) used in this study. We used the National Geodetic Survey's online datum conversion tool https://www.ngs.noaa.gov/cgi-bin/IGLD85/IGLD85.prl (accessed on 2 February 2022) for the conversion.

UAV Data and Processing
We conducted six UAV surveys over the study area in July 2017; July, August, September and October 2019; and July 2021 [58,59]. Three of those surveys are reported in this study ( Table 1). The photographs from each survey were processed applying structure from motion (SfM) techniques [60] using the Agisoft Photoscan software, version 1.4.2. We conducted five primary operations in the SfM workflow: (1) a sparse point cloud was generated, where photos are aligned, camera position and orientation of each photo estimated, and tie points used in the generation of a sparse point cloud; (2) ground control points were used for camera refinement, a step during which estimated point coordinates and camera parameters were adjusted by minimizing the sum of reprojection error and reference coordinate misalignment error, (3) a dense point cloud was extracted from depth maps where combined depth maps generated for each photo was transformed into partial dense point clouds and then merged into a dense point cloud, and (4) orthomosaics were generated by merging the original images projected on the object surface and in the selected projection. All SfM products were generated and exported in the NAD83 UTM 16N projection with a NAVD 88 vertical datum. Details on the type of sensors used, the number of images acquired and ground control points (GCP), and accuracy are provided in Table 1.

SAR Data
The Sentinel-1A satellite was launched on 3 April 2014 by the European Space Agency (ESA) to provide C-band continuity following the retirement of ERS-2 in 2010 and the termination of the Envisat mission in 2012 (ESA, 2016). The twin satellite Sentinel-1B was launched on 26 April 2016. It shares the same orbit as Sentinel-1A but with an orbital difference of 180 • . The short revisit time (twelve days on its own and six days with Sentinel-1B) and high spatial resolution (range direction: 5 m; azimuth direction: 20 m; Table 2) makes the Sentinel-1 mission a particularly effective InSAR mission. Over the study area (western lower Michigan), the Sentinel-1 mission collected data only with the Sentinel-1A satellite in ascending orbit (satellite moving north). We used Level-1 Single Look Complex (SLC) Sentinel-1A images (102 scenes; Figure 1a) acquired (March 2017-April 2021) in Interferometric Wide (IW) swath mode in ascending orbit (Track 121; Table 3 and Figure 1a). We downloaded the SAR data from the Alaska Satellite Facility SAR data service (https://asf.alaska.edu/, accessed on 5 June 2021). Table 2. Sentinel-1A ascending data parameters.

Parameters Value
Azimuth angle 301 • Incident angle 29-46 • Spatial resolution in slant range~5 m Spatial resolution in azimuth~20 m Wavelength 5.6 cm Satellite geometry Right Looking The radar is sensitive to displacement along the satellite LOS direction but not to movement orthogonal to the LOS [56,[58][59][60]. Thus, velocity components in the north or south directions will go undetected. The study area displays significant changes in local morphology and experiences temporal variations in surface cover (e.g., soil, snow, leaves), vegetative cover (e.g., shrubs, trees, grass), and soil moisture (e.g., wet or dry conditions). Those factors could induce variations in ground scattering properties between consecutive image acquisitions and pose challenges to using SAR imagery for monitoring displacement. Moreover, the available SAR acquisitions over the study area are limited to the ascending orbit, which hinders the analysis of displacements over all possible slope orientations and prevents the measurement of displacements along the three components [61].  Figure 2a. Figure 2b is a schematic cross-section along slope direction showing the bluff slope and its ablation and accumulation zones with respect to the satellite LOS incidence angle (θ~39.8 • ). Within the ablation zone, the downslope facilitates downward movement of earth material, a pattern indicative of a motion away from the satellite and represented by negative deformation velocities. In the accumulation zone, the displaced earth materials build up, producing a pattern indicative of a motion towards the satellite and represented by positive velocities. It is noteworthy to mention that the accumulation zone (bluff toe) is washed out during high LL periods. Hence, a motion towards the satellite (indicative of significant E-W movement) in the accumulation zone is not expected with ascending Sentinel-1 observations. ) and look direction (α~79.6°) with respect to the bluff's slope orientation (N80W-S80E) plotted over an ArcGIS pro base map. The solid red arrow represents the geographic north direction, the solid white arrow represents the satellite flight direction, and the solid blue arrow represents the satellite look direction (orthogonal to flight direction). (b) Schematic cross-section along slope direction shows the bluff slope and its ablation and accumulation zones with respect to the satellite LOS. The solid green lines represent the displacement vector T along the slip surface, and the black dashed arrows the displacement vector T along the LOS direction. The incidence angle is represented by θ (~40). Source of base map: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community.

InSAR Methodology
We selected the ISBAS method [62,63] to estimate ground deformation over the study area. Unlike the traditional Small Baseline Subset (SBAS) technique [25], pixels that are coherent for a subset of the investigated time period in addition to those deemed coherent throughout the entire period can be processed in ISBAS [62,64]. Therefore, spatial coverage obtained by ISBAS is usually better than that obtained by the traditional SBAS. The reliability of ISBAS application has been validated in a number of landslide monitoring studies [63][64][65]. Our earlier experiments with the classical SBAS [66,67] and persistent scatterer interferometry (PSI) [68] methods failed to produce results over the research area due to temporal decorrelation. Our field observations and high-resolution satellite and and look direction (α~79.6 • ) with respect to the bluff's slope orientation (N80W-S80E) plotted over an ArcGIS pro base map. The solid red arrow represents the geographic north direction, the solid white arrow represents the satellite flight direction, and the solid blue arrow represents the satellite look direction (orthogonal to flight direction). (b) Schematic cross-section along slope direction shows the bluff slope and its ablation and accumulation zones with respect to the satellite LOS. The solid green lines represent the displacement vector T along the slip surface, and the black dashed arrows the displacement vector T along the LOS direction. The incidence angle is represented by θ (~40). Source of base map: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community.

InSAR Methodology
We selected the ISBAS method [62,63] to estimate ground deformation over the study area. Unlike the traditional Small Baseline Subset (SBAS) technique [25], pixels that are coherent for a subset of the investigated time period in addition to those deemed coherent throughout the entire period can be processed in ISBAS [62,64]. Therefore, spatial coverage obtained by ISBAS is usually better than that obtained by the traditional SBAS. The reliability of ISBAS application has been validated in a number of landslide monitoring studies [63][64][65]. Our earlier experiments with the classical SBAS [66,67] and persistent scatterer interferometry (PSI) [68] methods failed to produce results over the research area due to temporal decorrelation. Our field observations and high-resolution satellite and UAV imagery analysis showed that the study area witnesses seasonal changes in the landscape and vegetative cover, including snowfall in winter and loss of leaves by deciduous plants and along its face as well. The classical SBAS method [66,67] requires the interferometric coherence value of a pixel to exceed a threshold (typically 0.3) for all interferograms considered during the observation period. This often results in coverage gaps in SBAS results where the seasonal changes in scattering properties result in low coherence over sections of the interferograms [62]. For those reasons, the ISBAS method was selected over the classical SBAS method and was applied using ENVI SARScape software (Version 5.6.0, Caslano, Switzerland). The improved ISBAS method considers each interferogram individually and selects and unwraps pixels with high coherence for each interferogram separately. A pixel is designated as a high-coherence pixel when the percentage of interferograms with high coherence at that particular pixel exceeds a user-defined threshold, 65% in all of the interferograms [69]. For periods where coherence drops below the threshold, the velocities are interpolated using measurements bracketing (i.e., before and after) those periods. This generally happens because of rainfall, seasonal snow and vegetation cover, large displacement, or bluff erosion [70]. The selected parameters and processing steps to estimate the LOS velocities are summarized in Figure 3.
The workflow of the adopted ISBAS processing from interferogram generation to geocoding is shown in Figure 3. In the first step, the connection graph is constructed to select Sentinel-1 image pairs. Guided by the temporal and spatial baseline thresholds (60 days, 45% critical baseline), the graph defines the combination of SAR image pairs for interferogram generation by selecting the reference and secondary image pairs.
Sentinel-1A Single Look Complex (SLC) ascending dataset (Track 121; Table 3) was used, and precise orbits were applied to refine the orbital information. The images were subset over the study area (Miami Park), and SLC image stacks were coregistered to a single reference image. We used a multi-looking factor of 4 × 1 to increase the interferograms' signal-to-noise ratio (SNR) (which provided a more reliable coherence estimation) and obtain squared pixels. This yielded a pixel size of 13.27. Differential interferograms were then created for image pairs having small temporal (<60 days) and perpendicular baselines (<120 m).
In the coregistration step, a local nonparametric shift estimate was calculated using a DEM and orbital information. A set of windows (64 × 64 in range and 128 × 128 in azimuth) was defined on the reference image to compute the cross-correlation function, the residual parametric shift was estimated, and proper shifts were applied. The Goldstein filter was used to smooth the interferograms before the final interferogram unwrapping [72]. Coherence thresholds of 0.3 were used for phase unwrapping. The 2D phase unwrapping was performed based on the Minimum Cost Flow (MCF) algorithm with a regular grid covering the data extent [73].
A third-degree polynomial was applied to the selected 109 reference points to estimate and remove the phase constant and ramp during the refinement and reflattening steps ( Figure 3). The reference points were selected outside the actively moving landslide area on highly coherent pixels (i.e., coherence more than 0.8), where unwrapped values showed no residuals and phase jumps. The selected points were mainly located in the permanent settlement area of Miami Park and South Haven ( Figure 1a) and in the immediate surroundings, where no displacement is expected. Finally, the pixel values were interpolated using a nearest neighborhood interpolation to preserve the localized deformation signals.
The phase unwrapping process was reapplied to the refined interferograms. The mean velocity field and possible topographic phase residuals were calculated in the first step of ISBAS inversion with a linear deformation model. Topographic phase residuals were subtracted from wrapped interferograms. In the second inversion step, the deformation time series was estimated by the singular value decomposition (SVD) inversion. During unwrapping and inversion steps, ISBAS methodology was applied using a subset of interferograms to estimate deformation velocity for each pixel with a constraint that a pixel should have reliable coherence in at least 65% of all the interferograms [69]. Finally, the generated displacement time series and mean displacement map were geocoded in a UTM reference system. Remote Sens. 2022, 14, x FOR PEER REVIEW 9 of 24 Figure 3. Workflow of the Intermittent Small Baseline Subset (ISBAS) processing for the bluff deformation analysis following procedures described in [71]. The workflow displays the parameters (dashed boxes), the processing steps (solid circles with black background), and the operations (solid rectangles with grey background). UTM refers to Universal Transverse Mercator.
Sentinel-1A Single Look Complex (SLC) ascending dataset (Track 121; Table 3) was used, and precise orbits were applied to refine the orbital information. The images were subset over the study area (Miami Park), and SLC image stacks were coregistered to a single reference image. We used a multi-looking factor of 4 × 1 to increase the interferograms' signal-to-noise ratio (SNR) (which provided a more reliable coherence estimation) and obtain squared pixels. This yielded a pixel size of 13.27. Differential interferograms were then created for image pairs having small temporal (<60 days) and perpendicular baselines (<120 m).
In the coregistration step, a local nonparametric shift estimate was calculated using a DEM and orbital information. A set of windows (64 × 64 in range and 128 × 128 in azimuth) was defined on the reference image to compute the cross-correlation function, the residual parametric shift was estimated, and proper shifts were applied. The Goldstein filter was used to smooth the interferograms before the final interferogram unwrapping [72]. Coherence thresholds of 0.3 were used for phase unwrapping. The 2D phase unwrapping was performed based on the Minimum Cost Flow (MCF) algorithm with a regular grid covering the data extent [73]. Workflow of the Intermittent Small Baseline Subset (ISBAS) processing for the bluff deformation analysis following procedures described in [71]. The workflow displays the parameters (dashed boxes), the processing steps (solid circles with black background), and the operations (solid rectangles with grey background). UTM refers to Universal Transverse Mercator.

Results
This section discusses our findings in three main areas. Firstly, the applied ISBAS technique and the introduced improvements upon its application to the standard SBAS techniques. Secondly, the distribution and magnitude of the slow deformation over the bluff and the factors causing it from the refined ISBAS application. Thirdly, the timing and distribution of landslides (rapid deformation) along the cliff and the factors controlling their occurrence from UAV surveys acquired in 2017, 2019, and 2021.

Refinement of SBAS Application
Meteorological and perpendicular baseline constraints were considered in selecting Sentinel-1 acquisitions for the ISBAS analysis. Figure 4a is a connection graph that displays all of the available Sentinel-1 images (green points) based on their acquisition dates and orbital separation (perpendicular baseline). The yellow point represents the reference image. Each interferometric pair (interferogram) is represented by two green dots (scenes) joined by a line. Figure 4b is a refined connection graph that displays a subset of the images and interferograms shown in Figure 4a. In selecting the subset dataset, we calculated the mean coherence for each of the interferograms. We also omitted the interferograms with low mean coherence values (<0.3), given the difficulty in retrieving phase information from such pixels [71]. The low coherence was mainly due to temporal decorrelation caused by landcover change related to snow accumulation (four scenes) or change in vegetation (five scenes). The refined ISBAS LOS results (Figure 4d) were retrieved by reprocessing the data stack (91 scenes and 229 interferograms; Table 3) after omitting 11 scenes whose overall coherence exceeded a threshold of 0.3. The mean temporal and perpendicular baselines for the Sentinel-1 ISBAS network were 48 and 47 m, respectively. Comparing the initial ( Figure 4c) and refined (Figure 4d) LOS velocities shows more coherent pixels in Figure 4d.  Figure 5a shows the LOS displacement velocity over the study area. The figure reveals stable conditions except for four areas outlined by polygons 1, 2, 3, and 4. Three of those anomalous areas (1, 3, and 4) are located on the steep slopes of the bluff face, and the fourth (Anomaly 2) lies on the plateau but extends westwards into the upper reach of the hill (Figure 5b-d). The average slopes in areas occupied by anomalies 1, 2, 3, and 4 are  Figure 5a shows the LOS displacement velocity over the study area. The figure reveals stable conditions except for four areas outlined by polygons 1, 2, 3, and 4. Three of those anomalous areas (1, 3, and 4) are located on the steep slopes of the bluff face, and the fourth (Anomaly 2) lies on the plateau but extends westwards into the upper reach of the hill (Figure 5b-d). The average slopes in areas occupied by anomalies 1, 2, 3, and 4 are 54 • , 10 • , 56 • , and 26 • , respectively. Relatively high LOS displacement velocities were observed over the area occupied by anomalies 1 (up to −21 mm/year) and 3 (up to −10 mm/year), the areas with the steepest slopes compared to those over areas with moderate slopes (anomaly 2: up to −8 mm/year) ( Table 4). During our field trips, we noticed that the slope over the area occupied by anomaly 1 and its extension eastwards over the plateau was recently modified to construct luxury residential buildings. Thus, caution should be exercised in attributing the observed subsidence to natural factors alone. We also observed layers of sand and clays along the bluff face over areas occupied by anomalies 2 and 4 ( Figure 1f). 54°, 10°, 56°, and 26°, respectively. Relatively high LOS displacement velocities were observed over the area occupied by anomalies 1 (up to −21 mm/year) and 3 (up to −10 mm/year), the areas with the steepest slopes compared to those over areas with moderate slopes (anomaly 2: up to −8 mm/year) ( Table 4). During our field trips, we noticed that the slope over the area occupied by anomaly 1 and its extension eastwards over the plateau was recently modified to construct luxury residential buildings. Thus, caution should be exercised in attributing the observed subsidence to natural factors alone. We also observed layers of sand and clays along the bluff face over areas occupied by anomalies 2 and 4 ( Figure 1f).   We conducted field visits and UAV baseline acquisitions in the summer of 2017 to identify the factors that could cause the observed slow displacement detected from the ISBAS analysis. In Figure 6, we display field and UAV-related observations that could bear on the factors driving the observed displacement, namely scarps (Figure 6b), groundwater discharge locations on the slope (Figure 6c), and cracks ( Figure 6d). Figure 6b shows curvilinear scarps with a vertical displacement of 0.9 to 1.3 m extending in an east-west direction along the upper reach of the slope. Figure 6c shows groundwater discharge at the interface of a sandy soil layer (thickness 0.9 to 1.2 m) with an underlying clay-rich layer (thickness: 0.72 m). Figure 6d shows a well-developed curvilinear crack (length: 19 m) towards the upper reach of the bluff face and multiple smaller, yet visible, subparallel curvilinear cracks downslope from it (Figure 6d; red arrows). We verified the features mentioned above (scarps, groundwater discharge, and cracks) in the field.

Slow Deformation, Spatial and Temporal Distribution, Controlling Factors
The slow displacement over anomalies 2, 3, and 4, the curvilinear cracks and scarps, and groundwater discharge along the bluff are consistent with findings from earlier field experiments in the study area in which networks of piezometers, clinometers, and pole and cable systems were installed to assess the magnitude of, and the factors controlling, bluff deformation [18,19,56]. They noticed the following: (1) landslides occur along sections of the bluff where sand layers alternate with clay units, a setting where perched aquifers are likely to occur, (2) slow and limited vertical displacements (cm/year) and cracks, some of which develop into scarps with progressive deformation, precede landslides, and (3) deformation occurs mainly during the spring season when thawing of the bluff face occurs. Those observations, together with their findings from dewatering experiments where removal of perched groundwater was found to stabilize bluffs, led Chase and Kehew [19] to hypothesize the following mechanism for the failure of bluffs. In the winter season, bluff faces freeze and block groundwater discharge creating elevated pore pressures in soils. Thawing of the saturated bluff sections occurs in the spring season, allowing discharge of excess groundwater and causing slumping of the bluff slope.
Given the similarities in the geologic setting (alternating sand and clay layers, groundwater discharge, cracks, mini-scarps) and deformation rates and characteristics that we report over anomalies 2, 3, and 4 and those reported by Chase and Kehew [19], we suggest that the InSAR-based slow deformation is mainly related to the phenomenon described by Chase and Kehew [74]. If this were true, one would expect to observe slow subsidence rates and pronounced seasonal cycles in the ISBAS LOS displacement velocities over areas occupied by anomalies 2, 3, and 4, but not over the interleaving areas. The former (slow subsidence rates) is due to soil creep preceding landslide deformation, whereas the latter (seasonal cycle) is related to the freeze and thawing of the bluff face. During the winter season, the frozen bluff face blocks groundwater discharge, pore pressure increases in soils, and so does the volume of groundwater behind the frozen bluff face bringing the bluff face closer to the satellite. The opposite occurs in the spring, where thawing commences the discharge of the once blocked groundwater behind the frozen bluff face, causing the bluff face to retreat away from the satellite. Indeed, we observe those features in the ISBAS   Figure 8 shows mosaics of oblique photographs over the study area extracted from our UAV surveys in 2017, 2019, and 2021, and Figure 9 displays enlargements for the regions outlined by the red boxes in Figure 8 (anomalies 2 and 4). Figures 8 and 9 show the distribution of the coastline, bluff, and plateau in 2017, 2018, and 2021. Comparison of those features allows the detection of landslides along the bluff and constrains their timing. If and when a landslide occurs along the cliff, the vegetative cover is replaced by a mass of rock, debris, or earth material that moves down a slope. Inspection of Figure   In contrast, the interleaving areas along the bluff (stable slope) do not show longterm trends (up to −2 mm/year) and exhibit weak seasonal cycles (4 mm). None of the inland areas on top of the plateau (e.g., plateau1, plateau2, anomaly 3) show seasonal cycles (2, 1, and 2 mm, respectively). Two of the three areas do not show long-term subsidence (plateau 1: up to −1, plateau 2: up to −2, respectively), yet anomaly 3 shows high subsidence rates (up to −10 mm/year). The long-term subsidence of anomaly 3 is less likely to be related to soil creep that precedes landslides, given its location on the top of the plateau. As described early, we observe groundwater discharge along the bluff adjacent to anomaly 3 (Figure 6c), which could indicate sand migration along the groundwater flow direction and cause subsidence in the area.

Rapid Deformation: Spatial and Temporal Distribution
We interpret these findings to indicate that the applied ISBAS methodology can be used to identify areas where long-term slow LOS displacement away from the satellite and strong seasonal deformation cycles occur. In our case, those are regions (areas cov-ered by anomalies 2 and 4) where freezing and thawing of groundwater are causing seasonal land deformation. The following section demonstrates that the identified anomalies over the bluff (anomalies 2 and 4) are more susceptible to landslide development (rapid catastrophic deformation). Figure 8 shows mosaics of oblique photographs over the study area extracted from our UAV surveys in 2017, 2019, and 2021, and Figure 9 displays enlargements for the regions outlined by the red boxes in Figure 8 (anomalies 2 and 4). Figures 8 and 9 show the distribution of the coastline, bluff, and plateau in 2017, 2018, and 2021. Comparison of those features allows the detection of landslides along the bluff and constrains their timing. If and when a landslide occurs along the cliff, the vegetative cover is replaced by a mass of rock, debris, or earth material that moves down a slope. Inspection of Figure 8  In the previous section, we suggested that areas along the bluff witnessing prolonged LOS deformation trend away from the satellite and strong seasonal cycles (anomalies 2 and 4) are susceptible to landslide development. A comparison of the 2017 and 2019 acquisitions in Figure 9 shows that this is the case. The red arrows on both images point to landslide locations (areas L1, L2, and L3) where contiguous wide (1370, 490, and 180 m 2 ) vegetated areas were replaced by bare soil. With one exception of limited size and distribution (L4), we did not detect any landslides outside of areas covered by anomalies 2 and 4. The black arrow on the 2017 and 2019 acquisitions points to this landslide's location and the red arrows to three large landslides within anomalies 2 and 4. Those findings suggest that by identifying the regions witnessing prolonged slow LOS deformation away from the satellite and strong seasonal cycles, we identify bluff sections that have been weakened by the groundwater discharge and most likely to develop landslides. This suggestion is corroborated by the analysis of the 2017 and 2019 images but less so by the 2021 acquisition, where many landslides occurred in areas that did not show the early prolonged slow deformation observed over areas occupied by anomalies 2 and 4. In the previous section, we suggested that areas along the bluff witnessing prolonged LOS deformation trend away from the satellite and strong seasonal cycles (anomalies 2 and 4) are susceptible to landslide development. A comparison of the 2017 and 2019 acquisitions in Figure 9 shows that this is the case. The red arrows on both images point to landslide locations (areas L1, L2, and L3) where contiguous wide (1370, 490, and 180 m 2 ) vegetated areas were replaced by bare soil. With one exception of limited size and distribution (L4), we did not detect any landslides outside of areas covered by anomalies 2 and 4. The black arrow on the 2017 and 2019 acquisitions points to this landslide's location and the red arrows to three large landslides within anomalies 2 and 4. Those findings suggest that by identifying the regions witnessing prolonged slow LOS deformation away from the satellite and strong seasonal cycles, we identify bluff sections that have been weakened by the groundwater discharge and most likely to develop landslides. This suggestion is corroborated by the analysis of the 2017 and 2019 images but less so by the 2021 acquisition, where many landslides occurred in areas that did not show the early prolonged slow deformation observed over areas occupied by anomalies 2 and 4. the red arrows to three large landslides within anomalies 2 and 4. Those findings suggest that by identifying the regions witnessing prolonged slow LOS deformation away from the satellite and strong seasonal cycles, we identify bluff sections that have been weakened by the groundwater discharge and most likely to develop landslides. This suggestion is corroborated by the analysis of the 2017 and 2019 images but less so by the 2021 acquisition, where many landslides occurred in areas that did not show the early prolonged slow deformation observed over areas occupied by anomalies 2 and 4.   Figure 10b shows progressive encroachment of the waterfront towards the bluff, reaching its toe by 2019, a setting conducive to bluff toe erosion and landslide development. The LLs fluctuate throughout the year as well. They are at their lowest in February, when the precipitation is locked as snow accumulation on land, and at their highest levels in July. At this time, much of the snow would have melted, and the runoff had made it to the lake. Figure 10b shows progressive encroachment of the waterfront towards the bluff, reaching its toe by 2019, a setting conducive to bluff toe erosion and landslide development. The LLs fluctuate throughout the year as well. They are at their lowest in February, when the precipitation is locked as snow accumulation on land, and at their highest levels in July. At this time, much of the snow would have melted, and the runoff had made it to the lake.

Discussion
Applying standard radar interferometric methods to extract LOS velocities over the Miami Park bluff has been challenging. As described before, the study area witnesses significant seasonal land-cover variations related to snow accumulation in the winter, loss of leaves in the autumn, and precipitation-related variations in moisture content, all of which cause decorrelation between Sentinel-1 scenes. We applied two main refinements over the standard SBAS applications to minimize landcover change-related decorrelations within the study area. We excluded all pairs of scenes whose mean coherence was below a threshold of 0.3. The majority of those scenes are winter scenes that experienced significant landcover variations. We selected the ISBAS method over the classical SBAS, a technique in which a pixel is considered a high-coherence pixel when the percentage of interferograms with high coherence at that particular pixel exceeds a user-defined threshold. Next, we show how the ISBAS method is beneficial, given the nature and style of deformation in the study area.
Our findings indicate that we have two main deformation styles over the study area's bluffs. The first is prolonged slow LOS deformation (mm/year) away from the satellite process related to the freezing and thawing of groundwater. This process does not occur along the entire length of the bluff but is restricted to specific sections where the local geology and groundwater hydrology play a role. For example, the thawing and freezing that destabilizes the bluff's face could be localized and intensified locally by the high clay content, soil saturation, perched aquifers, preferred pathways for groundwater flow, and thickness of loose sediments. The thawing and freezing of the bluff face is also affected by year-to-year climatic variations (e.g., snow accumulation, snowmelt, precipitation, and temperature), which control pore pressure and groundwater table along the bluff and contribute to its instability [13,14]. We suggest that the observed slow LOS deformation away from the satellite over those bluff sections is related to recurring cycles of prolonged slow soil creep-related deformation followed by landslides. We successfully identified two sections of the bluff (areas covered by anomalies 2 and 4) and one over the plateau and extending toward and over the cliff crest (area covered by anomaly 3) that are witnessing this style of deformation. LOS deformation velocities of up to −8 and −18 mm/year were measured for regions covered by anomalies 2 and 4, respectively (Table 4). Our suggestion is supported by the following observations over those two areas: (1) the presence of small scarps and well-developed curvilinear crack at the upper reach of the bluff face (Figure 6b,d), and (2) the slow LOS deformation away from the satellite (Figure 7) followed by landslide development (Figures 8 and 9; Section 4.3).
The ISBAS can potentially detect the prolonged slow creep-related deformation over a particular pixel yet does not yield an incoherent return once a landslide occurs. In other words, it can identify pixels/areas undergoing recurring cycles of prolonged slow soil creep-related deformation that precedes landslide development. Identifying such areas with standard PSI and SBAS techniques is more challenging.
The second process occurs mainly due to a rise in LLs, landward progression of the shoreline, and erosion of the bluff toe and its failure. The nature and style of deformation resulting from the freezing and thawing of groundwater at the bluff face are quite different from those resulting from the rise in LLs. The landslides resulting from the thaw and freeze of groundwater, but not the rise in LLs, are usually preceded by prolonged (multiple years-long) slow deformation detected on ISBAS LOS velocity products. We could not rule out that slow deformation could have preceded LL rise-related landslides. However, if this was the case, the slow deformation was likely short in duration (days to weeks) because it escaped detection by the ISBAS velocity product. The distribution of deformation associated with the freeze and thaw is selective. It depends on the local geology and bluff groundwater hydrology, whereas the deformation related to the rise in LL depends much less on those settings. Inspection of Figures 8a and 10a shows that this is the case. The LL peaked in 2019 and 2020, and so did the bluff failures; by 2021, the bluff failures associated with the peak LL extended across the entire length of the study area (Figure 8a). Our adopted ISBAS applications enable the extraction of reliable deformation rates along the bluff over areas experiencing significant land-cover change. It predicts the occurrence of landslides by identifying the prolonged slow deformation that often precedes slope failures caused by freeze and thaw of the bluff face. The methodology has its limitations. Landslides preceded by short-lived, slow LOS deformation cannot be predicted, as is the case with the ones caused by a rise in LL. The Sentinel-1 mission did not collect descending dataset over the study area, which prevented us from obtaining an additional SAR data set with a different look angle which could have constrained the velocities in the vertical direction. The launch of the NASA-ISRO SAR (NISAR) mission in 2023 could potentially enable measurement of displacement with both ascending and descending geometries. In addition, the NISAR mission's longer wavelength (L-band) will complement the Sentinel-1 mission (C-band) and provide enhanced penetration of the vegetative and snow cover. Another limitation of the study is the lack of recent ground measurements (inclinometers and global positioning systems). However, installing inclinometers on the studied bluffs remains a challenging task due to the recent high LLs and rates of bluff failures. Future studies could also investigate the interplay between (1) the soilvegetation-atmosphere interactions which can induce significant pore pressure variations at shallow depths and cause slope movements, especially along the clay-rich segments of the bluff [75][76][77], and (2) the infiltration of rainwater and snowmelt and bedrock exfiltration that can provide the local trigger for landslides [78,79].

Conclusions
Over the past few years, many landslides have been reported along Lake Michigan's bluff, with the recently reported rise (2014-2020) in LLs. Miami Park in southwest Michigan was selected as a test site to assess the style of deformation along the bluff, its spatial and temporal variations, and to identify the factors controlling the observed deformation. The applications of standard PSI and SBAS methods in Miami Park were hindered by the extensive land-cover changes throughout the year and the rapid catastrophic deformation associated with landslides. Both factors yield incoherent pixels and unresolved LOS solutions over large sections of the bluff. The use of the ISBAS application and the refinements that were introduced minimized the decorrelation across large sectors of the bluff and enabled the extraction of meaningful LOS deformation velocities over those areas. We extracted LOS bluff deformation rates using a stack of Sentinel-1A radar imagery in ascending acquisition geometry acquired between 2017 and 2021. LOS deformation velocities amounting to up to −18 mm/year were observed over three areas; the subsidence was corroborated by field and UAV observations that revealed the presence of curvilinear cracks and miniature scarps. Strong seasonal cycles (up to 8 mm) were observed over the bluff sections experiencing prolonged LOS deformation away from the satellite. We attribute the slow and prolonged subsidence rates to soil creep preceding landslide deformation and the seasonal cycle deformation to the freezing and thawing of the bluff face and associated water saturation. The proposed freeze and thaw-related deformation is supported by the confinement of the prolonged LOS subsidence, the strong seasonal deformation, groundwater discharge, cracks, and miniature scarps to these areas, and their absence in the stable interleaving areas. We demonstrated that the combined use of InSAR time-series analysis and temporal UAV surveys could identify the locations witnessing prolonged slow LOS subsidence along the bluff, namely bluff sections weakened by groundwater discharge and most likely to develop landslides.
We identified a second style of deformation along the bluff caused by erosion of the bluff toe by rising LLs. This failure is not preceded by prolonged slow deformation as with the freeze and thaw-related landslides and is less dependent on the local hydrology and geology. Inspection of temporal UAV images and our field observations revealed a few landslide occurrences between 2017 and 2019 acquisitions and many more between the 2019 and 2021 acquisitions at PLL. Our methodology can successfully predict the occurrence of landslides caused by freeze and thaw of the bluff face by identifying prolonged LOS subsidence that precedes slope failures. However, it cannot predict landslides without prior deformation or landslides preceded by short-lived deformation resulting from the rise in LLs.  Data Availability Statement: All the Sentinel-1 SAR dataset used in this research are publicly available on https://asf.alaska.edu/ (accessed on 25 May 2022). The LL data can be accessed at https://tidesandcurrents.noaa.gov/waterlevels.html?id=9087031 (accessed on 25 May 2022), and UAV datasets are available in part at https://www.usgs.gov/products/data/data-releases (accessed on 25 May 2022) [58,59] and through the senior author.