Introduction

The effects of climate change have the potential to impact slope stability (e.g., Sidle and Burt 2009; Crozier 2010; Gariano and Guzzetti 2016). Adverse impacts are expected to be greatest at high altitudes and northerly latitudes where degradation of permafrost in rock and soil, debuttressing of slopes as a result of glacial retreat, and changes in ocean ice-cover (and therefore less protection for coastal bluffs from waves and currents) are likely to increase the susceptibility of slopes to landslides (e.g., Haeberli and Beniston 1998; Gruber and Haeberli 2007; Coe and Godt 2012; Huggel et al. 2012; Geertsema et al. 2013; Deline et al. 2015a; Cloutier et al. 2016). In the USA, the greatest increases in air temperature and precipitation from climate change (e.g., van Oldenborgh et al. 2013; Walsh et al. 2014; Wuebbles et al. 2014), and thus the most dramatic adverse impacts to slope stability, are expected in Alaska (Coe 2016). But, to date, conclusive evidence of changes in the frequency, magnitude, or mobility of landslides in Alaska has been elusive (e.g., Huggel et al. 2012). This situation is due, at least in part, to a paucity of systematic records on historical landslides that can be compared to data from new landslides, which is a problem in alpine and glacial environments in many parts of the world (e.g., Schneider et al. 2011; Wood et al. 2015).

Glacier Bay National Park and Preserve in southeast Alaska (Fig. 1) has attracted the attention of landslide researchers in the last few years because of several large rock avalanches (e.g., Mt. La Perouse, Fig. 2) that were initially detected seismically (e.g., Geertsema 2012; Ekström and Stark 2013; Petley 2014, 2016a). Following these rock avalanches, we asked ourselves a simple question: How do these recent events compare to older events in the historical record? To answer this question, we used Landsat imagery to compile an inventory of rock avalanches in Glacier Bay for the period from 1984 through 2016 (Bessette-Kirton and Coe 2016). Landsat imagery provided a consistent medium that allowed us to systematically detect and map avalanches that covered areas larger than about 10 ha (100,000 m2, 0.1 km2). Although the 30-m Landsat pixel resolution is not ideal for the identification of small landslides, it is adequate to discern the large rock avalanches typical of Alaska (e.g., Post 1967; Jibson et al. 2006). Higher-resolution imagery of the area is available (e.g., 15-m Landsat 7 and 8 panchromatic band imagery, 0.5-m DigitalGlobe imagery), but, to avoid biasing our inventory toward more rock avalanches being mapped from high-resolution imagery, we used 30-m Landsat imagery for the entire 33-year period.

Fig. 1
figure 1

Map of the Glacier Bay study area with rock avalanches listed in Table 1

Fig. 2
figure 2

