Spaceborne InSAR mapping of landslides and subsidence in rapidly deglaciating terrain, Glacier Bay National Park and Preserve and vicinity, Alaska and British Columbia Remote Sensing of Environment

The Glacier Bay area in southeastern Alaska and British Columbia, encompassing Glacier Bay National Park and Preserve, has experienced rapid glacier retreat since the end of the Little Ice Age in the mid-1800s. The impact that rapid deglaciation has had on the slope stability of valley walls and on the sedimentation of fans and deltas adjacent to fjords and inlets is an ongoing research topic. Using 3-year (2018 – 2020) Sentinel-1 datasets, and an automated time-series persistent scatterer interferometric synthetic aperture radar (PSInSAR) processing method, we detected landslides or delta subsidence at 27 sites within a vast 180 × 180 km remote region encompassing Glacier Bay proper. Most of the sites that we identified had not been previously identified. We categorized the hazard source areas that we identified into three general types:1) slow-moving landslides on steep rocky slopes not near ( > 2 km away from) present-day glacier termini (e.g., Tidal Inlet landslide), 2) slow- moving landslides directly adjacent to ( < 2 km away from), and associated with glacier thinning and retreat, and 3) subsidence of glacial outwash fan deltas. In categories 1 and 2, we observed 22 landslides moving at velocities ranging from 0.5 to 4 cm/yr. In category 3, we detected five fan deltas subsiding at velocities ranging from 0.5 to 6 cm/yr. Within our measurement error, these velocities were consistent during the monitoring period. Because acceleration was not observed, the issuance of warnings of imminent rapid failure is not warranted, however, continued remote monitoring is warranted. Our interferometric synthetic aperture radar (InSAR) results could be combined with other data sets including field observations, subaerial and submarine landslide inventories, bedrock fabric mapping from newly available light detection and ranging (lidar) data, and geologic maps to produce an inherent susceptibility map for landslides in bedrock and fan deltas. This map could be used to forecast susceptibility for both earthquake and climatically induced landslides.


Introduction
Glaciers, a climate-sensitive water resource, are shrinking rapidly and play a prominent role in increasing natural hazards, including landslides (e.g., Immerzeel et al., 2020;Wolken et al., 2021). In the early 21st century, worldwide, alpine glaciers have lost more mass than ice sheets in Greenland and Antarctica (Zemp et al., 2019). Of these, Alaska has the largest negative glacier mass changes among all polar, arctic, and alpine regions (Zemp et al., 2019;Ciracì et al., 2020;Hugonnet et al., 2021). Additionally, permafrost in Alaska is expected to degrade with warming temperatures (e.g., Pastick et al., 2015). The combination of receding glaciers and degrading permafrost has generated growing concerns about landslides induced by warming temperatures from climate change (Stevens et al., 2018;Hock et al., 2019;Wolken et al., 2021;Allen et al., 2022). As average winter temperatures rise closer to 0 • C ( Monahan and Fisichelli, 2014), a change in precipitation from snow to rain is also a potential factor that could increase the occurrence of landslides (e.g., Pavelsky et al., 2012). For example, in the next decade, winter temperatures for Gustavus, Alaska, the gateway town to Glacier Bay National Park and Preserve, are projected to change from 3 months with mean temperatures below freezing, to 1 to 2 months with mean temperatures below freezing (SNAP, 2022). An associated increase in precipitation is also expected (SNAP, 2022). Thus, over the long-term, precipitation will increase and occur more frequently as rain, snow elevations will rise, and snow will melt earlier and more often, which could lead to prolonged increases in ground water pore pressures and an increase in landslides (e.g., Baselt and Heinze, 2021). The loss of ice impacts rock slopes in a variety of ways (e.g., Evans and Clague, 1994;McColl, 2012) including causing a reduction in cohesion and an increase in pore-water pressure in fractures in rocks (e.g., Gruber and Haeberli, 2007), and debutressing (e.g., Deline et al., 2021) and/or thermomechanically damaging (Grämiger et al., 2017) rock slopes when it occurs in the form of receding glaciers. Rock-slope failures in mountainous regions can trigger cascading hazard events, particularly when they enter lakes or fiords, where outburst floods, debris flows, and tsunamis can be generated (e.g., Evans and Delaney, 2015;Haeberli et al., 2017;Higman et al., 2018).
Recent examples of cascading hazards that have attracted attention include tsunamis triggered by landslides entering fjords in Greenland (Svennevig et al., 2020) and Alaska (Higman et al., 2018); long runout debris flows triggered by rock/ice avalanches in Peru (Petley, 2020) and India (Shugar et al., 2021); and the threat of tsunamis from landslides in Prince William Sound, Alaska (e.g., Dai et al., 2020;Barnhart et al., 2021). In the United States, a sequence of recent landslide events (e.g., Taan Fjord landslide and tsunami, Higman et al., 2018;Lamplugh rock avalanche, Bessette-Kirton et al., 2018) and damaging and threating landslides (e.g., Pretty Rocks landslide, National Park Service, 2022a; Barry Arm landslides, Dai et al., 2020) in Alaska has resulted in the state becoming a focal point for research related to cascading hazards initiated by slope failures possibly caused by warming temperatures. One of the focus areas has been Glacier Bay National Park and Preserve (GBNPP) in southeast Alaska. GBNPP is visited by tens to hundreds of thousands of people in cruise ships and tour boats each year (about 64,000 in 2021 during the global COVID-19 pandemic, National Park Service, 2022b). Recent large, rapid landslides in GBNPP (e.g., Bessette-Kirton and Coe, 2016) have raised concerns that landslides may be increasing in size and mobility, possibly due to degrading permafrost and retreating glaciers from warming temperatures (Geertsema et al., 2013;Coe et al., 2018). One particularly obvious landslide within GBNPP, is the Tidal Inlet landslide (Wieczorek et al., 2007). This landslide is perched above the Tidal Inlet fjord and could cause a tsunami that would adversely impact a cruise ship route if the landslide failed rapidly into the water (Wieczorek et al., 2007).
Sedimentation rates on outwash plains, fans, and deltas along the margins of Glacier Bay are among the highest in the world (Cowan et al., 2010). This rapid accumulation of loose material may make fans and deltas especially susceptible to failure during an earthquake , which are common in GBNPP because of its proximity to a transform tectonic plate boundary. Submarine failures of fans and deltas in other coastal areas of Alaska have generated tsunamis that have caused destruction during earthquakes. For example, much of the damage caused by the 1964 Anchorage earthquake (Hansen, 1966) was from failures of fan deltas (e.g., Lee et al., 2006;Haeussler et al., 2014). Thus, GBNPP is potentially highly vulnerable to mass movements and cascading natural hazards that could adversely impact the public.
The ability to identify and track movement of mountain slopes, as well as any acceleration of such precursor movement prior to rapid failures, is a goal of landslide researchers and practitioners tasked with providing actionable information to local decision makers (e.g., Loew et al., 2017;Intrieri et al., 2018;Kristensen et al., 2021). For scientific and practical purposes, capturing and continuously monitoring landslides at a regional scale could be implemented by using remote sensing datasets. Accomplishing this task with remote sensing methods however, is difficult in steep terrain that is variably covered with snow and ice and is plagued by frequent cloud cover. These environmental and atmospheric conditions limit the use of synthetic aperture radar (SAR) (due to decorrelation) and optical imagery (due to poor visibility) that are commonly used to detect and track landslide movement under more favorable conditions (e.g., Intrieri et al., 2018).
In this paper, we use automated persistent scatterer interferometric synthetic aperture radar (PSInSAR) processing of Sentinel-1 synthetic aperture radar (SAR) data to detect and characterize very slow (16 mm/ yr to 1.6 m/yr, IUGS, 1995) and extremely slow (<16 mm/yr, IUGS, 1995) movement of mountain slopes, outwash plains, and fan deltas in GBNPP and its vicinity from 2018 to 2020. Throughout the paper, we use the general term "slow moving" to describe landslides that are moving at velocities <1.6 m/yr. Our results highlight active landslides and subsiding deltas where more detailed remote or in-situ monitoring could take place. Some of these areas of slow movement have not been previously noticed, and because of their proximity to drainages and fjords, they have the potential to generate cascading hazards that could impact humans far from the source areas. In addition to identifying the potential hazard source areas, we present evidence from recent (2019-2020) light detection and ranging (lidar) data and fieldwork in the summer of 2021 that support our hazard interpretations for some source areas. Finally, we explore possible regional conditions that lead to the identified hazards, and how our results could be used for continued remote monitoring and landslide hazard and risk assessments.

