Persistent Scatterer Interferometric Synthetic Aperture Radar

We used PSInSAR to identify ground deformation in the Prince William Sound region using Sentinel-1 satellite images acquired during snow-free months (typically June–October) from 2016 through 2022. Persistent scatterer interferometric synthetic aperture radar quantifies phase differences at locations that exhibit consistent reflection and high coherence (in other words, there is little loss or change in radar reflectance) through a long period of time. These stable natural and artificial reflectors are known as persistent scatterers (PSs; for example, Ferretti and others, 2001; Hooper and others, 2004). Unlike traditional differential interferometric synthetic aperture radar (InSAR), which measures differences in the phase of returning waves between two SAR images, the use of PSs throughout a long temporal series of InSAR images in the PSInSAR technique allows for the extraction of useful phase information on a sparse grid of targets. This method can correct baseline-induced artifacts, digital elevation model (DEM) errors, and atmospheric disturbances through spatiotemporal regression and filtering. This permits the creation of time series of subtle (on the order of 10s of millimeters) ground deformation (Ferretti and others, 2001; Hooper and others, 2004; Hooper, 2008; Lu and Kim, 2021) and is effective in identifying landslides in paraglacial regions (Kim and others, 2022).

Differential interferograms throughout the study area were generated and stacked using methods described in Kim and others (2022). Stacked interferograms were then processed using the Stanford method for persistent scatterers PSInSAR technique (Hooper and others, 2004; Hooper, 2008). Persistent scatterers candidates were initially chosen based on pixels with lower values of variations in amplitude using the amplitude dispersion index (Ferretti and others, 2001), which corresponds to scatterers with a higher signal-to-noise ratio and thus better quality of phase. The PS probability was then estimated by iterative phase analysis (eq. 1) based upon the phase decorrelation noise, ${\gamma }_x$, which describes the variance in the residuals between observed and estimated phase at the PS candidate pixel (Hooper and others, 2004):

The iterative process continues, and the solution converges, until ${\gamma }_x$ is no longer decreasing (Hooper, 2008; Qu and others, 2022). The resulting phase values of the PS candidates then go through processing steps that consist of DEM error correction, phase unwrapping, atmospheric filtering, and geocoding while maintaining the original resolution of SAR acquisitions (for Sentinel-1, 5 m by 20 m in range and azimuth, respectively) based on the geometry between the satellite orbit and the elevation of each pixel. This results in a time series of phase differences (in other words, ground deformation) for each PS.

