Radar coherence and NDVI ratios as landslide early warning indicators

Abstract. The catastrophic failure of the Mud Creek landslide on California’s Big Sur Coast on 20 May 2017 highlighted once again how difficult it is to detect a landslide’s transition from slow moving to catastrophically unstable. Automatic detection methods that rely on InSAR displacement measurements to detect precursory acceleration are available but can be plagued by imaging geometry complexities and tedious processing algorithms. Here, we present a novel approach for assessing landslide stability by using relative interferometric coherence from Sentinel-1 and Normalized Difference Vegetation Index (NDVI) from 5 Sentinel-2. Our method computes the ratio of mean interferometric coherence or NDVI on the unstable slope relative to that of the surrounding hillslope. We show that the coherence ratio of the Mud Creek landslide dropped by 50% when the slide began to accelerate five months prior to its catastrophic failure in 2017. Coincidentally, the NDVI ratio began a near-linear decline. In contrast, the landslide accelerated during the rainy seasons of 2015 and 2016, but neither of those accelerations resulted in a drop of the radar coherence ratio. This suggests that radar coherence and NDVI ratios may be able to aid in both the early 10 detection of landslides and indicate whether an acceleration critically threatens the stability of a slope.

from interferometric synthetic aperture radar (InSAR) data is not without challenge due to imaging geometry complexities and tedious processing algorithms.
More commonly, landslides are mapped after they occur. This practice is particularly important for organizing rescue efforts 25 in response to large numbers of landslides triggered by earthquakes or tropical storms. To this end, landslides are frequently mapped by hand from high-resolution optical images (i.e., Roback et al. (2018)), a process that is tedious and time consuming.
To speed up the production of such landslide inventories, automated and semi-automated procedures have also developed (Mondini et al., 2011). Because landslides frequently damage the vegetation cover, many of the (semi-)automated methods draw on the Normalized Difference Vegetation Index (NDVI; Tucker (1979); Rosenthal and Blanchard (1985)), which can be 30 calculated from red and near-infrared bands of multispectral optical images.
The obvious drawback of approaches that rely on optical imagery is the need for cloud-free imagery, which may not be available immediately after a landslide triggering event (particularly if it was rainfall induced). In contrast, synthetic aperture radar (SAR) has the capability to see through clouds. Due to the increased availability of SAR data (e.g., freely available Sentinel-1 imagery from the European Space Agency), radar-coherence based techniques have recently been explored for the 35 purpose of mapping landslides and damage to infrastructure (Burrows et al., 2019;Yun et al., 2015). Radar coherence is a measure of the similarity of a target's scattering properties between two radar acquisitions (Zebker and Villasenor, 1992). For InSAR applications, where displacements are calculated based on changes of the radar phase, coherence is the primary indicator of data quality. A reduction in radar coherence indicates that either the surface properties of the target have changed or that the imaging geometry has shifted substantially. Because imaging geometries for overlapping Sentinel-1 scenes are generally 40 highly stable, coherence changes are predominantly temporal in nature, making radar coherence a good measure for detecting changes at the Earth's surface. Coherence based landslide mapping has been achieved by classifying a coherence map based on either an absolute coherence threshold, the difference between a pre-event and a co-event coherence map, or the difference between a co-event and a post-event coherence map (Burrows et al., 2019;Yun et al., 2015).
InSAR has also been applied to detect precursory acceleration of the 20 May 2017 Mud Creek landslide in California. 45 Handwerger et al. (2019) showed that the seasonal accelerations of the Mud Creek landslide could be tracked throughout the period for which radar data is available, and that a larger speed-up occurred in the months prior to the failure. An acceleration of a landslide body is considered the best indicator for a future failure; however, determining how much acceleration is indicative of impending failure remains unknown. As with the Mud Creek landslide, unstable slopes frequently accelerate seasonally, or in response to temperature changes or precipitation, without failing catastrophically (Iverson, 2000;Segui et al., 2020;50 Handwerger et al., 2013).
Here, we introduce two novel measures for assessing landslide stability: Radar coherence and NDVI ratio. We evaluate the value of these measures to provide information about the impending failure of the Mud Creek landslide.
The Santa Lucia Mountains rise abruptly from the Pacific Ocean on California's Big Sur Coast, about 150 miles south of San 55 Francisco. Formed in a transpression zone of the San Gregorio-Hosgri Fault System known as the Big Sur Bend, the crest of this rugged, high-relief mountain range rises to over 1700 m asl and is never more than 18 km from the coast (Johnson et al., 2018). Geologically, Miocene marine sediments, Mesozoic to Precambrian granitic and metamorphic rocks, as well as Cretaceous-Jurassic marine sedimentary and metasedimentary rocks make up the Santa Lucia Mountains (Graham and Dickinson, 1978). The Franciscan Mélange that dominates the geology near Mud Creek consists of mesozoic graywacke 60 sandstones, highly sheared argillite shales, metamorphosed greenstones, and conglomerates, and is well known for its highly variable, but generally low, rock strength (Medley and Zekkos, 2011;California Geologic Survey). The Mélange is overlain by unconsolidated, clay rich regolith.
The Big Sur region lies in a Mediterranean climate, where yearly precipitation averages around 1 myr −1 and typically falls between November and April. The total yearly precipitation depends strongly on the storm and drought cycles controlled by 65 the El Niño Southern Oscillation (ENSO). Following a multi-year drought, the winter of 2016/2017 brought an extraordinary number of intense atmospheric-river driven storms to California, resulting in the state's wettest year on record (Swain et al., 2018).
The Mud Creek slide occurred on 20 May 2017, after two weeks with minimal rain. It destroyed almost half a mile of Highway 1, a vital transportation corridor for both tourism and local residents and the only direct connection between Carmel 70 Highlands on the north end of Big Sur and San Simeon to the south. The landslide potential of this particular stretch of road became apparent only in the months preceding the slide, as accumulating debris on the road began to require near daily maintenance (Warrick et al., 2019). Upon recognizing the threat of an imminent landslide, the California Department of Transportation (CalTrans) evacuated all personnel and construction material from the site. Following the catastrophic failure, CalTrans commenced a $54 million project to construct a new road over the landslide deposit. The highway was reopened on 18 July 2018, 75 more than a year after the slide.