Geographical and geologic setting
Our study area, located between latitudes 58 and 60 • north, includes most of GBNPP in southeast Alaska and parts of Tatshenshini-Alsek Provincial Park in British Columbia (Fig. 1). It covers an area of ~180 × 180 km that consists of hundreds of kilometers of southeastern Alaska coastline; rugged, glaciated mountains in multiple mountain ranges; mountain, valley, and tidewater glaciers; outwash fan deltas; temperate rainforest; and glacial lakes, inlets, and fjords. The study area is bounded by the Takhinsha Mountains and the Chilkat Range to the east, the Pacific Ocean to the west, the Alsek Ranges to the north, and Cross Sound and Icy Strait to the south (Fig. 1). Most glaciers in the Glacier Bay region have had negative mass balance over the last 50 years. That is, glaciers have retreated and their equilibrium-line altitudes (ELAs) have risen from the mid-1970s to present (Pelto, 1987;Johnson et al., 2013). The total ice-covered area in GBNPP has decreased 15% in the last 50 years (Loso et al., 2014). Consequently, some tidewater glaciers have converted to land-terminating glaciers (Loso et al., 2014;National Park Service, 2021a). As glaciers retreated, terminal moraines were deposited. Glacial meltwater transported glacial sediments downstream and deposited them on outwash plains and fan deltas at junctions with fjords.
The study area is at an active tectonic plate boundary between the Pacific and North American Plates. This plate boundary is marked by the strike-slip Fairweather fault, which accommodates dextral slip of about 43 mm/yr (Elliott et al., 2010). The Fairweather Range contains some of the highest coastal mountains in North America (e.g., Mount Fairweather is the highest at 4671 m). Some of the lowest latitude tidewater glaciers in the world are sustained by moist air blown from the Gulf of Alaska, which rises over high mountain ranges and provides precipitation to the area in the form of snow and rain (National Park Service, 2015). The landscape of the study area is a direct result of cycles of glacial advance and retreat. Glacier Bay was completely filled with ice as recently as 1750 (Powell, 1980). Glacial retreat since 1750 has created Glacier Bay proper (a fjord with numerous arms and inlets) and caused rapid uplift (up to 3 cm/yr) from viscoelastic rebound (Larsen et al., 2005).
The geology of the region is dominated by accreted terranes (Fig. S1; Brew et al., 1978;Brew, 2008;Smart et al., 1996;Wilson et al., 2015). The Yakutat terrane, which is west of the Fairweather fault, is the youngest of all the terranes in the region and consists of a Late Cretaceous accretionary complex overlain by Eocene marine and continental clastic rocks (Nelson and Colpron, 2007). The Chugach terrane, which is bounded by the Fairweather fault in the west and the Border Ranges fault in the east (Fig. S1), consists of sedimentary and metamorphic rocks that formed ~70-140 million years ago. The Alexander terrane, which is east of the Border Ranges fault, is the oldest terrane in the region and is composed of ~500 million year-old rocks (e.g., Silurian turbidites and limestone) formed by a volcanic arc system, as well as other sedimentary rocks, and a 90-120 million year-old Cretaceous granite (Nelson and Colpron, 2007;National Park Service, 2021b).
The Glacier Bay area is susceptible to rock slope failures triggered by earthquakes and climatic conditions (Coe et al., 2019). Rocks in the area have been weakened by tectonic stresses and earthquake deformation, and by damage from repeated cycles of glaciation (e.g., Grämiger et al., 2017). Warming from climate change is likely making the terrain more susceptible to landslides because of degrading mountain permafrost and glacial thinning and retreat, which causes bedrock rebound and can alter pore water pressures in valley walls (e.g., Oestreicher et al., 2021). Cascading hazards are possible (e.g., landslide-generated tsunamis), if subaerial or submarine landslides are large and fast moving when they enter or displace water. These tsunamis can impact areas far afield from landslide source areas. One precedent for a landslide generated tsunami in GBNPP is the 30 M m 3 subaerial rockslide and subsequent tsunami triggered by a M 7.8 to 8.3 earthquake near Lituya Bay in 1958. This landslide generated a giant wave runup (524 m in height) and resulted in multiple casualties near the mouth of Lituya Bay about 11 km away from the landslide source (Miller, 1960;Mader and Gittings, 2002). The narrow geometries of inlets and fjords in the GBNPP area would likely exacerbate tsunami impacts (Wieczorek et al., 2007).