February 16, 2014 rock avalanche from Mount La Perouse (avalanche #3 in Table 1). Relief from head scarp to toe is about 1770 m. Photo taken on March 6, 2014 by Marten Geertsema

In this paper, we analyze and assess our rock-avalanche inventory (Bessette-Kirton and Coe 2016) for changes in frequency, magnitude, and mobility over the 33-year period of record. Because the Glacier Bay area is tectonically active (e.g., Elliott et al. 2010), and earthquakes can trigger rock avalanches (e.g., McSaveney 1978; Jibson et al. 2006), part of our work focused on searching available historical earthquake catalogs to determine if any of the rock avalanches in our inventory were triggered by earthquakes. The discussion section of this paper presents some possible reasons for the pattern of observed changes. This paper expands on a preliminary analysis of the rock-avalanche inventory contained in Coe et al. (2017).

Geologic setting

The study area is a glaciated peninsula located along the transform boundary of the Pacific and North American plates. The plate boundary is marked by the strike-slip Fairweather fault (Pflaker et al. 1978, Fig. 1) which accommodates a large portion of the relative motion between the plates through dextral slip of about ~ 43 mm/year (Elliott et al. 2010). Elevations in the area range from sea level up to 4670 m at Mount Fairweather at the US/Canadian border (Fig. 1). Much of the study area is covered by ice, either in the form of valley glaciers or from the Brady ice field. In the last 250 years, the area has been undergoing rapid uplift (up to 30 mm/year) caused primarily by viscoelastic rebound from post-Little Ice Age (LIA) deglaciation (Larsen et al. 2005). Glacier Bay proper (a fjord, Fig. 1) is a result of post-LIA deglaciation (e.g., Connor et al. 2009).

Bedrock in the study area is exposed at mountain peaks, flanks, and ridge lines and in coastal zones around the periphery of the study area. The bedrock geology is complex and consists of Paleozoic and Mesozoic accretionary terranes containing Tertiary sedimentary, intrusive, and volcanic rocks (Brew et al. 1978; Wilson et al. 2015). Bedrock units that crop out in the ice-covered portion of the study area are as follows: the Triassic to Cretaceous Kelp Bay Group consisting of phyllite, quartzite, greenschist, greenstone, graywacke, and graywacke semischist; the Cretaceous Chugach accretionary complex consisting of flysch, graywacke, and volcanic rocks; and unnamed Tertiary granites, gabbros, and diorites (Wilson et al. 2015). Gruber (2012a, b) indicates a high likelihood that permafrost is present in high-elevation areas in the study area (see Fig. 3 and Gruber (2012b) for map data).

Fig. 3
figure 3

Permafrost zonation map for the Glacier Bay study area (modeled zonation from Gruber (2012b). Head scarps of the three largest rock avalanches in the period 1984–2016 are labeled

Rock slides and rock avalanches in Alaska have long been of interest to scientists and land managers because of their ability to cause glaciers to advance (e.g., Field 1968; McSaveney 1975) and generate hazardous tsunamis (e.g., Petley 2016b). In the Glacier Bay region, a large, earthquake-triggered rock slide in 1958 entered Lituya Bay on the west side of the study area (Fig. 1) and produced a tsunami that killed two people (Miller 1960). At Tidal Inlet, just outside the study area to the east (Fig. 1), a large active rock slide has the potential to generate waves greater than 10 m in the western arm of Glacier Bay near the mouth of Tidal Inlet (Wieczorek et al. 2007). Within the last 5 years, documentation of rock avalanches has begun because of potential links to climate change (e.g., Geertsema 2012; Geertsema et al. 2013).

Methods

For our inventory mapping, we used Landsat imagery available from the U.S. Geological Survey (USGS) EarthExplorer website (USGS 2017a). Specifically, we used short-wave infrared bands (band 5, Landsat 5 and 7, and band 6, Landsat 8) for mapping (Bessette-Kirton and Coe 2016) because short-wave infrared wavelengths penetrate thin clouds and discriminate changes in ground moisture content (see USGS 2017b). We found that the short-wave infrared bands were most effective for identifying rock avalanches on glaciated terrain in the oceanic climate of southeast Alaska. Individual short-wave infrared band images and multiband false-color images were used to identify and delineate rock avalanches. Band combinations 453 (Landsat 5/7) and 564 (Landsat 8) were found to be most effective for differentiating between snow/ice and bare ground.

Landsat images acquired from late spring (May) to early fall (October), with minimal cloud cover, were most effective for mapping. To initially identify rock avalanches, we examined a minimum of one image acquired from May to October for each calendar year from September, 1984 to September, 2016 (Bessette-Kirton and Coe 2016). We then acquired subsequent imagery to refine the area covered and the range in possible dates of occurrence for each of the initially identified rock avalanches (Table 1). We were able to determine the exact dates of three recent events from other sources (Table 1).

Table 1 Rock-avalanche data from Bessette-Kirton and Coe (2016), listed from oldest to youngest. Rock avalanche #5 is the Lituya Mountain avalanche (Geertsema 2012), avalanche #3 is the Mount La Perouse avalanche (Petley 2014), and avalanche # 20 is the Lamplugh Glacier avalanche (Petley 2016a; Bessette-Kirton 2017). All rock avalanches terminated on glaciers. Aspect for each rock avalanche was measured from a single 5-m DEM cell and may or may not be reflective of the overall rock-avalanche travel direction

We identified and mapped rock avalanches (e.g., Fig. 4) by locating areas with (1) high contrast compared to surrounding, undisturbed snow and ice; (2) different spectral signatures between successive Landsat images; and (3) lobate forms typical of rock-avalanche deposits. For each mapped rock avalanche, we used a 5-m digital elevation model (DEM) from the USGS National Elevation Dataset (USGS 2017c) to record the maximum travel distance (length, L), avalanche relief (height, H, the change in elevation between start (head scarp) and end (toe) points measured along a curvilinear centerline), and slope direction (aspect) of head scarps.

Fig. 4
figure 4

False-color composite Landsat images showing a the smallest (0.11 km2) rock avalanche mapped (Landsat bands 453, avalanche 9, Table 1) and b the largest (22.19 km2) rock avalanche mapped (Landsat bands 564, avalanche 20, Table 1, the June 28, 2016 rock avalanche onto the Lamplugh Glacier). Each rock avalanche is outlined with a white polygon. Note the difference in scale between a and b

Our inventory may not record all of the rock avalanches that have occurred since 1984. We make this statement because it is possible that snow and/or cloud cover prevented us from identifying and mapping some small (~< 0.5 km2) rock avalanches, especially those that might have occurred in winter months at elevations above the snowline. We are reasonably confident that we did not miss rock avalanches larger than ~ 0.5 km2 because deposits from those rock avalanches are visible in Landsat images for multiple years after they occur. Additionally, in several cases (as noted in Bessette-Kirton and Coe 2016), rock-avalanche head scarps were not clearly visible in the imagery because of shadowing, or a lack of contrast between ground disturbed by the rock avalanche, and undisturbed snow/ice or bare rock. In these cases, the head scarp location was approximated based on the source area topography in the 5-m DEM.

As estimates of magnitude and mobility are needed for hazard and risk assessments, we used our inventory, and the methods described in Coe et al. (2004), to estimate the annual exceedance probabilities for specific magnitude and mobility values. These estimates provide initial magnitude/frequency and mobility/frequency relations for the study area. However, because our inventory is likely a minimum record, our estimated recurrence intervals are probably too large and our annual exceedance probabilities are probably too small.

Because we are interested in the possible effects that climate change has on rock-avalanche characteristics, we assessed whether or not any of the rock avalanches in our inventory could have been triggered by earthquakes. To do this assessment, we used our head scarp locations (Table 1) to search all 27 catalogs of earthquakes available through the USGS earthquake hazards website (USGS 2017d). To conduct our search, we used guidance provided in Keefer (1984) and Jibson (2013). Specifically, for the range in dates of occurrence for each rock avalanche (Bessette-Kirton and Coe 2016 and Table 1), we searched for earthquakes that were within the magnitude and maximum distance criteria provided by Keefer (1984) and representing earthquakes below the “disrupted landslides” maximum threshold line in Fig. 5. Disrupted landslides have experienced movement in the past and may, therefore, have established failure surfaces. We were conservative in our use of the disrupted landslides threshold because these types of landslides are most sensitive to triggering by earthquakes (for example, at a distance of 100 km from an earthquake epicenter, a disrupted landslide could be triggered by a smaller magnitude earthquake than the one required to trigger a coherent landslide).

Fig. 5
figure 5

Variation in earthquake epicentral distance from landslides for earthquakes of varying magnitudes. Magnitude is moment magnitude for M ≥ 7.5 and surface magnitude for M < 7.5. For each rock-avalanche head scarp location and date range given in Table 1, we searched for possible earthquakes with distances and magnitudes in the shaded area. Modified after Keefer (1984) and Jibson (2013)

Results

We examined 104 Landsat images and mapped 24 rock avalanches (Bessette-Kirton and Coe 2016, Table 1, Figs. 1 and 3). We found no earthquakes that met the magnitude and distance criteria shown in Fig. 5. In other words, we found no earthquakes likely to have triggered any of the mapped rock avalanches. All of the mapped rock avalanches were in the northern half of the study area and originated in probable rock-permafrost zones (Fig. 3) at elevations between 1134 and 3323 m above sea level (asl, Fig. 6a). Unlike some other recent rock avalanches in Alaska which initiated from mountain slopes near the termini of valley glaciers (e.g., the Icy Bay rock avalanche on October 17, 2015; Petley 2016b; Rozell 2016), all rock avalanches in our inventory initiated from high-mountain ridges, flanks, and peaks well away from valley glacier termini. Seventy-five percent of rock-avalanche source areas had slope directions (aspects) between 340° and 90° (i.e., north- to northeast-facing slopes, Fig. 6b).

Fig. 6
figure 6

Variation in a elevations and b aspects for rock-avalanche head scarps sampled from a 5-m DEM. Aspects are shown on a 360° Rose diagram. Numbers on circles are percentages of total

Further analysis of avalanche characteristics indicated that the mean-recurrence interval between rock avalanches for the 33-year period of record was 505 days (Fig. 7a), but that there were three distinct clusters of rock-avalanche activity where the recurrence interval for successive events was less than the mean (Fig. 7a). These three clusters are between ~September 7, 1984 (the beginning of the period of record), and July 20, 1986; June 15, 1994, and June 27, 1995; and June 11, 2012, to present (i.e., the last cluster had not ended as of November 2016). The first two clusters (clusters 1 and 2, Fig. 7b) have five rock avalanches each and no exceptionally large (> 5 km2) avalanches. The third cluster (cluster 3, Fig. 7a) contains 10 rock avalanches and the largest rock avalanches in the period of record including the rock avalanches on June 11, 2012 at Lituya Mountain; on February 16, 2014 at Mount La Perouse; and on June 28, 2016 at Lamplugh Glacier (Table 1, Fig. 7b). The sizes of these three recent rock avalanches exceed the next largest avalanche in the period of record (avalanche 17, Table 1) by factors of 1.5, 2.1, and 5.9, respectively (Table 1, Fig. 7b). Additionally, the head scarp elevations for two (Lituya and La Perouse) of these three avalanches were the highest in the period of record (Fig. 6). The large size of rock avalanches in the third cluster corresponds with increased rock-avalanche mobility, with the third cluster having the lowest mean H/L value of the three clusters (Fig. 7c). The correlation between magnitude (in this case, area, rather than volume) and mobility has long been recognized by landslide researchers and exists for a full range of landslide magnitudes (e.g., Corominas 1996).

Fig. 7
figure 7

Diagram showing rock-avalanche characteristics and annual probabilities. Median dates were used for avalanches with a range of possible dates of occurrence (see Table 1). The Poisson probability model was used to estimate annual probabilities using the methods described in Coe et al. (2004). a Recurrence intervals for the 33-year period of record. Recurrence intervals were plotted at the median date of the second rock avalanche in the sequence. For example, a recurrence interval calculated from rock avalanches with median dates of January 1, 2012, and January 1, 2013 was plotted at the January 1, 2013 date. b Magnitude (size, expressed as total area including the avalanche scar, travel path, and deposit) for the 33-year period of record. c Height/length for the 33-year period of record. d Magnitude/frequency relations for 33- and 28-year periods of record. e Mobility/frequency relations for 33- and 28-year periods of record

To illustrate the impact that the last cluster (cluster 3, ongoing) of large, highly mobile rock avalanches has had on hazard estimates in the study area, we established two sets of magnitude/frequency and mobility/frequency relations. The first set (Fig. 7d, e) is based on all of the data in the 33-year period, and the second set (Fig. 7d, e) is based only on data from the 28-year period prior to the start of the last cluster (i.e., June 11, 2012). The difference in the magnitude/frequency relations for the two periods of record indicates that annual probability of exceeding various avalanche sizes has increased between 3 and 14% (Fig. 7d). For example, the annual probability of having an avalanche larger than 1 km2 increased from 13% for the 28-year period of record to 24% for the 33-year period of record. The differences in mobility/frequency relations for the two periods indicate that the annual probability of having mobility values from 0.2 to 0.7 (lower H/L values indicate higher mobility) has increased between 15 and 17% (Fig. 7e). These estimates of probability may be overly influenced by clustering that distorts the record because the record length is relatively short. Continued systematic inventory mapping in the study area is needed to fully understand the significance of clustering. However, with the data that we currently have available, rock avalanches in the third cluster are clearly unusual because of their large size and low H/L values.

Comparison to other historical rock avalanches on glaciers

Direct comparisons of our data to other data sets from glaciated regions of Alaska and western Canada are somewhat difficult because studies often do not include consistent geometric and morphometric information, and most previous studies are for earthquake-triggered rock avalanches. Ideally, we would compare our data with a large amount of pre-1984 data from within our study area. This approach would allow us to evaluate rock-avalanche characteristics over a time span longer than 33 years, while holding other variables, such as uplift, tectonic deformation, structure, and lithology in source areas, constant. However, from existing literature, we have only been able to identify four pre-1984, rock avalanches within the study area (Table 2). The largest rock avalanche of the four is the 1965(?) Fairweather-glacier avalanche with an area of 8.5 km2 and an H/L value of 0.33. Compared to the 10 rock avalanches in our most recent cluster (cluster 3), the area of the Fairweather rock avalanche ranks as the second highest, but the H/L is tied for eighth (i.e., 7 of the 10 more recent avalanches are more mobile, see Table 1). The 2016 Lamplugh Glacier rock avalanche has an area that is about 2.6× larger (22.19 km2) and an H/L value less than half (0.15) of the Fairweather rock avalanche.

Table 2 Data for other historical rock avalanches within the study area. We presume that these rock avalanches were not triggered by earthquakes because the original sources do not identify earthquake triggers. No data are available for cells that are blank

The mobility and area of the Lamplugh Glacier rock avalanche are also exceptional when compared to older data sets from outside of the study area. For example, the H/L value for the Lamplugh Glacier rock avalanche is smaller than any of the 16 H/L values listed by Evans and Clague (1999) and Geertsema and Cruden (2014) for rock avalanches on glaciers in western British Columbia and the larger Canadian Cordillera, respectively. The area of the Lamplugh Glacier rock avalanche is at least 2× greater than any of the ~ 75 rock avalanches triggered by the M 9.2 Anchorage, Alaska earthquake of March 27, 1964 (Post 1967), and any of the ~ 7 rock avalanches triggered by the M 7.9 Denali fault, Alaska earthquake of November 3, 2002 (Jibson et al. 2006). The largest of the avalanches triggered by the Denali earthquake (the McGinnis Glacier North rock avalanche, Jibson et al. 2006) had an area of 10.61 km2 and H/L equal to the Lamplugh rock avalanche (0.15). The rest of the Denali avalanches all had areas less than 7 km2 and H/L values of 0.16 to 0.23. Also, the aspects of source areas for the rock avalanches triggered by the Anchorage and Denali earthquakes (Fig. 8) are much more broadly distributed compared to the dominant north to northeast directions of rock-avalanche source areas in our Glacier Bay data set (Fig. 6).

Fig. 8
figure 8

Rose diagram showing aspect data for source areas of 113 rock avalanches triggered by the M 9.2 Anchorage and the M 7.9 Denali earthquakes. Numbers on circles are percentages of total. All aspects were measured from maps presented in Post (1967) and Jibson et al. (2006). Compare with our Glacier Bay dataset in Fig. 6

In comparison to the worldwide data set of 56 rock avalanches compiled by Deline et al. (2015b), the Lamplugh Glacier rock avalanche ranks first in area and tied for the fourth smallest H/L. This statement is qualified by the fact that the Deline et al. (2015b) compilation includes some of the rock avalanches from publications described in the previous two paragraphs, includes 21 rock avalanches without documented areas, and gives area as the surface area of the deposits only. Our area for the Lamplugh Glacier rock avalanche (22.19 km2) includes both the source (~ 1 km2) and deposit (~ 21 km2) areas. The 1970 Nevados Huascarán rock avalanche in Peru had a similar area (22.00 km2) to the Lamplugh Glacier rock avalanche. The next largest rock avalanche listed by Deline et al. (2015b) is the 10.21-km2 (deposit only) McGinnis Glacier North rock avalanche of Jibson et al. (2006). In terms of H/L values, the 1967 Steinholtsjokull Glacier rock avalanche in Iceland and the 2002 Kolka-Karmadon rock avalanche in Russia had H/L values smaller than the Lamplugh Glacier rock avalanche.

Discussion

Our results indicate that the probability of large (> 5 km2), highly mobile (H/L < 0.3) rock avalanches has increased in the Glacier Bay study area since June, 2012. The question is: why did large rock avalanches become more common? An additional, related question is: what factor or factors controlled the timing of rock-avalanche clusters? Because we have excluded earthquakes as the cause of any rock avalanches in our inventory, an exploration of other factors that could induce changes in the effective stress in rock outcrops in glaciated environments is needed. Four possible factors include debuttressing of rock slopes from thinning of glaciers, increased precipitation, changes in regional crustal strain, and degradation of rock permafrost.

A recent analysis of glacial thinning in Glacier Bay National Park indicates that thinning between 1995 and 2011 was most pronounced in glacial ablation zones and near glacier termini (Loso et al. 2014) at elevations lower than our rock-avalanche source areas. Additionally, a longer-term study of glacier surface elevation changes from 1948 to 2000 showed both increases and decreases in elevation (from ~− 2 to ~+ 2 m) in accumulation zones near rock-avalanche source areas (Larsen et al., 2007), which suggests to us that changes in snowfall and glacier thicknesses were probably not factors that influenced the timing or size of rock avalanches. Recent work on the Brady Glacier in the southeastern part of the study area indicated that average glacier equilibrium-line altitude measured from 2003 to 2011 (745 m) was about 145 m higher than for measurement periods from 1961 to 1962 and 1976 to 1984 (~ 600 m, Pelto et al. 2013), which suggests increasing air temperatures at progressively higher elevations.

Variations in air temperatures and precipitation in Alaska can be influenced by the El Niño Southern Oscillation (ENSO) and the Pacific Decadal Oscillation (PDO) ocean-atmosphere phenomena. The ENSO causes climatic variations in the Americas on seasonal to annual time scales. The PDO is an ENSO-like phenomenon that can persist for decades and preferentially impacts climatic variations in the North Pacific Ocean and North America (Mantua et al. 1997). Two indices have been developed to provide a numerical measure of the warm and cool phases of ENSO and PDO, the Multivariate ENSO Index (MEI, Wolter and Timlin 2011) and the PDO index (Mantua et al. 1997). To investigate the possible influence of the ENSO and PDO phenomena on the timing of our rock avalanches, we compared data from both indices with our rock-avalanche clusters (Fig. 9). A visual inspection of Fig. 9 indicates that there is not a consistent pattern between either of the phenomena and the timing of rock-avalanche clusters. For example, the first cluster occurred during a cool ENSO phase and a warm PDO phase, whereas the second cluster occurred during a warm ENSO phase and a cool PDO phase. The third cluster occurred during both cool and warm phases of both phenomena.

Fig. 9
figure 9

Diagram showing the timing of rock-avalanche clusters with respect to the El Niño Southern Oscillation (Multivariate ENSO Index, MEI) and the Pacific Decadal Oscillation (PDO) index data from 1984 to 2017. Question marks signify that the start and end times for the first and last clusters (respectively) are not known

In regard to the possible influence of changing regional crustal strain on rock-avalanche occurrence, McSaveney et al. (2015) hypothesized that an observed increase in rock-avalanche frequency in the Southern Alps of New Zealand correlates with elastic strain that has been accumulating in the crust since the last great earthquake. The hypothesis is that relatively loose, fractured rock near the surface is breaking off (as rock avalanches) more frequently because the fit with underlying rock is becoming progressively worse as elastic strain accumulates in the crust. In the Glacier Bay region, with ongoing post-LIA uplift and regional-tectonic strain, elastic strain is undoubtedly accumulating. However, two characteristics of our rock avalanches do not fit well with accumulating elastic strain being the primary control on rock-avalanche occurrence. Specifically, the north to northeast slope aspects of 75% of rock-avalanche source areas, and the predominant spring and summer dates of occurrence (see Table 1), do not fit well with the accumulating strain hypothesis of McSaveney et al. (2015). We would expect both of these characteristics to be more broadly distributed if rock-avalanche occurrence was primarily controlled by accumulating strain. Regarding the apparent strong seasonal control on rock-avalanche occurrence, we acknowledge that this observation may be biased because we may be missing small rock avalanches that occurred in the fall and winter seasons, and we may have misdated a few large avalanches that occurred in the fall and winter because they became snow-covered and were not identified until spring and summer when the snow melted off them. Accumulating elastic strain is likely a conditioning factor that weakens rock masses in the Glacier Bay region, but a different factor probably controls the timing and size of rock avalanches.

We hypothesize that degradation of rock permafrost is the primary factor that controlled the timing and size of rock avalanches in the Glacier Bay region. Existing literature on this topic highlights permafrost degradation caused by warming air temperatures as the most likely cause of changes in rock-avalanche and rock-fall activity in many areas of the world (e.g., Schär et al. 2004; Geertsema et al. 2006; Harris et al. 2009; Deline et al. 2011; Ravanel and Deline 2011; Fischer et al. 2013), but attributing specific events to permafrost degradation is difficult, especially without monitoring data at the precise locations of rock-slope failures (Allen et al. 2011; Huggel et al. 2012). Even so, recent work on the degradation of soil permafrost in Alaska suggests that the areal extent of near-surface permafrost will decline by 16 to 24% by the end of the twenty-first century (Jafarov et al. 2012; Pastick et al. 2015).

Gaining a better understanding of the processes by which permafrost in steep-rock slopes degrades and leads to rock failures is a topic of current research (e.g., Hasler et al. 2012; Krautblatter et al. 2012). Existing work indicates that subsurface conduction and advection of heat can degrade permafrost in rock (e.g., Gruber and Haeberli 2007; Noetzli and Gruber 2009). In particular, advection through the movement of water in fractures can create pathways where deep (tens of meters) thawing can occur along fractures leading to an increase in effective stress, a reduction in rock-mass strength, and destabilization of large volumes of rock (e.g., Gruber and Haeberli 2007; Fischer et al. 2010; Hasler et al. 2011). Within this context, our observation that rock avalanches preferentially initiated from north- to northeast-facing slopes (Fig. 6b) was unexpected. Work by Gruber et al. (2004) suggests that slopes with these aspects likely have a larger extent of rock permafrost (compared to south-facing slopes) and may be preferentially susceptible to failure during warm months of the year. Gruber et al. (2004) found that north-facing slopes in the European Alps experienced a greater number of rock failures than south-facing slopes during the summer of the 2003 heat wave. They attributed this observation to south-facing slopes regularly experiencing larger subsurface temperature fluctuations than north-facing slopes. Thus, the anomalously warm air temperatures of 2003 did not impact them as much as north-facing slopes, which had more extensive rock permafrost and were not calibrated to anomalously warm air temperatures. Indeed, in our study area, an examination of the dates of occurrence in Table 1 suggests that rock avalanches may predominantly occur during the spring and summer months.

There are no monitoring sites where rock temperatures and hydrologic conditions are being recorded in our study area. We therefore have to extrapolate from available air temperature data to make interpretations. Permafrost degradation is highly sensitive to air temperature (e.g., Deline et al. 2015a), and there has been some previous success in identifying air temperature patterns prior to large rock avalanches in Alaska and the European Alps (Huggel et al. 2010). Huggel et al. (2010) found that there were unusually warm temperatures (well above freezing and up to 8.5 °C above long-term average temperatures) in the weeks or days before avalanches, and that prior to some of the slope failures, there were sudden drops in temperature (usually below freezing) after warm periods in the hours to days before the avalanches.

In the vicinity of our study area, there are few meteorological stations, and those that do exist tend to have relatively short periods of record and are located at elevations at or near sea level. In contrast to the elevations of these stations, the mean elevation of the head scarps of the rock avalanches in our inventory (see data in Bessette-Kirton and Coe 2016) is 1968 m asl. Given this situation, we analyzed meteorological data from two National Atmospheric and Oceanic Administration (NOAA) stations that have relatively long periods of records and relatively close proximity to our study area. The two stations were the Glacier Bay station (station ID: COOP 503294) near Gustavus (Fig. 1) and the Juneau International Airport station (station ID: COOP 504100), respectively, located about 75 and 150 km southeast of the center of the study area. The Glacier Bay station is at an elevation of 4.2 m asl and has a record from 1966 to 2016. But the data are incomplete with about 35 years missing from the annual temperature summary records provided by the National Climatic Data Center (NCDC). The Juneau station is at an elevation of 1.5 m and has records from 1942 to 2016. Data from the Juneau station are more complete than Glacier Bay, with only 8 years missing from annual summary data.

Data from the Juneau station show an increase in mean annual air temperature (MAAT) of about 2 °C during the period of record (Fig. 10a). Projecting the Juneau data to an elevation of 1968 m asl using a moist adiabatic lapse rate of − 0.005 °C/m (e.g., Hinzman et al. 1998) suggests that the present-day MAAT at 1968 m in the Juneau area should be well below freezing (0 °C), which is consistent with a general lack of rock avalanches being reported at high elevations in the mountainous, continental interior area near Juneau. The same pattern is observed with mean annual maximum air temperature (MAMT) at Juneau (Fig. 10b), although the long-term trend of MAMT at 1968 m asl is approaching 0 °C and, as of 2015, was within 1 °C of reaching 0 °C (Fig. 10b).

Fig. 10
figure 10

Annual air temperature data from the Juneau International Airport and Glacier Bay stations. Data at head scarp elevations were estimated using a moist adiabatic lapse rate of − 0.005 °C per m. Best-fit trendlines are shown. a Mean annual air temperature at Juneau. b Mean annual maximum air temperature at Juneau. c Mean annual air temperature at Glacier Bay. d Mean annual maximum air temperature at Glacier Bay. Trendline suggests that the mean annual maximum temperature exceeded 0 °C in about 2010

At the Glacier Bay station, the climatic pattern is similar to Juneau (Fig. 10c, d), with one important exception. At Glacier Bay, the MAMT trendline for data projected to 1968 m asl exceeded 0 °C in about 2010. The exceptionally large rock avalanches began in 2012 with Lituya Mountain (Table 1, Fig. 7, the third cluster of activity in Fig. 7a). Of course, the sparseness and projection (to an elevation of 1968 m) of the data from the Glacier Bay station make the position and slope of the trendline, and the timing of the 0 °C crossover, inexact and conjectural. Even so, existing literature (e.g., Davies et al. 2001; Krautblatter et al. 2013) indicates that strength of ice in rock can decrease with warming between rock temperatures of − 5 and 0 °C. Therefore, the established warming trends in this area of Alaska are likely more important than the position and timing of maximum air temperatures rising above 0 °C.

Some additional climatological data for the study area are available from the North American Freezing Level Tracker (NAFLT), which was developed by the Western Regional Climate Center (WRCC) and provides historical data on the elevation of the freezing level (0 °C) in North America (see data at WRCC 2017). NAFLT identifies the elevation at which 0 °C is first encountered when progressing downwards from the upper atmosphere. Air temperatures are taken from mean daily temperature profiles derived from four, 6-h temperature averages from a National Centers for Environmental Protection (NCEP)/National Center for Atmospheric Research (NCAR) Global Reanalysis product (Kistler et al. 2001) with a grid spacing of 2.5° and data availability from 1948 to present. Freezing elevations at specific geographic locations are interpolated from the 2.5° grid. Average elevations of freezing levels for varying lengths of time can be computed and compared to a long-term mean value.

For the Glacier Bay study area, we evaluated a variety of different time periods ranging from 1 to 24 months. However, we found mid-winter (January) and mid-summer (July) months most useful for comparison with the timing of our rock-avalanche clusters (Fig. 11) because these months contained the lowest and highest freezing elevations for individual years. A visual examination of Fig. 11a reveals a correlation between January freezing elevations and rock-avalanche clusters. All three clusters occurred during periods with two or more consecutive years of January freezing elevations above the long-term average. In other words, the rock-avalanche clusters occurred during multiyear periods of time when mid-winter freezing temperatures were at elevations higher than normal. This same pattern holds for the two historical rock avalanches that occurred in 1961 (Fig. 11 and Table 2). There is not a similar systematic pattern for July freezing elevations (Fig. 11b).

Fig. 11
figure 11

Diagram showing the timing of rock-avalanche clusters with respect to data from the North American Freezing Level Tracker (NAFLT) for the approximate center of the northern half of the study area (latitude 58.74° N, longitude 137.27° W). Mean elevations of 0 °C (gray horizontal lines) are from the period 1981–2010. Question marks signify that the start and end times for the first and last clusters (respectively) are not known. a Elevations of 0 °C for the month of January. b Elevations of 0 °C for the month of July

Figure 11 also reveals that, beginning in 1980, there are a greater number of years (21 of 37, 57%) with January elevations above average compared to years with January elevations below average (16 of 37, 43%). This pattern does not exist for post-1980 July elevations where the number of years with freezing elevations above and below average is nearly evenly split (18 above average, 19 below average). Additionally, the freezing elevation for January 2014 (the approximate center of the third cluster) is the highest on record. Also, the only two periods with five consecutive years of above average freezing elevations correspond with the first (1984–1986) and the last (2012–2016) rock-avalanche clusters. These results suggest to us that the timing of rock-avalanche clusters, and the increase in rock-avalanche magnitude and mobility, is determined more by winter temperatures than by summer temperatures. Presumably, warmer than normal winter temperatures weaken rock masses that then fail when temperatures approach or exceed 0 °C in the late spring and summer.

These results fit with existing literature that indicates that trends in air temperature increases in Alaska are more pronounced in winter (Arendt et al. 2009; Wendler and Shulski 2009). In regard to the third rock-avalanche cluster from 2012 to 2016, the US 2016 annual climate report (NOAA 2017) states that each of the 3 years from 2014 to 2016 were the three warmest on record for the State of Alaska, with 2016 being the warmest, 2014 the second warmest, and 2015 tied for third place (with 2002). The average winter temperature for Alaska for the December 2015 to February 2016 period was 10.6 °F (5.94 °C) above average and set a new record (NOAA 2017). Walsh et al. (2017) note that the exceptional warmth in the 2015/2016 winter was most pronounced for the January through April time period. The record-breaking, multiyear and winter-season warmth is the most likely explanation for increase in rock-avalanche size in the third cluster. Given that these recent large rock avalanches likely have deep (> 10 m) failure surfaces (e.g., see volume estimates for the Lituya Mountain and Lamplugh rock avalanches in Geertsema (2012) and Bessette-Kirton (2017), respectively), prolonged warming may be needed for the transmission of heat and/or water to the depth of the failure surfaces. This statement is supported by Krautblatter et al. (2013) who proposed a conceptual model based on experimental results which suggests that rock-mass strength decreases as rock-permafrost slopes warm and that these changes can occur over time frames from months to millennia. Regardless of our imprecise understanding of the exact processes responsible for rock-slope failures due to warming temperatures, the observed increase in large, highly mobile rock avalanches in the Glacier Bay study area indicates that a critical temperature threshold was most likely exceeded, at least for the third cluster from 2012 to 2016.

Our interpretation is that the third rock-avalanche cluster was primarily caused by rock-permafrost degradation from the cumulative effect of long-term warming and the prolonged, record-breaking warmth of 2014–2016. We recognize and acknowledge, however, that accumulating elastic strain, glacial thinning, and increased precipitation could be contributing factors that conditioned the slopes for failure. An application of the principle of Occam’s razor, which asserts that the simplest explanation for a phenomenon is more likely to be accurate than more complicated explanations, suggests that these rising and record-breaking temperatures should be considered the primary driving force for the observed increases in rock-avalanche size and mobility in the Glacier Bay region. Even so, continued monitoring of rock-avalanche activity in the study area should be done to critically assess the long-term reliability of the interpretations and conclusions from this study.

Conclusions

Our study of 24 rock avalanches which occurred in the western part of Glacier Bay National Park and Preserve between 1984 and 2016 indicates that the probability of occurrence of large (areas > 5 km2), highly mobile (H/L < 0.3) rock avalanches has increased by at least a factor of 2 in the last 5 years (2012–2016). These large rock avalanches began to occur approximately 2 years after the long-term trend in mean annual maximum air temperature is thought to have exceeded 0 °C. Possibly more important, however, most of these large rock avalanches occurred during a multiple-year period of record-breaking winter and spring air temperatures. We suggest that large rock avalanches in the study area may have been caused by rock-permafrost degradation. However, we recognize that other factors, such as accumulating elastic strain, glacial thinning, and increased precipitation, may also play an important role in preconditioning the slopes for failure during warm temperatures.

The observations and results presented here provide baseline information on rock-avalanche characteristics that can be compared to future avalanche characteristics. Rock-avalanche activity in the study area should continue to be monitored to critically assess the long-term reliability of the interpretations and conclusions from this study. We acknowledge that the mapped rock-avalanche polygons that form the basis of this paper may be different, and less accurate, than polygons that could be mapped from higher-resolution satellite, aerial, and handheld imagery. We specifically chose not to use the highest-resolution imagery because we desired a long-term inventory that was unbiased by increasing image resolution. New mapping could create an inventory that fully utilizes recently available high-resolution imagery.