Radar data
In this study we us SAR data from ESA's Sentinel-1 satellite to perform a traditional InSAR displacement analysis and to compute the coherence radar. SAR images contain measurements of backscatter amplitude and radar phase for each point on 80 the ground. Radar backscatter is the amount of energy reflected back to the sensor from each ground point. It is influenced by the geometry of the target, its surface roughness, and dielectric properties. InSAR uses the change in the phase of radar waves between measurements to quantify changes at the Earth's surface. Two SAR images acquired by the same satellite of a given area at different times can be processed into interferograms: images that represent the phase difference [0,2π] between the two acquisitions at each point. With knowledge of the radar wavelength, these phase changes can be converted to surface displacements. Because InSAR is sensitive to deformation only in the instrument's line of sight (LOS), three independent measurements are required to obtain the true 3-D motion of a target. In absence of independent measurement, LOS deformations can be projected onto the downslope direction by assuming that the primary motion of a landslide follows gravity.
Alongside any computed interferogram is a coherence (γ) image that serves as the primary quality indicator for InSAR data. It is a measure of how similar the ground properties are at the time of the radar acquisitions (Scott et al., 2017), and 90 is computed from the local phase variance in the interferogram. A loss of coherence can be due to spatial (e.g., a change of satellite viewing geometry) or temporal changes (e.g., change in surface properties like vegetation growth). On a vegetated slope like Mud Creek, coherence is predominantly controlled by three factors: A change in the dielectric properties of the soil due to a change in soil moisture (Molan and Lu, 2020), changes to the geometric characteristics of the surface due to growth or removal of vegetation (e.g., Ruescas et al. (2010)), or displacement of the target that is larger than half a radar wavelength For our analyses, we used JPL's InSAR Scientific Computing Environment (ISCE; Rosen et al. (2012)) and processed 51 Sentinel-1 images from the ascending orbit (track 35) and 64 from the descending orbit (track 42). SAR images for our study site were available back to April 2015, with images typically acquired every 12 days. We obtained Single Look Complex (SLC) images for tracks 35 (ascending) and 42 (descending) from the Alaska Satellite Facility (ASF DAAC: https://search.asf.alaska. 100 edu/#/). We allowed for a maximum time difference of 48 days between images and produced 172 interferograms from the ascending orbit and 208 from the descending orbit (see Table 1). Because of a gap in data availability from the ascending orbit between late 2015 and early 2016, we increased the permissible time difference between images in that period to 6 months.
This allowed us to produce a connected network of interferograms for both orbits. We multi-looked the images with two looks in range and one look in azimuth, and removed the topographic phase with a 10 m DEM obtained from the USGS' National
Images with low overall coherence are typically omitted from InSAR analyses. Because our area of interest is small relative to the size of the full interferogram, mean image coherence is a poor indicator for the data quality in our area of interest.
Rather than manually identifying low quality images, we calculated the mean coherence within our area of interest for each interferogram and only retained images with a mean coherence above a defined threshold. We then computed time series of 110 displacement, radar coherence ratio and amplitude ratio from all the retained images.