Methods
There are challenges in capturing ground deformation in the Glacier Bay region from spaceborne SAR data and interferometric SAR (InSAR) methods. Most of the region is covered by snow from November to May which causes InSAR observations to lose coherence. Rapidly moving landslides (e.g., in our study, movement >~ 2.8 cm between SAR data sets) could also not be detected by InSAR observations. This lack of detection is because the interferometric phases are constrained within 2π modulo (phase periodic with a 2π radian period), and rapid surface motion exceeding a half radar wavelength will result in a large spatial gradient of phase changes and thereby loss of coherence . Additionally, throughout the year, much of the region is covered by moving glaciers, which create problems for coherence and for  (2018)(2019)(2020). Ground deformation includes slow-moving landslides distal to present-day glacial termini (orange placemarks); slow-moving landslides associated with ice withdrawal or movement of debris along the edge of glaciers (cyan placemarks); and subsiding fan deltas near the fronts of glaciers (purple placemarks). The LOS (line-of-sight) deformation rates including deformation areas highlighted in Figs. 2, 4, and 6 range from − 50 to 50 mm/yr. Positive values mean that distances from the satellite to the ground surface have become shorter and negative values mean that distances have become longer. Orange line shows the boundary of GBNPP. Light blue line is the boundary between the United States and Canada. Red and navy lines are Quaternary (Koehler et al., 2012) and neotectonic (Brew and Ford, 1985;Plafker et al., 1994;Koehler et al., 2013) faults in the area, respectively. Black triangles indicate the locations of GPS stations within our study area. See the KMZ file in the supplemental files for polygons outlining the approximate area deforming at each placemark. The KMZ file can also be accessed via USGS data release (Kim et al., 2022).https://smu.box.com/s /cvjqmbc70yxk74ump0qhegaexmi6owth (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) accurate interpretation of the causes of observed deformation. SAR scenes acquired during the 5-month, relatively snow-free period from June to October can have strong atmospheric influences which make subtle ground deformation difficult to decipher. A few InSAR pairs are not enough to identify deformation from low coherent atmospherecontaminated pixels. Multi-temporal InSAR methods are better suited for that purpose (Berardino et al., 2002;Ferretti et al., 2001). Therefore, for our work, we applied the PSInSAR method to Sentinel-1 datasets acquired from June to October, with a 12-day interval between acquisitions. PSInSAR method exploits the phase measurements at persistent scatterers (PSs) that represent point-like responses from radar targets that are characterized by high coherence and phase stability over the entire period of observation. This technique can overcome decorrelation and atmospheric anomalies through spatiotemporal filtering, in contrast to the conventional InSAR processing which depends on a limited number of interferograms (Hooper et al., 2004;Hooper, 2008;Lu and Kim, 2021).
Sentinel-1 images are available for the Glacier Bay region from 2018 to present (2022) from the European Satellite Agency (ESA) Copernicus Open Access Hub and the Alaska Satellite Facility (ASF)'s Vertex platform. Given the defined area of interest (Glacier Bay), all Level-1 Sentinel-1 SLC (single look complex) datasets for summers (~June--October) from 2018 to 2020 were downloaded from the data archives and processed together (i.e., they were not processed as separate years) during PSInSAR processing. Only the bursts of the Sentinel-1 Interferometric Wide (IW) swath that include the interest area were selected, merged, and processed. We then used the enhanced spectral diversity (ESD) method to generate a precisely co-registered stack (Prats-Iraola et al., 2012). Sentinel-1 SAR data included descending (P145) and ascending (P50) tracks that cover almost the entirety of GBNPP. All preprocessing of data, from data transfer to co-registration, was automatically implemented. The co-registered stack was separated into multiple patches with 3000 × 2000 pixels. Differential interferograms were generated using the Copernicus digital elevation model (DEM) with global coverage at 30 m resolution (GLO-30). The patch interferograms were processed using the PSInSAR method at full spatial resolution, with processing steps that included PSs selections, DEM error correction, phase unwrapping, and atmospheric filtering (Hooper et al., 2004;Hooper, 2008). Interferograms over Glacier Bay obtained during summer months can have strong tropospheric effects (vertical stratification), so regression analysis between elevation and unwrapped phases was also applied to reduce such effects (Bekaert et al., 2015). The patch processing has pros and cons. Because large areas of ~180 × 180 km were processed in parallel, processing time can be greatly reduced using patch processing (from weeks to a day or so). Also, iterative PSInSAR processing to decipher error/noise components through spatiotemporal filtering and to extract deformation components (Hooper et al., 2004;Hooper, 2008) were readily applied to each patch. However, there can be discontinuities between patch outputs. Thus, to minimize such discontinuities, overlap between patches was averaged in the merge process. Calculation of deformation rates and identification of landslides was carried out on the merged outputs. There are five Global Positioning System (GPS) stations within our study area (black triangles in Figs. 1 and S2) that could potentially be used to compare with our PSInSAR measurements. However, two stations (ALSC, BEUT) have become inactive, and are not currently providing GPS measurements, one station (AB43) is located on an isolated island, "Cape Spencer", which does not have neighboring PS points, and the other station (QUIC) was missing year-long measurements during our study period. Consequently, only a single GPS station (BMCP) is available for comparing our PSInSAR measurements with a reference point on the ground. Comparison between the two is in good agreement, with the InSAR line-of-sight (LOS) measurements falling within tighter clusters on the same trend as the GPS measurements ( Fig. S3(a), (b)). Based on geographical characteristics of narrow inlets, steep slopes, and valleys, and validation results at only a single GPS station ( Fig. S3; see supplementary describing how GPS measurements are converted to those in the LOS direction for comparing GPS and InSAR measurements), we instead used the averaged value of each patch as a reference value for phase unwrapping and calculation of deformation rates, while excluding outliers of 20% in the upper and lower bounds.