Ascending satellite scenes for snow-free months for the years 2016–21 and descending satellite scenes for snow-free months for the years 2016–22 for Sentinel-1 tracks P014 (ascending, 162 scenes between June 9, 2016, and October 24, 2022), P065 (descending, 102 scenes between June 1, 2016, and October 27, 2021), P094 (ascending, 128 scenes between June 3, 2016, and October 29, 2021), and P160 (descending, 114 scenes between August 8, 2018, and October 22, 2022) were used in PSInSAR processing and are available from Alaska Satellite Facility’s vertex platform (https://search.asf.alaska.edu/) and the European Space Agency’s Copernicus Open Access Hub (https://scihub.copernicus.eu/). All movement described is one dimensional in the line of sight (LOS) of the satellite, where positive values indicate that distances between the satellite and the ground surface have become shorter, and negative values indicate that distances have become longer. It is assumed that the rate of deformation during snow-covered (and therefore decorrelated) periods follows that during the previous epoch, and each epoch and its following snow-free season is bridged by InSAR observations to minimize the effect of the winter period of no data.

Landslide Delineation and Confirmation

Landslide presence and precise geospatial extent inherently have some degree of uncertainty because of the partially qualitative nature of landslide identification (Guzzetti and others, 1999). To minimize this uncertainty, we visually assessed the PSInSAR results to identify contiguous areas greater than 50,000 square meters (m²) that exhibited deformation during the analysis period. This area threshold corresponds to the typical landslide minimum size found in an optical landslide inventory of northwest Prince William Sound (Higman and others, 2023), and allows us to reduce the inclusion of false positives because of other deformation phenomena in the area of interest (for example, localized deformation because of glacial processes). The identification of deforming slopes was confirmed using two or more SAR tracks. Once these preliminary areas were defined, we use optical satellite imagery and digital elevation data to identify features indicative of slope instability processes, including scarps, tension cracks, and signs of surface activity (for example, rockfall). Optical satellite imagery with 0.5-m resolution was sourced from Maxar’s WorldView (Maxar, 2023). Digital elevation data sources were obtained from the Alaska Division of Geological & Geophysical Surveys’ Elevation Portal (Alaska Division of Geological & Geophysical Surveys, 2023), including the statewide 5-m interferometric synthetic aperture radar DEM and local airborne light detection and ranging (lidar) datasets.

For each identified instability, a speculative perimeter polygon was manually delineated to calculate landslide area. The perimeter was drawn to circumscribe the PSInSAR points that showed deformation; exceptions to this approach occurred when clear lineations associated with slope deformation (for example, antiscarps, normal scarps, tension cracks, strike-slip faults) were visible in aerial imagery, or if the landslide has been previously delineated using high resolution imagery (for example, Coe and others, 2021). Where possible, the upslope extent of the landslide was defined based on the highest altitude normal scarps, the lateral extents were defined along shear zones, and the lower extent was assumed to be within some area of the slope. Owing to a lack of high-resolution bathymetric data, only subaerial portions of the landslide were considered. Resulting landslide polygons represent only a plausible extent of the landslide boundary, and therefore, landslide delineation and calculated landslide area are estimates that should be revised if additional data become available.

Landslide Volume

Because of uncertainty in the depth and geometry of the landslide failure surfaces, a range of landslide volumes were approximated using logarithmic spirals and volume-area (V-A) relations. Logarithmic spirals were calculated using a python package digger (Barnhart and others, 2023), using methods described in Barnhart and others (2021). Work at Barry Arm (Dai and others, 2020; Schaefer and others, 2023) suggests that the basal failure surface is likely curved, which suggests that a logarithmic spiral geometry is valid to use for landslides in flysch. In polar coordinates (r, θ), a logarithmic spiral is defined as follows:

The elevations and locations of the landslide headscarp and toe imply an average ground surface slope of Φ. Equation 2 and the location of the headscarp and toe along the logarithmic spiral line can be rearranged to generate the following relation (eq. 3) among Φ, the headscarp angle α, and the toescarp angle β:

Additionally, we calculated an estimated volume (V, in cubic meters [m³]) of the landslide based on surface area (A, in m²) using the V-A scaling power law for bedrock landslides defined in equation 4 as follows:

All volume calculations consider only the scenario in which the entire delineated subaerial extent fails catastrophically; in other words, no partial failure scenarios were considered. Thus, none of the reported volume estimates should be considered the smallest possible failure for a given location. Additionally, it is plausible that landslide failure could mobilize adjacent ground not included in the defined extent, which would increase the total volume of failed material. Landslide volume is an important input into the subsequent tsunamigenic plausibility assessment, but volume determination based on aerial extent and V-A scaling and modeling is imprecise and replete with uncertainty. Using a range of potential landslide volumes reflects this uncertainty.

Tsunamigenic Plausibility Assessment

Because of the remoteness of the study area and lack of substantive infrastructure contained within, the primary concern in this initial assessment of landslide hazard is the potential for catastrophic failure of a given landslide into waterbodies, and the corresponding tsunami generation. Given the preliminary nature of this report, we relied upon simple, empirically derived estimates of landslide mobility to characterize the tsunamigenic potential of identified instabilities using the method described in Scheidegger (1973). In other locations, this method has been accepted as a conservative approach for estimating the potential runout and tsunami generation of identified landslides (Hermanns and others, 2012). Two steps are required for this method. First, the potential reach of the landslide was determined using the scaling relations (eq. 5) derived by Scheidegger (1973), where a mobility estimate (F) can be calculated as follows:

The mobility index F is inversely related to landslide travel distance, so a low value of F indicates a long travel distance, and a high value of F indicates a short travel distance. Mobility index F was calculated for all three volume methods described in the “Landslide Volume” section.

Second, we compared F to the mobility required for the landslide to reach the nearest shoreline. Accurate calculation of F_shore requires two processing steps: (1) calculation of the height of the “top” of the landslide, and direction and length of the fall line from the calculated top to the bottom of the landslide (H_slide, L_slide; fig. 2), and (2) calculation of the height and length of the runout path from the landslide to the nearest water body (H_path, L_path; fig. 2). To determine the maximum height used in the calculation, we first identified the slope aspect of the area-weighted center (or centroid) of the landslide polygon to define bearing (compass direction) of the fall line during landslide failure. The fall line (A–A' in figure 2), was then defined as a straight line from the top of the landslide, through the center of mass, and terminating at the bottom of the landslide below the center of mass on that same bearing. The variables H_slide and L_slide were then defined as the elevation difference and planimetric distance from A–A', respectively (fig. 2). Elevation and slope aspect were calculated using the 5-m interferometric synthetic aperture radar-derived DEM available for Alaska (Alaska Division of Geological & Geophysical Surveys, 2023).

Once the fall line direction, distance, and height were determined for the landslide polygon, we calculated the variables H_path and L_path along the path from the bottom of the landslide to the water body (A'–A” in figure 2). Here, we used an eight flow direction flow routing algorithm (Tarboton and others, 1997) to define the location of flow path from the bottom of the landslide to the nearest water body. Application of the algorithm on the original 5-m DEM resulted in an unnecessarily complex flowpath. Therefore, we iterated through several length-scale and smoothing parameter values and resampled the input DEM to 50-m and smoothed it using a mean filter within a circular neighborhood a length-scale (radius) of 10 pixels. These values preserved the broad-scale topographic influence on flow direction while reducing the unnecessary complexity of the defined flow path from the landslide to the nearest water body.

Inland water bodies included in this analysis were initially derived from the USGS National Hydrography Dataset Plus High-Resolution dataset (USGS, 2023). We include all waterbodies >0.5 km² in our analysis. Coastal shorelines were defined using National Oceanic and Atmospheric Administration (NOAA) Center for Environmental Information Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) shorelines data (NOAA, 2023), which were manually revised to better reflect current shoreline locations in areas of rapid glacial recession or advancement in the study period (2016–22). We then calculated the values of H_path and L_path as the elevation difference (using the original 5-m DEM) and planimetric distance from A'–A” (fig. 2), where A” is located at the point at which the flowpath intersects the closest large waterbody (for example, a lake, fiord, or ocean). F_shore is then calculated as:

In the cases where the bottom of the landslide is located directly above a waterbody (for example, B–B' in fig. 2), it is not necessary to calculate a flowpath and instead define F_shore as H_slide / L_slide. The name of the nearest waterbody shoreline for each landslide is provided in table 1.1. Like mobility index F, F_shore is inversely related to landslide travel distance, where a low value of F_shore indicates a long travel distance and a high value of F_shore indicates a short travel distance. We then assigned a “plausible” tsunami potential rating to those landslides where F_shore is greater than or equal to F and an “implausible” tsunami potential rating for those landslides where F_shore is less than F. If any one of the three volumetric methods resulted in a plausible tsunami potential rating, we adopted a conservative approach and assigned the landslide as being plausibly tsunamigenic (refer to table 1.1). Again, this reflects the uncertainty in the area and volume calculations required for this analysis.

Landslide Inventory and Implications

The hazard associated with a possible tsunami likely varies substantially across landslides. Of the 43 deforming landslides, we identify 11 landslides that meet the defined mobility criteria for plausible tsunami generation. Landslide characteristics (for example, size, material properties, failure mechanism, runout behavior), waterbody type (for example, lake, fjord, open ocean), submarine characteristics (for example, water depth, submarine landforms), and landslide runout surface (for example, land, glacier ice, or direct entry into the water) all may affect the resulting tsunami characteristics (for example, wave height). However, a comparison based purely on size suggests that the inventory developed here contains two landslides with tsunamigenic potential (Yale D, Yale F; table 1, app. 1) that are of comparable size or larger than the Barry Arm A landslide, which has the potential to generate a hazardous wave of regional significance (Dai and others, 2020; Barnhart and others, 2021).

In comparing our landslide inventory to those in similar environments, such as Glacier Bay in Alaska (Kim and others, 2022), we find that the rates of landslide velocity are similar, on the order of 4 cm or less per year. These observed landslide velocities classify them as very slow (16 mm/yr to 1.6 m/yr) to extremely slow (less than 16 mm/yr) according to the International Geotechnical Society’s UNESCO Working Party on World Landslide Inventory (WP/WLI, 1995). Although average velocity is useful in defining destructive potential, it does not account for potential acceleration and the speed during catastrophic failure. Additionally, like in Glacier Bay, many of the landslides are near present-day (as of April 2023) glacier termini. Thus, landslide and landslide-generated tsunami hazards may continue to evolve as glaciers continue to thin and retreat.

Limitations

The data presented in this report are intended to provide a preliminary means of assessing rates of deformation within the study timeframe and to assess the potential hazard of the landslides identified using PSInSAR methods. Specifically, we use these results to identify landslides that show evidence of very slow to extremely slow movement and which may generate a tsunami should they fail catastrophically. The methods used for this analysis are not without limitation. An exhaustive compilation of all instabilities in the region, and confirmation of deforming areas as landslides, is challenging because of cloud and vegetation cover, poor lighting conditions, satellite imaging geometry, and the lack of high-resolution topography and bathymetry data. Most of the Prince William Sound region is covered in snow from November to May, which causes C-band Sentinel-1 InSAR observations to lose coherence. Additionally, much of the region is covered by moving glaciers, which cause decorrelation due to movement and melt processes, reducing the number of potential persistent scatterers. Moderately to rapidly moving landslides cannot be detected by InSAR observations because interferometric phases are constrained within 2π modulo (phase periodic with a 2π radian period). In other words, surface motion exceeding a half radar wavelength (for example, movement greater than about 2.8 centimeters between SAR datasets within a given pixel for the Sentinel-1 satellite) results in a large spatial gradient of phase changes and thereby loss of coherence (Lu and Kim, 2021). To correct unwrapped InSAR phase values, the phase variations between neighboring pixels should be less than π. This can result in a biased estimate of landslide velocity. Thus, the velocities presented herein are only accurate within the limits of the C-band PSInSAR method; for example, deformation measured using InSAR and pixel offset tracking of the Barry Arm A landslide indicates that the landslide was moving at a rate of 1–10 meters per year between 2016 and 2020 (Dai and others, 2020; Schaefer and others, 2023), which is not reflected in PSInSAR results herein. Other methods such as lidar differencing, longer-wavelength InSAR, or pixel offset tracking could provide more spatial measurements throughout a given landslide and thus provide a more accurate assessment of landslide velocities. The best methodological approaches for a given landslide is dependent on several landslide characteristics (for example, velocity, geometry, landslide type), which may be elucidated with further investigations. This additionally highlights the need for multi-sensor, multi-method approach for developing landslide inventories and their velocities in these regions.

The viewing geometry of the Sentinel-1 SAR satellite also causes displacement measurements to be affected by slope gradient, slope aspect, and geometric distortion (Natijne and others, 2022); north- and south-facing slopes have poor sensitivity because near-polar-orbiting SAR satellites are measuring distance from a right angle of satellite moving direction. Vegetation can further obscure C-band InSAR detection of landslides; longer-wavelength InSAR observations are more suitable for revealing landslides beneath vegetation and (or) snow cover (Lu and Kim, 2021; Xu and others, 2023). Thus, we consider this inventory to be more complete on vegetation- and snow-free slopes, and on east- and west-facing slopes.

Regardless of satellite orbit direction (ascending or descending track), coherence in the Prince William Sound region decreases with greater temporal separation between interferometric pairs. The number of high coherence points exceeding 0.3, which is often used as threshold of a reliable interferometric phase (Hanssen, 2001), is quickly reduced even during summer months, although multiples of a year (for example, 1, 2, 3, and so forth years for temporal baselines) tend to keep coherence because of the seasonal similarity in the area (fig. 7). Therefore, the dissipating coherence is a major obstacle in making timely, reliable InSAR measurements and identifying and monitoring landslides in this area.

The empirical nature of the V-A scaling and the mobility analysis results in uncertainty in the estimates of runout potential. Landslide volumes were estimated using V-A scaling and runout relations defined in other geographic, geologic, and climatic settings, which may significantly differ in lithology, structure, and origin. To ensure accurate size and volume estimates, more detailed analysis of the morphology and kinematics of the landslides included in this study is required. Additionally, estimates of landslide mobility required to reach the shoreline follow the shortest distance from the bottom of the landslide to the shoreline, which may not reflect the effect of topography (for example, intermediary ridgelines) on the potential travel path of a failed landslide. All mobility estimates are also based upon observed runout behavior of landslides on dry land. Landslides that run out on glaciers, of which there are several potential examples in this study, may exhibit runout characteristics that differ significantly from those forecasted using the methods included in this report (for example, Bottino and others, 2002). Specific analyses of the runout behavior of large rock avalanches in this environmental setting could better constrain the potential for catastrophic failures to enter the water and generate a tsunami.

Prioritizing Future Analyses

This report only depicts the plausibility of a given landslide entering a waterbody during catastrophic failure and makes no attempt to estimate the potential wave height, travel time, runup, inundation, current velocities, duration, or resonance of a landslide-generated tsunami. Additional modeling (for example, George and others, 2017; Barnhart and others, 2021) could be performed to better constrain the potential effect of a future catastrophic failure on communities, mariners, recreationalists, and the natural and cultural resources in, on, and around Prince William Sound. In addition, hazards presented by non-tsunamigenic landslides can be considerable, such as the potential for a landslide dam and subsequent dam-break flood. Persons traveling through, recreating, or performing commercial or subsistence activities should be aware of potential landslide hazards anytime in or below steep, potentially unstable terrain.

This report does not ascribe any likelihood to the catastrophic failure of these landslides. The determination of a likelihood or a probability estimate for any landslide is a difficult proposition. Slope stability and the potential for failure can be affected by many intrinsic and external factors, including material properties, weathering and alteration, changes in stress and strength conditions, hydrometeorological effects, and strong ground motion during a local or regional earthquake. In addition, catastrophic failure of a landslide is also possible without easily detectible precursory movement. Therefore, it must be assumed that landslide hazards are always present in the identified locations. Future investigations into the controls on slope stability, potential triggering conditions, material properties, geologic structure, and monitoring technologies may improve our understanding of the potential for catastrophic failure.

We conclude this discussion with a qualitative prioritization of the identified landslides for future analyses into “high,” “moderate,” and “low” categories (table 2). We assigned a “high” priority ranking to 11 landslides where a tsunami was determined to be a plausible outcome from rapid catastrophic failure because of the potential far-reaching nature of a tsunami. Within this category, landslides above connected water bodies may (for example, Yale D) may be expected to have farther-reaching effects than those above isolated water bodies (for example, Big Johnstone). However, because of boat traffic in these isolated water bodies and recreational activity along their shorelines, they remain in the “high” priority ranking. “Moderate” priorities were assigned to those landslides that did not meet the defined plausibility criteria but were (1) close to the numerical value defining plausibility (within 0.1), (2) proximal to known recreation areas, infrastructure, and (or) city boundaries, and (or) (3) have a glacier composing most of the runout surface which could increase runout travel distance. In addition to those 11 landslides identified as having tsunamigenic plausibility, we assigned a “moderate” priority to 6 landslides. The rationale for this assignment can be found in table 2. The remaining 26 landslides were assigned a “low” priority ranking for future analyses given the lack of known risk should these landslides fail catastrophically. It is important to note that further investigations may reveal characteristics of a given landslides that would modify these categories.