Displacement
To produce a displacement time series we used only data from the descending track and removed all interferograms with a mean coherence of less than 0.35. This produced a fully connected time series with 193 interferograms. We selected a point just south of the landslide as our stable reference region and computed a time series using the NSBAS method implemented in 115 JPL's Generic InSAR Analysis Toolbox (GIAnT; Agram et al. (2013)). We retrieved surface displacements by projecting the measured line of sight deformation at each point onto the fall line. To do this we defined the unit line of sight vector a to point from the satellite to the target and the unit fall line b to point from target down slope along the steepest gradient. We describe both vectors in a polar coordinate system where θ describes a positive counter-clockwise rotation from x, and φ describes the angle from positive up (Fig. 2). The full slope parallel deformation D t could then be retrieved as where LOS is the measured line of sight deformation and δ is the angle between a and b, which is defined as:

Amplitude and coherence ratios
To construct a time series of coherence and amplitude evolution, we filtered out all interferograms with mean coherence of less 125 than 0.5 in our area of interest (Fig. 3). This higher threshold was necessary to detect any differences between the unstable part of the slope and the surrounding area. We retained 132 interferograms from the ascending track and 141 interferograms from the descending track after applying this filter criterion. ISCE computes coherence using a 5 x 5 pixel triangular weighted window. For signals s 1 and s 2 , coherence is given by 130 where * indicates the complex conjugate (Jung et al., 2016).
We created a time series of radar coherence ratio (C R ) by calculating the ratio between the mean coherence over the slide (the area that ultimately failed) and the mean coherence over the surrounding slope (termed Reference Slope, see Figure 4) as: Equivalently, we computed the amplitude time series as where γ Slide and σ 0 Slide , and γ Ref Slope and σ 0 Ref Slope are the mean values of radar coherence and backscatter over the slide and reference slope, respectively. The reference slope was determined from optical images and a DEM to include the surrounding hillslope with similar aspect and vegetation pattern. Figure 4 shows the evolution of the coherence on the slide and surrounding hillslope for four interferograms acquired in the fall of 2016 and spring 2017.  7 https://doi.org/10.5194/nhess-2020-227 Preprint. Discussion started: 27 July 2020 c Author(s) 2020. CC BY 4.0 License.

Optical data
To investigate whether coherence changes were driven by changes in vegetation cover, we computed a time series of Normalized Difference Vegetation Index (NDVI), a measure of vegetation productivity that can be derived from optical satellite imagery. NDVI was computed from bands 4 (red) and 8 (near infrared) of ESA's Sentinel-2 (10 m resolution) as By design, NDVI values vary between -1 and 1, where denser or more productive vegetation results in more positive values and sparse vegetation or bare ground results in low or negative values. We calculated the NDVI for all available cloud-free images back to January 2016 and then computed our time series in the same manner that we did for the amplitude and coherence ratios:

Amplitude & coherence ratios
The results of coherence, amplitude, and NDVI ratios are shown in Figures 6 and 7 The amplitude time series on both the ascending and descending tracks show some seasonal trends, but there is no clear difference between 2017 and earlier years. It appears that there is little difference between the backscatter received from the slide body and the surrounding slope in the descending track. Conversely, the backscatter received from the slide body on the 175 ascending track is higher than that of the surrounding slope.

NDVI
The average NDVI ratio of 0.8 prior to 2017 indicates that the Mud Creek slope is generally less densely vegetated than the surrounding hillside. However, starting in early 2017, the gradual decline of the NDVI ratio indicates that vegetation cover within the slide area was degrading. In addition to the ratio time-series, we computed the mean NDVI for the slide area and 180 the reference slope separately to ensure that the evolution of the NDVI ratio is in fact controlled by a decrease of vegetation productivity on the slide and not an increase of productivity on the reference slope (Fig 8).

Discussion
Identifying slopes prone to catastrophic failure, or forecasting when a known landslide will fail, remains a largely unsolved challenge. Five months before its final failure, in conjunction with some of the largest rainfall events measured during Califor-190 nia's 2016/2017 winter, the coherence ratio over the Mud Creek landslide dropped by 50%. At the same time, the NDVI ratio began a near linear decline (see Figure 6). These data raise the question of what such observations reveal about the driving processes and/or whether they can be valuable indicators of impending landslides.
Even where the location of a landslide is known, reliable predictions of the anticipated failure time have only been achieved in a few cases (Intrieri et al., 2019;Loew et al., 2016), and typically require ground-based measurements of displacement 195 rates. As in many cases, the deformation of the Mud Creek landslide was only investigated after the failure (Handwerger et al., 2019(Handwerger et al., , 2015. While this practice provides important insights into landslide mechanics, it does not provide advance warning that can be used to mitigate damage. As Handwerger et al. (2019) show, unwrapping errors make measuring the true surface velocities particularly difficult for fast moving and accelerating landslides. Because radar measurements only reveal phase changes in modulo 2π, a displacement of π looks identical to one of 3π. Unwrapping errors can be reduced 200 by subtracting the mean landslide velocity prior to the phase unwrapping. This can reveal the missing phase cycles and help recover more of the true deformation (Handwerger et al., 2015(Handwerger et al., , 2019. However, this approach requires assessing each landslide individually, and does not lend itself to automatic processing at large scales. We did not apply a correction to subtract the mean landslide velocity and are therefore not able to recover the full deformation. In particular, we cannot recover the acceleration in 2016 due to unwrapping errors that mask the true acceleration that year. Nevertheless, our ∼0.6 m cumulative surface 205 displacement measurements are in good agreement with the ∼0.8 m reported by Handwerger et al. (2019), and the temporal displacement patters beyond 2016 are nearly identical. Thus acceleration prior to the failure is therefore detectable, even if a landslide velocity model is not applied. This suggests that InSAR-derived surface displacements can be used to identify slope instabilities; however, the significant impact of the unwrapping errors highlight that such an approach is not without challenges and may not yield reliable results.

210
Our time series of coherence ratios indicates that a more systematic analysis of these data might provide valuable information for landslide detection and warning. The premise is that general environmental and atmospheric changes (e.g., vegetation cycles, meteorologic storms, ionospheric disturbances) affecting InSAR coherence should not vary between a landslide and its surrounding slopes. If they do, then we should be able to link that difference to changes on the ground. The observed loss of coherence is likely to be caused by one or more factors influencing coherence: Changes to the dielectric constant of the 215 ground, the local surface geometry, or an acceleration of the slope, and we currently can not parse out the contributions of the different factors in detail. However, we can draw on amplitude, NDVI ratios and displacement to shed some light on the temporal evolution of the coherence ratio: Changes of soil moisture alter the dielectric constant of soil, which can lead to a loss of coherence. Moister soils are typically associated with higher radar backscatter (Oldak et al., 2003;Paloscia et al., 2013). The amplitude time series indicates that average backscatter from the slide body was slightly higher (ascending orbit) 220 or equal (descending orbit) to that of the surrounding slope. However, there is no discernible trend that sets 2017 apart from previous years. We interpret these data to show that the Mud Creek slide may have been slightly wetter than the surrounding slope, but that this difference is not significant, and does not seem to persist into 2017. The NDVI ratio began its near linear decline concurrent with an observed the reduction in the coherence ratio. However, the NDVI decline is gradual, indicating that the change of vegetation cover -and therefore the surface geometry -was not sudden. This observation is supported 225 by the reported, continuous sloughing off of vegetation and debris during spring 2017 (Warrick et al., 2019). The removal of vegetation decreases root cohesion (e.g., Schmidt et al. (2001)), promoting additional surface erosion, which would also contribute to an altered surface geometry. Finally, high displacement rates observed in early 2017 can also have contributed to the loss of coherence by surpassing the amount that can be unambiguously determined with InSAR. Combined, these data suggest that, initially, the stark coherence loss may have been driven by the acceleration of the slope in early 2017. Following 230 the initial acceleration, the degradation of the vegetation cover throughout the spring of 2017 may be responsible for and/or have contributed to coherence ratio remaining low until the final failure.
The discrepancies between the low coherence area, the area that exhibits measurable surface displacements and the area that ultimately failed indicate that neither displacement nor vegetation degradation can fully explain the low coherence values; since neither of the two changes were observed in these peripheral areas of the low coherence region (Fig. 5). This observation 235 highlights the need to estimate the individual contributions of different factors to the coherence drop, especially their temporal evolution and interplay.
A main disadvantage of relying on surface displacements remains the viewing geometry. As in the the case of the ascending track at Mud Creek, unsuitable viewing geometries can significantly impede displacement measurements. In contrast, we note that the coherence data from the ascending track -the track that is not well suited to measure surface displacements -provide the 240 more continuous time series of coherence ratio during the spring 2017 acceleration. Therefore, a systematic analysis of radar coherence may complement traditional InSAR measurements in valuable ways. The radar coherence ratio remained steady during the acceleration phases in 2015 and 2016 which were followed by renewed stabilization of the slide. In 2017, during the acceleration period that preceded catastrophic failure, the coherence ratio experienced a step decrease five months prior to failure, while the NDVI ratio showed a linear decline; it is unclear whether the momentary NDVI drop in 2016 was due to the 245 slope acceleration that year. It is evident that NDVI ratios are limited to vegetated slopes, and a more in-depth analysis of NDVI ratio and coherence ratio at multiple landslide sites is necessary to assess their full value. Nonetheless, both the coherence and NDVI approach are simple to compute and could be applied at large scales. Such techniques would be especially beneficial in areas where InSAR-derived pre-failure acceleration is difficult to resolve.

250
Radar coherence has long been used to assess areas of damage after natural catastrophes, but the value of radar coherence as a possible predictor for impending landslides has not yet been studied. In this study we showed that time series of radar coherence ratio and NDVI ratio may be able to distinguish a non critical acceleration of a landslide from one leading to failure many months in advance. Comparatively easy to compute, radar coherence ratios have the potential to be generated at large spatial scales to identify unstable slopes. We believe that the encouraging initial results presented here motivate further 255 investigations of these new parameters.
Code availability. https://github.com/mjacqu/MudCreek Author contributions. Mylène Jacquemart designed the study, processed and analyzed the data, and wrote the manuscript. Kristy Tiampo supervised the work.
Competing interests. Funding for this work came from an NASA's Earth and Space Science Fellowship (M. Jacquemart). The authors declare 260 no competing interests.