Results and interpretations of processes
Using the PSInSAR method, we identified and confirmed ground deformation in the Glacier Bay region using the descending and ascending tracks, respectively (see Fig. 1 and Fig. S2 for results from descending and ascending tracks, respectively). Because the descending track (P145) has better coherence over most potential landslide areas and contains more PS points than the ascending track (P50), we used the PSInSAR results from the descending track as a primary means to identify ground deformation. We then used results from the independent ascending track to confirm the deformation. In locations where ascending track PS points were too sparse for verification, we used optical satellite images in Google Earth for verification. Moreover, when we observed a strong linearity in time-series measurements from the descending track, we marked the locations as active deformation despite an absence of PS points from an ascending track. The areas of ground deformation that we identified can be categorized as follows: 1) slowmoving landslides distal (> 2 km away from) to present-day glacial termini (orange placemarks in Fig. 1), 2) slow-moving landslides associated with ice withdrawal or movement of debris at or near glacier margins (< 2 km away from glacier termini, cyan placemarks in Fig. 1), and 3) subsiding glacial outwash plains and deltas (purple placemarks in Fig. 1). Because >20 deformation areas were identified, in the following sections we highlight only a few representative examples of deformation from each category. In some of these cases, we provide background information needed to place observed movement into known geological and historical context, as well as our interpretive insights into possible factors responsible for the observed movement.

Slow-moving landslides distal to (> 2 km away from) present-day glacial termini
Slow-moving landslides not near present-day glacial termini may be deforming because of seasonal snow melt, rainfall, or mountain permafrost degradation due to the warming effects of climate change. Weakening of rocks by tectonic processes, permafrost degradation, and cycles of glacial erosion and debuttressing (meaning glacial unloading or loss of lateral glacial support, e.g., Deline et al., 2021) has likely made parts of the Glacier Bay region susceptible to landslides. Using PSInSAR, we identified four slow-moving landslides in the Glacier Bay region distal to present-day glacial termini (Fig. 1).

Tidal Inlet landslide
The most concerning landslide in this category is in Tidal Inlet on the east side of the West Arm of Glacier Bay (Fig. 1). It is the most concerning landslide because a rapid failure of this landslide into the water of the narrow inlet could produce a wave that could detrimentally impact marine activities in the West Arm (Geist et al., 2003;Wieczorek et al., 2007;National Park Service, 2018). Wieczorek et al. (2007) measured movement of 3-4 cm/yr at the Tidal Inlet landslide during a 2year period (2002)(2003)(2004) using repeat GPS campaigns. Since 2004, landslide movement has not been monitored in Tidal Inlet. Thus, any deceleration or acceleration of possible movement since 2007 remains unknown, with the latter being a potential indicator of imminent catastrophic failure (e.g., Loew et al., 2017). The landslide is underlain by the late Silurian to middle Devonian Pyramid Peak Limestone (Dsp; Rossman, 1963;Fig. 2) and the late Silurian Tidal Formation (Stg, Fig. 2), which is exposed near the shoreline of Tidal Inlet. Most of the active landslide has been mapped as weathered Quaternary surficial deposits (Qs, Fig. 2(a); Wieczorek et al., 2007). The landslide faces south-southwest and has a slope angle that ranges from 30 to 40 • . The headscarp has a slope angle of ~45 • , and the upper portion of the landslide between the headscarp and a major internal scarp contains rotational blocks and escarpments that are covered with thick vegetation ( Fig. 2(a), (b)). These observations indicate that the landslide is unstable, but its potential for rapid failure is unknown.
PSInSAR results from both descending ( Fig. 2(b)) and ascending ( Fig.  S4(a)) tracks show that the bare ground downslope from the major internal scarp is the most active area of the Tidal Inlet landslide with − 14.3 ± 1.5 mm/yr and − 9.5 ± 0.7 mm/yr of movement (~2.5 cm of total movement for ~2.5 years), in the LOS directions for the descending track ( Fig. 2(c)) and ascending track ( Fig. S4(b)), respectively. Negative deformation rates shown in Fig. 2(b) indicate that the most dominant direction of landslide motion is downward (vertically) and southward (horizontally), based on SAR imaging geometry. An examination of observed movement with respect to landslide structures mapped from 2019 lidar data ( Fig. 3(a)) indicates that the most active area is downslope from the major downhill facing internal fault scarp that "bites" into a prominent vegetation-covered topographic bench within the landslide body ( Fig. 3(a), (b)). Most other landslide structures within the bench are uphill-facing normal-fault scarps, which are typical of slow ridge-spreading movement. Such movement is often found along the flanks of valleys formerly occupied by glaciers. Ridge-spreading movement and uphill facing scarps have been documented in areas with both very recent glacier retreat (e.g., within the last 20 years, Dai et al., 2020;Coe et al., 2021), as well as in areas where glaciers have been gone from the valley for 10,000 years or more (e.g., Varnes et al., 2000). Results from multiple GPS campaigns above the major downhill facing scarp and near other internal uphill facing scarps in 2002, 2003, and 2004 showed horizontal southerly movements of up to 79 ± 15 mm between July 2002 and August 2004, for a mean rate of horizontal motion of ~40 mm/yr (Wieczorek et al., 2007). However, Wieczorek et al. (2007) did not detect vertical motion within their reported limits of uncertainty. Our results show that the center of the landslide near the major downhill facing scarp is still creeping slowly. We could not use PSInSAR to measure movements between the major internal scarp and headscarp ( Fig. 2(a)) because of the steep slope angle and thick vegetation. Other parts of the slope adjacent to the Tidal Inlet landslide, which are cut by extensive uphill facing scarps (Fig. 3(a)), showed no detectable movement between 2018 and 2020. In addition, deformation rate in vertical (Fig. S4(c)) and east (Fig. S4(d)) directions and their uncertainties (Fig. S4(e), (f), respectively) can be estimated based on three-dimensional (3D) surface deformation and PSInSAR measurements using geometric parameters (incidence and heading angle) (Samsonov and d'Oreye, 2012;Kang et al., 2021), although deformation in north-south direction cannot be properly calculated due to large  Fig. 2(b)). As stated in Fig. 1, positive values mean that distances along look direction (arrow in Fig. 2(b)) from the satellite to the ground surface have become shorter and negative values mean that distances have become longer. uncertainty from near-polar orbits of the SAR satellite and its right-only look direction (Wright et al., 2004). However, the estimated deformation rate in vertical and east-west directions does not add much more information for discriminating landslides on the slope, which was the primary purpose of our study. Therefore, our subsequent landslide identification in GBNP (i.e., described in the text that follows) relied on examining the PSInSAR results from descending track and validation with those results from the ascending track.

Buckwell W4 landslide
Another example of a slow-moving landslide distal to a present-day glacier terminus occurs about 1 km southwest of Buckwell W4 Peak in Canada (59.42 • N, 136.85 • W), the highest mountain (2721 m) in the Alsek Ranges (Fig. S4(g), (h)). The region is part of Tatshenshini-Alsek Park, British Columbia, Canada, with surface water flowing toward the Tkope River. This landslide had not been previously identified and occurs near the contact between Cenozoic intrusive and Paleozoic-Mesozoic sedimentary rocks (Chorlton, 2007), which are exposed at the site. More detailed geological information is not available due to the site's remote location. The measured peak LOS deformation at the location of the white dot in Fig. S4(g), (h)) was − 31 ± 4.2 mm/yr (descending track) and 12 ± 4.8 mm/yr (ascending track). A deformation time series is shown in Fig. S4(i)). Because the sign of measured deformation is opposite from two different look directions, downslope movement from the steep (~20-40 • ) slope is dominant in the landslide.

Slow-moving landslides at or near glacier margins (<2 km away from glacier termini)
Tidewater and lake-terminating glaciers terminate in the ocean or glacial lakes and discharge ice into the water. These glaciers are particularly sensitive to climate change. Warming temperatures can trigger their retreat, which once initiated, can be sustained through a series of positive feedbacks between retreat, acceleration, and dynamic thinning, and can persist for long time periods (>~20 years; Meier and Post, 1987;Post and Motyka, 1995;Pfeffer, 2007). From PSInSAR results, we observed slow moving landslides on valley walls at or near (<2 km) the termini of three glaciers: the lake terminating Melbern Glacier (Fig. 4(a)) and Grand Plateau Glacier (Fig. 4(b)), and the tidewater Gilman Glacier (Fig. 4(c)). At representative locations on these slopes shown with white dots in Fig. 4(a-c)), the measured LOS deformation was − 35.3 ± 4.1, − 29.8 ± 3.7, and − 29.9 ± 2.5 mm/yr, for Melbern, Grand Plateau, and Gilman glaciers, respectively. In general, these deformation rates are 2 to 3 times greater than the deformation rates of the Tidal Inlet landslide, which was last near a glacier terminus in the 1870s (Powell, 1980). However, glacier-distal landslides such as the Tidal Inlet landslide may have had periods of faster deformation in the past when they were adjacent to retreating glacial termini. A time series of slope movement at Melbern, Grand Plateau, and Gilman glaciers is given in Fig. S5(a).

Melbern landslide
Melbern Glacier (Fig. 4(a)) is in Canada but is contiguous with the much larger Grand Pacific Glacier to the south, which terminates at the north end of Tarr Inlet in GBNPP. Melbern Glacier terminates in Lake Melbern and has undergone extreme retreat for the past 40 years. Since the 1970s, Melbern Glacier retreated about 15 km to its present position and thinned about 300-600 m . Consequently, Melbern Glacier and the two glaciers to the north (Konamoxt and Tikke Glaciers) separated, an ice dam between Melbern Glacier and Konamoxt Glacier melted, and continued glacier melting resulted in the formation of Lake Melbern, which did not exist prior to the mid-1970s. SAR intensity images (ERS-1/2, ENVISAT, ALOS, Sentinel-1) acquired from the early 1990s to present are helpful for delineating the main terminus locations of Melbern Glacier through time. From 1991 (cyan line in Fig. 4(a)) to 2020 (gray lines for other years in Fig. 4(a)), the glacier terminus retreated >5 km and Lake Melbern continuously expanded ( Fig. S5(b)). Our observations of ground movement on the valley confirms a report from Clague and Evans (1994) about "destabilization of formerly ice-covered, steep rock slopes", although our 3-year record is not sufficient to interpret its cause. The slope instability may be associated with rock damage caused by repeated cycles of glaciation and rapid deglaciation (Cossart et al., 2008), as well as glacial debuttressing caused by glacial thinning during recent ice retreat, and additional rock damage caused by rapid glacial retreat.

Grand Plateau landslide
Grand Plateau Glacier has two major distributary tongues: one located in Alsek Lake that discharges water through the Alsek River into the Pacific Ocean, and the other in Grand Plateau Lake. Glacier termini at both locations have retreated substantially (Loso et al., 2021): the terminus in Alsek Lake retreated as much as 5 km from 1996 to present and the terminus at Grand Plateau Lake retreated as much as 4.5 km for the same period ( Fig. S5(b)). As this glacial retreat continues, Loso et al. (2021) predicted that Alsek and Grand Plateau Lakes will connect in the future, and the riverine ecosystem, landscape, and river path will consequently change. Laser altimetry measurements from 2017 to 2020 indicate thinning of termini at rates of up to 10 m/yr (Loso et al., 2021). Similar to Melbern Glacier, such rapid glacial thinning can cause debuttressing effects on valley walls along the route of glacier outflow. We identified multiple locations of slow deformation along the flanks of the retreating glacier (Figs. 1, 4(b)), and those four neighboring landslides have similar rate (10-30 mm/yr) of downslope movement. Structures caused by movement of these landslides are not prominent in a relatively low resolution 5-m IfSAR (Interferometric SAR) DEM but appear to be mostly downhill-facing scarps.

Gilman landslide
Gilman Glacier is different from the rapidly retreating Melbern and Grand Plateau Glaciers because it is a tidewater glacier that merges with Johns Hopkins Glacier and calves into Johns Hopkins Inlet in GBNPP (Figs. 1, 4(c)). Both Johns Hopkins and Gilman Glaciers have a recent positive mass balance, though Gilman Glacier has a smaller net value (Johnson, 2012). Johns Hopkins Glacier is slowly advancing and the terminus of Gilman Glacier has repeatedly retreated and advanced in the last 64 years (McNabb and Hock, 2014). As such, the tongues of both glaciers continue to coalesce and separate through time. Although there was a report that lower reaches of Johns Hopkins Glacier were thinning based on IceBridge measurements (NASA, 2021), and there is no direct measurement on the thickness of the Gilman Glacier, it seems likely that it is neither thinning nor thickening greatly, based on the observed slightly positive mass balance measurements (Johnson, 2012;Johnson et al., 2013). Therefore, the cause of the landslide that we identified near the terminus of Gilman Glacier (Fig. 4(c)) is more likely related to erosion and rock damage to the valley wall caused by repeated cycles of glaciation and movement of the glacier.
Landslide structures mapped from 2019 lidar data (Fig. 4(e)) are reflective of observed landslide movement (compare Fig. 4(c), (e)). Unlike structures at the Tidal Inlet landslide, all landslide fault scarps at the Gilman Glacier landslide are downhill facing. The most prominent downhill facing scarp forms the primary landslide headscarp and trends from southwest to northeast, ending at the ridgeline at an elevation of about 1140 m above sea level (Fig. 4(e)). Secondary scarps indicate that the landslide is complex. Many of the scarps probably reflect the presence of pre-landslide geological structures (dikes, foliation, fractures, etc.) that have influenced the location and style of subsequent landslide movement. The lack of any detectable upward movement or landslide thrust faults at the base of the slope suggests that there has either been very little movement of the landslide since inception, or that the landslide toe occurs at elevations lower than the surface of the glacier or was eroded by the glacier.

Subsidence in outwash plains and deltas
At least three glaciers (Carroll, Lituya, Brady, Fig. 5) have experienced substantial retreat for last 35 years and have extensive outwash deltas ( Fig. S6(b)). The debris covered Lituya Glacier has had the most retreat during this period (see Fig. S6(b)). In 1984, it was a tidewater glacier and was calving into Gilbert Inlet. By 2004, the glacier had transformed into a land-terminating glacier with a calving front that was completely disconnected from Gilbert Inlet. The now grounded glacier is still retreating, but its pace has slowed. From 1984 to 2004, Lituya Glacier retreated about 2 km, but has not visibly retreated since 2004 (Fig. S6). Based on previous observations, both Carroll and Brady Glaciers were tidewater glaciers in the early 1900s and late 1800s (Carlson et al., 1999 andPelto et al., 2013, respectively) before transitioning to land-terminating glaciers.
PSInSAR results reveal that the interior outwash plains of the three glaciers are undergoing rapid subsidence at locations shown with white dots in Fig. 5(a-c). Measured LOS deformation rates are − 56.2 ± 3.8 mm/yr at Carroll Glacier ( Fig. 5(a)), − 24.2 ± 3.7 mm/yr at Lituya Glacier ( Fig. 5(b)), and − 15.7 ± 1.8 mm/yr at Brady Glacier (Fig. 5(c). For time-series estimates of movement at each site, see Fig. S6(a)). The outwash plains contain sediment that ranges from boulder-to clay-sized particles (Fig. 6). One hypothesis about the cause of subsidence in the plain is that the moist, loose sediments consolidate in time, whereas pore water in soil layers diffuses away to downstream regions of the outwash plains (Terzaghi, 1925). The geomechanical process occurring in the outwash plain may be similar to soil compaction that can be seen in tailings dams (Hu et al., 2017), reclaimed lands (Wu et al., 2020), or regions with groundwater-withdrawal (Qu et al., 2015). Another hypothesis is that ice detached from the glacier (called "dead ice"; Schomacker, 2008) is rapidly buried in the outwash plain and gradually melts leading to subsidence of the surficial sediments. The process is identical to that which forms kettle lakes that are often found in the glacial outwash (Bennett and Glasser, 2009). Although the subsidence can result from a combination of two hypothesized processes, one process may dominate at each of the three outwash plains. The outwash plain at Carroll Glacier contains a kettle lake (Fig. 6(a)). There is an active meltwater stream along the western margin of the outwash plain, but on the east side of the plain, there is field and topographic evidence that subsidence and back tilting of sediments has caused former stream paths to now drain in an upslope (northerly) direction ( Fig. 6(a), (b)). Additionally, there are curvilinear normal fault scarps up to about 30 cm in height that are visible in sediments in the field and in high-resolution lidar data on the southeast side of the kettle lake ( Fig. 6(a), (c)). These scarps mimic the general shape of the lake and are the direct result of discrete deformation caused by subsidence. Therefore, melting of buried ice is likely the dominant process causing subsidence of the outwash plain at Carroll Glacier ( Fig. 5(a)). In contrast, outwash plains at Lituya and Brady Glaciers have well developed meltwater streams and sediment consolidation appears to be the more dominant process resulting in subsidence ( Fig. 5(b), (c)) but the exact cause of the subsidence warrants further investigation by field surveys. The lateral margins of outwash plains show evidence for an increase in elevation of ~5-8 mm/yr, probably from deposition of sediments along meltwater streams (Fig. 5  (a-c)). However, in more interior portions of the plains, subsidence from consolidation and melting ice has completely overwhelmed this growth from sediment deposition.

Discussion of implications
Our results identified new landslides and fan delta subsidence and establish baseline rates of movement for the Glacier Bay region. We detected 22 landslides in the study area, with the smallest landslide having an area of about 54,900 m 2 , and the largest having an area of about 3,375,200 m 2 (see data in Kim et al., 2022). The landslides moved at velocities ranging from 0.5 to 4 cm/yr between 2018 and 2020 ( Fig. 7  (a), (b)). Landslides proximal (within <2 km) to glacial termini and distal (> 2 km away from) to glacier termini have similar and overlapping velocities (Fig. 7). All the velocity measurements are in the LOS direction, and it is difficult to fully resolve 3D deformation from sparsely distributed, only right-looking observations. Although the LOS measurements cannot fully restore the actual slope motions, it is still sufficient to identify landslide motions. All of the observed landslide velocities are classified as very slow (16 mm/year to 1.6 m/year) to extremely slow (<16 mm/yr) by the International Union of Geological Sciences (IUGS) Working Group on Landslides (IUGS, 1995). During the 2018 to 2020 time period, none of the landslides accelerated, but instead moved at consistent velocities, at least when measured on an annual cycle and accounting for our range in measurement error (Fig. 7).
Observations of areas of identified landslide movement in Google Earth suggest that the landslides are predominantly failures within bedrock, although field investigations would be useful to verify this interpretation. Geologic units and headscarp elevations extracted from existing regional geologic maps and DEMs reveal that 68% of all landslides occur within flysch and flysch containing volcanic rocks (Fig. 8a) and that 45% of all landslide headscarps occur between elevations of 600 to 800 m above sea level (Fig. 8b). This elevation band roughly corresponds with the lower elevation limit of mountain permafrost in GBNPP (~500 m above sea level) as mapped by Gruber (2012), suggesting that permafrost degradation may play a role in landslide occurrence in GBNPP (as previously hypothesized by Coe et al., 2018). Within 2 km of glacier termini, the number of landslides in flysch and flysch with volcanic rocks is about 78% (Fig. 8c), suggesting that flysch may be particularly susceptible to slow landslide movement during glacial thinning and retreat. In support of this statement, the recently discovered ~500 Mm 3 Barry Arm landslide (Dai et al., 2020;Coe et al., 2021) in Prince William Sound, Alaska is also in a flysch. For the four landslides in GBNPP that are >2 km away from glacier termini, there is more geologic diversity, with two (50%) occurring in limestones and marbles, one (25%) in flysch, and one (25%) in granodiorite (Fig. 8e). The landslides with the highest elevations (1400 to 1600 m above sea level) in the study, also are >2 km away from glacier termini (compare Fig. 8(b), (d), and (f)). Five fan deltas in the study area subsided at velocities ranging from 0.5 to 6 cm/yr (Fig. 7(c)). As with landslides, subsidence velocity did not accelerate substantially during the monitoring period, but instead maintained a consistent slow rate of speed within our error bounds (Fig. 7). From a hazard point of view, these results necessitate answers to the following questions: 1) are the measured movement velocities cause for urgent action (e.g., issuance of a warning) by land managers? 2) should these areas be monitored in the  . 7. LOS speed (absolute value of velocity) of identified deformation in the study area, which was averaged within a radius of 100 m of each location in slow-moving landslides (a) distal to present-day glacial termini, and (b) proximal to present-day glacial termini, associated with ice withdrawal or movement of debris along the edge of glaciers, and (c) subsidence of outwash plains and deltas. Note that due to the large uncertainty for estimated yearly velocity (mm/yr) in a single year (observation for ~ 5 months a year), multi-year velocity was calculated. See the KMZ file in the supplemental files for polygons outlining the approximate area deforming at each placemark. The KMZ file can also be accessed via USGS data release (Kim et al., 2022). future, and, if so, how? and 3) how chould the results be incorporated into subaerial landslide hazard maps?
Addressing Question 1, in general, slow steady movement is cause for continued watchfulness, but not for alarm or issuance of a warning that a catastrophic, fast-moving landslide is imminent. Regarding the 22 landslides that we have identified, there are many examples in the literature of slow-moving bedrock landslides, that have not and may never fail catastrophically (e.g., Cascade landslide in Washington (Hu et al., 2016); Three Bears landslide in California (Liu et al., 2019); slowmoving landslides on the U.S. West Coast ). At least through June 2022, the Tidal Inlet landslide is a good example of one of these landslides. Work by Wieczorek et al. (2007) suggests that it has probably been slowly moving since the early 1900s and our results indicate that ongoing movement is also slow (Fig. 2(b)). For warning purposes, landslide acceleration is one of the main indicators that a catastrophic rapid failure could be imminent (e.g., Hermanns et al., 2013;Loew et al., 2017). However, there are also examples of landslides that have periods of acceleration without catastrophic failures. For example, the Rattlesnake Hills landslide in the Miocene Columbia River Basalt in Washington State accelerated up to about 7 cm/day during a 3month period in late 2017, then slowed to steady movement, and then decelerated (Washington Department of Natural Resources, 2019). Another example is the Veslemannen rockslide in Norway, which has had at least 16 episodes of rapid acceleration during 5 years of continuous monitoring before catastrophically failing in 2019 (Kristensen et al., 2021). There are geotechnical guidelines available that can be (c) and (d) characteristics for landslides <2 km away from glacier termini; (e) and (f) characteristics for landslides >2 km away from glacier termini. used to assess the potential for rapid failure of landslides that are slowly moving (e.g., Fell et al., 2007), but these guidelines require knowledge of controlling factors (e.g., the orientation of bedding or foliation, mechanism of failure, and geometry and roughness of the slip surface) that generally require field investigations. For the sites that we have identified, field investigations of each site, even if limited in scope, would be useful in determining these factors. However, based on the data that we currently have available, all of the landslides that we detected in the Glacier Bay region are slowly creeping downslope, but are not accelerating. Therefore, continued remote monitoring of these areas is warranted, but the issuance of a warning or warnings of imminent rapid failure of the landslides is not currently warranted.
The issuance of a warning or warnings for the five subsiding fan deltas are also not currently warranted. Geotechnical guidelines such as those presented by Fell et al. (2007) might also be applied to failures in fan delta sediment. However, we are not aware of previous work that has quantified subsidence of fan deltas in recently glaciated environments, and the relation between steady subsidence not related to a specific trigger, and the potential for landslide failures or liquefaction of fan deltas during an earthquake (for example, liquefaction failure in the Fraser River delta, Canada (Chillarige et al., 1997)) is unclear. Recent work by Avdievitch and Coe (2022) did not reveal any submarine or subaerial landslide scars on the fronts of fan deltas where we detected subsidence, but this is likely because the deposition of outwash sediment rapidly fills and covers these scars. More research is needed on the hazard and risk posed by subsiding fan deltas. For example, a hypothesis that could be tested by the occurrence of a large (e.g., M > 7) earthquake in the Glacier Bay region is: Fan deltas that are subsiding during tectonically quiescent time periods are more susceptible to failure during large earthquakes than are fan deltas that are not subsiding during quiescent periods.
In terms of Question 2, regarding the need for continued monitoring, the answer is "yes", in areas where we have baseline data, continued monitoring would be warranted. Ground-based monitoring generally is not feasible to cover all the landslides in the area because of the difficulty and expense of accessing and maintaining stations in remote wilderness locations, but continued satellite monitoring and targeted ground-based studies are warranted, especially for landslides that have potential to enter bodies of water and generate tsunamis. Costeffective operational monitoring of landslides with satellite data may be best achieved in regions with numerous landslides that pose a threat to the public. Currently, no academic, private, or government group conducts routine remote, regional landslide monitoring in the U.S., but government groups in other parts of the northern hemisphere (e.g., Greenland and Norway) have begun some routine monitoring for landslide movement (e.g., Lauknes et al., 2010 in Norway;and Svennevig et al., 2019 in Greenland). The tsunami threat posed by the Barry Arm landslide in Alaska prompted the U.S. Geological Survey to begin systematic InSAR monitoring of the site (Schaefer et al., 2020). Additionally, as we demonstrate in this paper, there are satellite data freely available that can be used for routine monitoring, and the expected 2023 launch of the joint U.S. National Aeronautics and Space Administration (NASA) and Indian Space Research Organisation (ISRO) Synthetic Aperture Radar Mission (NISAR) is expected to provide additional multiband SAR data for 70% of the Earth's ice-covered surfaces (see https://nisar.jpl.nasa.gov/). For efficiency, an automated PSInSAR method like the one demonstrated in this paper could be used to routinely process data to monitor large areas.
Use of an automated PSInSAR technique in steep terrain with snow and ice cover has limitations. First, there are, generally, sparse PS points on steep slopes and it is difficult to map the extent of displacements. Second, as the processing time periods become longer, the number of PS points diminishes due to the temporal decorrelation. After many trials, we found that, for PSInSAR processing, the use of Sentinel-1 datasets that were acquired over a period of three consecutive summers worked best for the Glacier Bay area. Because we only used data from summer months, we could be missing seasonal acceleration and deceleration patterns that would be useful for determining and isolating climatic, glacial, or permafrost variables that influence movement. Moreover, for tracking the long-term motion of landslides and potentially capturing short-term acceleration and deceleration, the use of long-term (>5 years) SAR images is necessary. The additional coverage provided by upcoming Sentinel-1C/D (Torres et al., 2021) and NISAR datasets should be helpful for this task. Third, fast-moving deformation (i.e., deformation much greater than 2π (~2.8 cm for Sentinel-1) between SAR data sets) cannot be measured due to phase decorrelation and the adverse characteristics of interferometric phases constrained within 2π modulo . Fourth, over heavier vegetated southern GBNPP, the coherence of the C-band InSAR is not maintained, even during summer months, and we could not identify landslides in this area. In addition, InSAR cannot measure the travel distance of landslides, rapid rock avalanches or debris flows, and other types of rapid slope collapses. In these cases, only the decorrelation information in the co-event InSAR image might provide a hint to the possible extent of landslide travel and inundation.
For Question 3, regarding how to incorporate InSAR created landslide inventories and monitoring results into regional landslide hazard maps, there are likely several different ways that this could be accomplished. However, clearly if a landslide is moving, it should be ranked as a high hazard zone on landslide susceptibility or probability maps. Ideally, InSAR results could be combined with other data sets (e.g., landslide inventories, lidar, geology) typically used to forecast landslide susceptibility. For areas with existing susceptibility maps, InSAR results could be used to test, validate, or refine the maps. The Norwegians have developed an approach to integrate geomorphic, structural, and movement data (Hermanns et al., 2013) into a quantitative landslide hazard assessment scheme. This, or a similar approach, could also work for areas in Alaska. In the Glacier Bay region, newly available lidar data are currently being used to create a map of the region that portrays bedrock structures (i.e., the bedrock fabric), landslide scarps in bedrock, and landslide deposits. In the near term, the InSAR results presented in this paper could be combined with the lidar derived inventory map, and field-based rock mass quality measurements and landslide and fan site visits, to produce an inherent susceptibility map for landslides in bedrock. By inherent susceptibility, we mean that the map could be used to forecast susceptibility for earthquake or climatically induced landslides.

Conclusions
Using medium-resolution, freely available Sentinel-1 SAR imagery and a time-series PSInSAR processing technique, we have demonstrated that we can detect landslides and subsidence using summer and early fall season (June-October) SAR imagery in the snow-and ice-covered Glacier Bay region in southeast Alaska. Using automated PSInSAR processing, we detected 22 slow moving landslides on rocky steep slopes both near and distal from glacier termini, and slow subsidence in five outwash fan deltas near the termini of glaciers. Velocities of landslides and subsidence were consistent (within our measurement error) and all <6 cm/yr. Continued remote monitoring of these areas is warranted. These InSAR results could be combined with other data sets to produce a predictive subaerial landslide susceptibility map for the Glacier Bay region. Given ongoing warming from climate change, systematic regional InSAR surveys of degrading cryospheric terrain, such as the one described here, could likely become an essential hazard and risk assessment tool in the near future.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.