Periglacial debris-flow initiation and susceptibility and glacier recession from imagery, airborne LiDAR, and ground-based mapping

The rapidly moving slurries of mud, rocks, water, and wood known as debris flows present hazards in mountainous areas around the world. The nature and severity of those hazards may be changing as the variables affecting the process itself, such as climate and glacier extent, and the human interaction with that process both change. Climate change may bring an intensification of storms, and increasing temperatures are responsible for the retreat of alpine glaciers (Oerlemans, 1994). Expansion of human settlements and other infrastructure into mountainous areas may increase communities’ exposure to potentially destructive debris flows. In the Cascade Range of the Pacific Northwest, proglacial areas on composite volcanoes are typically steep and mantled with ample unconsolidated material and are therefore prone to debris flows (Blodgett et al., 1996; O'Connor et al., 2001; Osterkamp et al., 1986; Sobieszczyk et al., 2009; Walder and Driedger, 1994).

There is concern among managers and policymakers that glacier retreat, rising snowline elevations, and more frequent and more intense storms are leading to greater debris-flow hazards. On Washington's Mount Rainier in August 2001, rapid melting of Kautz Glacier directed meltwater into the adjacent watershed of Van Trump Creek, incised a new channel, and led to initiation of a series of debris flows, which formed a new debris fan in the Nisqually River near the main road into Mount Rainier National Park (Vallance et al., 2002; Vallance et al., 2003). In November 2006, an “atmospheric river” event (Neiman et al., 2008), i.e., a storm track carrying warm, moist air from the tropics, produced heavy rainfall at high elevations with little antecedent snowpack and triggered debris flows and flooding that caused major damage to infrastructure, particularly on Mount Hood in Oregon and Mount Rainier in Washington. In this study, we use the data from debris flows that occurred on Mount Rainier in 2001–2006 to improve our understanding of the relevant processes and as the basis for an attempt to map and model debris-flow initiation susceptibilities, a preliminary step toward quantifying hazards on that mountain.

The destructive potential of landslides and debris flows has motivated previous attempts to assess the hazards posed by these mass-movement processes. Hazards from large lahars (debris flows originating on volcanoes) induced by active volcanism have been mapped for Mount Rainier based on an empirical relationship between volume and inundation area (Griswold and Iverson, 2008; Iverson et al., 1998), but such mapping provides little information regarding smaller meltwater- and storm-induced lahars. Models based on steady-state hydrology and the effects of pore pressures on slope stability provide maps of relative susceptibility to shallow rapid landsliding and associated debris flows (Dietrich et al., 1995; Miller and Burnett, 2007, 2008; Montgomery and Dietrich, 1994), and models based on transient hydrology form the basis for forecasting debris flow hazards due to individual storms (Chleborad et al., 2008; Iverson, 2000), but the periglacial debris flows on Mount Rainier appear to be largely initiated by overland flow of water that leads to distributed gullying and associated bank slumping rather than by localized landslides (Copeland, 2009; Vallance et al., 2002). Despite recent observations of debris flows on Mount Rainier (Vallance et al., 2002) and other areas with debris flows initiated by overland flow (McCoy et al., 2010), we still lack a method for determining relative hazards in areas prone to such events.

While field-based monitoring is leading to greater understanding of debris flows triggered by intense rain and overland flow (McCoy et al., 2010), advances in remote sensing provide opportunities for addressing the relative susceptibility to overland flow-induced events at the landscape scale. Many hazard and susceptibility mapping efforts in the past were somewhat limited by the coarse resolution of readily available digital topography. Evolving technology for airborne laser swath mapping, or light detection and ranging (LiDAR), evolving methods for production of high-resolution “bare-earth” digital topography, and advancing efforts to collect these data present new opportunities for hazard and susceptibility mapping.

In our overarching research question, we ask whether recent glacier retreat in the Cascades has exposed areas with relatively high susceptibility to debris-flow initiation, and if so, do these additional areas compose a large fraction of high-susceptibility areas on the steep composite volcanoes of the Cascades. More specifically, this paper explores whether data from remote sensing and field mapping might provide a basis for debris-flow susceptibility mapping, which is necessary in order to answer the larger question regarding increased hazard and risk associated with glacier retreat. Here, we map debris-flow initiation sites with field- and remote sensing–based methods, determine limits of debris-free glacier ice and characterize topography with remote sensing–based methods, and demonstrate preliminary data-based modeling of areas susceptible to debris-flow initiation.

With its summit at 4393 m, Mount Rainier is the tallest volcano in the Cascade Range of Washington, Oregon, and California. The upper slopes of Mount Rainier are composed of unconsolidated Quaternary-age volcaniclastic and morainal material in steep edifices. The slopes also contain numerous glaciers and perennial bodies of snow and ice, as well as debris-covered stagnant ice masses of unknown spatial extent (C. Driedger, 2008, personal commun.). Mount Rainier is drained by five major glacier meltwater-fed braided rivers (Fig. 1), the headwaters of which are located within canyons, often 300–900 m deep, formed by lateral moraines from Pleistocene glaciers below adjacent ridges (Crandell, 1971). Stream gradients on the upper flanks of the mountain, above tree line, are 0.13–0.15 and, at the confluence with major rivers near park boundaries, 0.019–0.078 (Crandell, 1971). Tree line is located at elevations in the range 1600 to 2000 m; forests cover over 56% of Mount Rainier National Park, and the majority of the forests are over 350 yr old; soils are mainly from colluvium and tephra and have been described as podzolics and lithosols; wildfire is the major disturbance at Mount Rainier followed by snow avalanches and lahars; and humans have inhabited and modified the landscape of Mount Rainier National Park since the end of the last ice age (Dryden, 1968; Hemstrom and Franklin, 1982).

Mount Rainier and the entire Cascade Range formed as a by-product of the subduction of the Juan de Fuca plate under the North American plate. Mount Rainier formed over highly eroded middle Tertiary age rocks of the Ohanapecosh, Stevens Ridge, and Fifes Peaks Formation that were intruded by the Tatoosh pluton (Fiske et al., 1963). The volcano was built through progressive layers of thin lava flows, breccias formed by violent steam explosions, mudflows, pyroclastic flows, and infrequent pumice and ash eruptions (Fiske et al., 1963). The volcano first erupted about half a million years ago, other significant eruptions have occurred 2500 and 1000 yr ago, and minor steam eruptions were recorded in the 1840s (Sisson, 1995). The interbedded structure of composite volcanoes leads to inherent instability, as materials of different compositions, internal frictions, and permeabilities are juxtaposed. Hydrothermal activity further weakens the edifice. A sector collapse 5000 yr ago produced the Osceola Mudflow (Scott and Vallance, 1993).

Mount Rainier National Park is located only 70 km east of the Seattle-Tacoma metropolitan area, a proximity that makes Mount Rainier one of the most potentially deadly volcanoes in the United States due to the potential for a future volcanic eruption or initiation of a large lahar (debris flow initiating on a volcano). Large debris flows may be triggered by volcanic or seismic activity, and failure of hydrothermally altered rock (Reid et al., 2001) results in a clay content sufficient to retard dilution and drainage of the flow and cause deposition outside of park boundaries, so that lahar deposits from Mount Rainier cover much of the Puget Lowlands (Hoblitt et al., 1995). Smaller debris flows caused by hydrologic events do not characteristically exit the park's boundaries (Vallance et al., 2003).

The Cascade Range is influenced by a maritime climate characterized by 2.5 m of precipitation annually, making it one of the wettest places in the United States. The climate of Mount Rainier from October through May is cool and dominated by precipitation, but the modified maritime climate also results in warm, dry summer months. The northeastern and eastern sides of the mountain receive less precipitation because they are in the rain shadow of the dominant winter storm winds from the southwest and west (Hemstrom and Franklin, 1982).

Mount Rainier holds ∼25 named glaciers (Fig. 1), which comprise the largest ice volume of any peak in the conterminous United States, 4.4 billion m³ (Driedger and Kennard, 1984). The vast majority of glaciers have retreated significantly since the end of the Little Ice Age (Driedger, 1993; Nylen, 2004). Nylen (2004) reports a loss of 18.5% of total glacier cover from 1913 to 1971, with south- and north-facing glaciers diminishing 26.5% and 17.5%, respectively. A few glaciers with source areas at or near the summit have retreated minimally, and a minority, including the Carbon, Cowlitz, Emmons, and Nisqually Glaciers, advanced during the 1970s and 1980s as a result of precipitation patterns in the 1960s (Driedger, 1993; Nylen, 2004), but even these glaciers are currently at their historic minima since the end of the Little Ice Age (P. Kennard, 2008, personal commun.). Short-term glacier advance is masked by a much larger, long-term trend in glacier retreat. Since the 1980s all glaciers have thinned, retreated, or fragmented (Nylen, 2004). Nylen (2004) reports that winter snowfall at the Paradise Ranger Station meteorological station on Mount Rainier (Fig. 1) decreased 5.1 cm-yr⁻¹ during the period 1954–1994, while there were no significant trends in summer temperatures during the same period.

Initiation sites of debris flows occurring in 2001–2005 were identified from the air, in the field, and from aerial photographs by previous researchers (P. Kennard, 2008, personal commun.; Vallance et al., 2002). Channels that were impacted by debris flows in 2006 were identified through field reconnaissance, and all channels with suspected or identified 2006 debris flows were inspected in the field with reference to aerial and oblique photos. Identification of debris-flow deposits during field inspection relied on the presence of boulder levees and reverse grading (Pierson, 2005; Wilford et al., 2004). Debris-flow initiation sites were located by examining the upstream reaches of channels impacted by a 2006 debris flow for geomorphic signatures of debris-flow initiation. Data utilized in conjunction with field reconnaissance included aerial photographs, satellite imagery, and LiDAR (Table 1). The LiDAR data were used to produce a bare-earth surface elevation model and images of pulse return intensity with mean relative accuracy (the internal consistency of the LiDAR-derived locations) of 0.11 m, mean absolute accuracy (vertical accuracy) of 3.7 cm, and average ground point densities of 0.73 m⁻² in priority drainages and 0.87 m⁻² in all other areas (Watershed Sciences, 2009). All glaciers assessed in this study fall within the LiDAR acquisition period of September 2007. The 2006–2007 water year had lower than average snowfall and early snowmelt (Ellinger, 2010), so it is likely that laser returns from ice in 2007 represent glacier ice and not remnant snow. The LiDAR points were used to create a digital elevation model and pulse return intensity map with nominal resolutions of ∼1 m.

Initiation sites were identified by field reconnaissance as accurately as possible in steep, unstable terrain and were mapped using a combination of the 2006 National Agriculture Imagery Program (NAIP) orthophoto, 2007 aerial photos, and 2007–2008 LiDAR depending on the snow-covered area present in images and the chronology of debris flows in a particular channel. The 1-m data were used for detailed analyses around debris-flow initiation and other candidate sites; slopes derived from elevations and areas of ice and surface water from LiDAR intensities were used to find gullies (steep-walled and steep-headed channels downstream of glacier ice and upstream of the perennial channel network) and other streams. The 1-m LiDAR data were also used to calculate site characteristics, such as channel slope, aspect, channel geometry, distance to glaciers, glacier area, glacier volume, glacier debris–covered area, and the upslope contributing area, for the upstream-most gully head potentially involved with debris-flow initiation in each of the six drainage basins experiencing debris flows in 2001–2006. The aforementioned channel slopes were calculated over the first 9 m downstream of the gully head.

Distances from initiation sites to historic glacier margins were calculated from glacier coverages from Nylen (2004), who synthesized information from historic topographic maps, aerial photographs, and field surveys to reconstruct the glacier extents on Mount Rainier for the years 1913, 1971, and 1994. Nylen (2004) analyzed glacier terminus position, area, debris-covered area, and volume. Nylen's 1913 glacier coverage is based on field mapping by the U.S. Geological Survey (USGS) in 1910, 1911, and 1913 (Nylen, 2004; USGS, 1955). The 1913 map was revised using aerial photogrammetry to create the 1971 and 1994 maps. Nylen considered only active glacier ice and not stagnant ice downslope of the active ice terminus. Where continuous ice masses feed two or more glacier lobes, Nylen divided the glaciers based on ice flow directions derived from ice surface contours. Changes in glacier areas may reflect retreat or advance or simply changes in glacier boundaries on ice surfaces due to changes in ice flow directions. Nylen (2004) calculated glacier volumes from radar surveys (Driedger and Kennard, 1984) and empirical estimation methods (Bahr et al., 1997; Mennis, 1997; Mennis and Fountain, 2001). Accuracy of the outlines was limited by less accurate field mapping in hazardous terrain, and errors in distinguishing debris-covered ice from debris-covered ice-free areas, seasonal snow from snow-covered glaciers, active glacier ice from stagnant glacier ice, and glacier ice from permanent bodies of snow or ice, both in the field and with aerial photographs.

We created 2008 glacier outlines for this study by editing Nylen's (2004) 1994 glacier outlines to match glacier margins visible on the 2006 NAIP orthophoto and 2007 LiDAR intensity data (Copeland, 2009). This process is also similar to how historic glacier terminus positions were derived in Nylen's work. Glacier areas from the 2008 coverage are not comparable to Nylen's data, because our omission of glacier fragments will tend to underestimate glacier areas, while our possibly mistaken inclusion of snow cover at high elevations will tend to overestimate glacier areas.

Debris-flow initiation mechanisms were inferred through field reconnaissance, image interpretation, and LiDAR DEM differencing. For the 2006 Pyramid Creek debris flow, we measured gully widths at 50-m intervals before and after debris-flow initiation on the summer 2006 NAIP orthophoto and summer 2008 LiDAR from the initiation site to the downstream extent of the erosion zone, as snow cover permitted. Snow cover in the images hindered such measurements in all other cases. We also analyzed images of channels impacted by debris flows in 2001, 2003, and 2005 for evidence of erosion contributing to debris-flow initiation. Multiple dates of LiDAR data only exist for Tahoma Creek below South Tahoma Glacier. In addition to the 2008 LiDAR coverage of all of Mount Rainier National Park, LiDAR data were acquired from a 2-km swath around Tahoma Creek from South Tahoma Glacier to the Nisqually River confluence. The differences between the 2008 and 2000 DEMs show zones of erosion and deposition during that time interval.

In order to analyze the entire LiDAR DEM at once, the original 1-m data were decimated to produce a DEM with 10-m discretization. For visualization, we also used decimation to produce a DEM with 38-m discretization (Fig. 1). The decimated 10-m data were used for the data model presented here.

We present here an initial attempt to model the currently identified Mount Rainier debris-flow initiation sites and thereby predict the relative probability of debris-flow initiation at similar locations across the DEM when exposed to the same conditions as are considered in the model. This work should be viewed as a demonstration of the data-modeling concept, rather than as a final word on the prediction of future debris-flow initiation sites, as the choice of variables and combinations of variables used in the model can lead to potentially different results. While the model incorporates variables hypothesized to influence the process of debris-flow initiation, the model is based on observed relationships among those variables and not on a process-based formulation, as in Miller and Burnett (2007).

We have chosen four simple variables that are reasonable from a priori theoretical and empirical standpoints for this first modeling attempt:

• (1) Logarithm of the minimum distance (m) from any pixel within the 1994–2008 retreat zones.

• (2) Ground slope, or the arctangent of the magnitude of the gradient vector calculated from central differences in the x- and y-directions (degrees).

• (3) Logarithm of the average simulated glacial meltwater flow volume (m³) originating from the retreat zones, described below.

• (4) Logarithm of the standard deviation of simulated glacial meltwater flows (m³).

This model formulation focuses on debris-flow initiation by meltwater, but it may also capture initiation by runoff from rainfall on impermeable glacier ice. This approach encapsulates topographic and hydrologic factors, but it does not incorporate any material property variation. Distance from the retreat zones was chosen because recently exposed areas may have unconsolidated substrates susceptible to headward migration of gully-head knickpoints, although the model does not explicitly account for any variations in material properties. Ground slope was chosen because debris flows require steep slopes. Average simulated meltwater flow volume originating from retreat zones was chosen because glacier meltwater might be responsible for debris-flow initiation; runoff from rainfall on glaciers might also be responsible, and simulated meltwater flow volume should also identify areas more likely to receive focused flow from supra-, intra-, and subglacial drainage networks. Standard deviation of simulated meltwater flow volume was chosen because locations with lower, but still substantial, likelihood of meltwater flow or rainfall runoff might be more likely to supply unconsolidated material to debris flows when they occur. The variables were evaluated at each point in the decimated DEM, and each point with valid values for each of the four variables and not covered by glacier ice in 2008 (i.e., the set of “valid” points or potential debris-flow initiation sites) was included in the model calculations (Fig. 1).

The above data model requires some measure of the mean and standard deviation of meltwater flows originating from the glacier recession areas. We used a branching meltwater flow algorithm that allows meltwater flow to both nearest and next-nearest neighbors and thereby favors simplicity by avoiding the necessity of pit-filling or cascade algorithms (Braun and Sambridge, 1997; Tucker et al., 2001):

• (1) Create a drainage network for the entire DEM with the criteria that water at a grid node can flow to any nearest or next-nearest downhill neighbor, that is, water at node (i, j) may flow to nodes (m, n), m ∈ [i – 2, i + 2], n ∈ [j – 2, j + 2], (m, n) ≠ (i, j).

• (2) Identify discrete meltwater flow sources, as described below.

• (3) Identify the set of points downstream of each source with the drainage network and a recursive descent algorithm.

• (4) Sort the set of downstream points in order of decreasing elevation.

• (5) Calculate the meltwater flows within the set from top to bottom. Source points are set to a specific value (described below). All other points are evaluated in sequence, summing the flows from each upstream point. Meltwater flows from point i to point j are calculated as:

  

  where k ranges over all downslope nearest and next-nearest neighbors of point i, and S_ij is downhill-positive slope between points i and j:

  

  where xⁱ is the two-dimensional location of point i.

• (6) Sum the meltwater flows at all points from all sources.

With the above procedure, 25 simulations were performed. In each simulation, 1000 meltwater flow sources were chosen at random points within the entire melt zone, and volumes for each source were assigned from the volume of ice at that pixel as calculated from Nylen's (2004) glacier thicknesses (Table 2). This value was then multiplied by the grid-cell size and adjusted for the oblique surface area of each cell by dividing by the z-component of the surface normal at each point to yield pixel volumes in cubic meters. Note that the resulting volumes are not intended to represent actual volumes but rather the relative meltwater flows from each glacier according to its recession.

Twelve debris flows had initiation sources at the heads of 17 gullies distributed over seven distinct initiation zones near the termini of glaciers in six drainages in 2001, 2003, 2005, and 2006 (P. Kennard, 2008, personal commun.; Vallance et al., 2002), and field work to map the debris-flow initiation sites was conducted during the summer of 2008 (Table 3). All debris flows were triggered by rainstorms, with the exception of a debris flow on August 14, 2001, from Kautz Glacier that flowed into Van Trump Creek (Vallance et al., 2002). On October 20, 2003, with 0.3-cm snow-water equivalent (SWE) at the Paradise Ranger Station, 7.1 cm of rain in 24 h (71 mm/day) again produced a debris flow in Van Trump Creek. On September 29, 2005, with 0.0 cm SWE at Paradise, 16.0 cm of rain fell in 48 h (80 mm/day) and resulted in four debris flows, including one in Van Trump Creek. With 0.0 cm SWE at Paradise, the storm of November 5–7, 2006, delivered 41.7 cm in 48 h (209 mm/day) and triggered debris flows in six drainage basins. Four of these basins, all on the south side of the mountain, had debris flows in 2005. All of the debris flows occurring in 2001–2006 initiated in 17 gullies above tree line in steep, unconsolidated Quaternary-age volcanic and morainal material downstream of glacier ice or permanent bodies of snow and ice (Figs. 1 and 2). Gullies are fed and presumably maintained by glacier meltwater or snowmelt and appear to predate the first recorded debris-flow initiation. Gullies that sourced debris flows in 2006 were identified on images from 1989, the date of the oldest available park-wide aerial photos. For some debris-flow tracks, multiple tributaries converge upstream of the last observed evidence of a debris flow, and all tributaries show evidence of erosion, possibly due to debris-flow initiation or from normal fluvial processes. In these cases, a single initiation site cannot be determined. It is possible that all or only some of the meltwater channels that converge to form the main channel could have contributed to debris-flow initiation.

The steep slopes of Mount Rainier also contain numerous gullies that have not previously been associated with debris-flow initiation. We identified all gully walls as pairs of linear steep-sloped and approximately parallel features downslope of glaciers or snowfields and upslope of the mapped stream network. Gullies have arcuate steep-sloped features, or headwalls, and both gully sides and headwalls are visible on slope maps (Fig. 3A). They are typically located below glacier ice or permanent snow that is visible on LiDAR intensity images or at stagnant, debris-covered glacier ice that is not visible on LiDAR intensity images but has been described by previous researchers (P. Kennard, 2008, personal commun.; Walder and Driedger, 1994). On LiDAR intensity images, gullies may appear dark due to the presence of glacier meltwater, as does debris-free glacial ice (Fig. 3B). Points indicating the head of each gully were placed immediately downslope of the gully headwall or where the steep gully walls converged (Figs. 3A and 3B).

Many gullies were located, but only a few initiated debris flows (Fig. 3C). Gullies with and without a recorded history of debris-flow initiation had similar widths, depths, gully wall slopes, and channel gradients, but elevations and upslope contributing areas are statistically different for debris-flow–producing gullies and gullies without debris flows. Debris-flow–producing gullies have greater elevations (median of 2182 m with debris flows versus 1805 m without) and greater contributing areas (median of 0.184 km² versus 0.012 km²; Fig. 3D).

One of the six debris flows in 2006 initiated near a proglacial accumulation of water in the headwaters of the West Fork White River west of Winthrop Glacier, at an elevation not accessible during summer fieldwork due to extensive snow cover. Field identification of the West Fork of the White River debris-flow channel led to identification of the initiation site near a lake adjacent to a fragment of Winthrop Glacier, and the slope at this initiation site is much less than at all other initiation sites (Table 3; Fig. 4). The 2006 West Fork White River debris flow may have initiated as the result of a perturbation of this pool of water. The possibility of a glacier or glacier lake outburst flood contributing to debris-flow initiation cannot be confirmed or disregarded, because available images are not of a sufficient temporal resolution to reveal changes in volume of the proglacial lake. The lake terminates at two gullies that connect to the lateral margin of Winthrop Glacier.

Based on image analysis and field surveys, all recent (since 2001) debris-flow initiation sites on Mount Rainier are located in gullies headed by glaciers, ice bodies, or permanent snowfields. Therefore, all debris-flow initiation sites are located within the extent of the historical glacier coverage and within the zone of the most recent period of glacier retreat for a particular glacier. Debris-flow initiation sites associated with small glaciers that have retreated and fragmented, such as Pyramid and Van Trump Glaciers, are located within the 1913 extent of glacier coverage (Nylen, 2004). Debris-flow initiation sites associated with large glaciers that have retreated, but remained predominantly single, massive bodies of ice, such as South Tahoma, Kautz, and Winthrop Glaciers, are located within the 1971 extent of glacier coverage (Nylen, 2004). Lastly, a debris flow initiated at a medium size glacier that has retreated minimally, Inter Glacier, where initiation sites are located within the 1994 extent of glacier coverage (Nylen, 2004).

Distance from the head of the gully that contributed to debris-flow initiation to the 2008 glacier or glacier remnant is zero for all initiation sites except those of South Tahoma, Kautz, and Winthrop Glaciers. The gullies that contribute to debris flows in Tahoma Creek are located on stagnant ice, over 100 m from the active terminus of South Tahoma Glacier, deposited by glacier retreat past a series of bedrock steps (C. Driedger, 2008, personal commun.). The large, debris-covered, stagnant ice terminus of Kautz Glacier (Nylen, 2004) is nearly 500 m downslope of the active glacier terminus and is incised by a channel that drains glacier meltwater and may have contributed to the 2006 debris flow in Kautz Creek. The debris flow that occurred in 2006 in the West Fork of the White River initiated below a fragment of the 1913 extent of Winthrop Glacier. The gully is separated from the fragment of Winthrop Glacier by a proglacial lake ∼200 m in length.

Debris flows are only located in areas exposed by glacier retreat in the past century. The three glacier termini that have retreated the farthest during the period 1913–1994 are those of South Tahoma, Nisqually, and Kautz Glaciers, with 2677, 2131, and 1589 m retreat along the midline of the glacier termini, respectively (Nylen, 2004). These glaciers have been associated with more debris flows during the period 1926–2006 than any other glaciers on Mount Rainier, or 28, eight, and six debris flows, respectively (including two, zero, and three, respectively, in 2001–2006; Table 3; Walder and Driedger, 1994). South-facing glaciers, such as South Tahoma, Nisqually, and Kautz, have been associated with debris flows more often than north-facing glaciers since 1926, and the average retreat of south-facing glaciers from 1913 to 2004 was 1755 m, compared to 530 m for north-facing glaciers (Nylen, 2004). Yet, it is important to note that since 1980, the termini of all glaciers on Mount Rainier have retreated, with the exception of Emmons, Winthrop, and Cowlitz Glaciers (Nylen, 2004). The termini of the majority of Mount Rainier glaciers have retreated during the period 1913–1994, and few glaciers have experienced debris flows. Glacier terminus position change, or cumulative retreat and advance, from 1913 to 1994 using data from Nylen (2004) for glaciers associated with debris flows and glaciers not associated with debris flows during the same period is not significantly different (Welch Two Sample t-test: t-statistic = −1.32, df = 3.65, p-value = 0.264).

Debris flows that initiated in 2006 are predominantly associated with small glaciers or glacier fragments less than 2.3 km² in size. The difference between the 2008 glacier area for glaciers associated with debris flows during the period 1994–2008 and glaciers not associated with debris flows during the same period is significant (Welch Two Sample t-test: t-statistic = −3.57, df = 19.7, p-value = 0.002). The average size of a glacier associated with a 2006 debris flow is 1.1 km², while the average size of all glaciers on Mount Rainier is 3.8 km².

All initiation sites from 2006 are located at elevations greater than 1937 m. From records dating back to 1926, we see that debris flows initiated predominantly at glaciers with termini located at lower elevations, such as those of South Tahoma, Kautz, and Nisqually Glaciers (Copeland, 2009; Walder and Driedger, 1994). Debris flows that occurred in 2001, 2003, 2005, and 2006 occurred at higher-elevation glacier termini such as Pyramid and Van Trump Glaciers, as well as fragments of glaciers located at high elevation, such as Winthrop Glacier (Figs. 1 and 2).

All debris flows, with the exception of the 2006 debris flow from Winthrop Glacier, initiated from gullies with channel slopes ≥25° and ground slopes ≥14° (Table 3). Average gradient for a 100-m radius circle surrounding each debris-flow initiation site is 32°. The average debris-flow channel gradient at initiation for the first ∼9 m of the longitudinal profile of the debris-flow track is 39°. Both types of calculated gradients are important for debris-flow initiation because channel gradient will aid sediment mobilization and the steepness of surrounding terrain will aid overland flow concentration. Nearly all proglacial areas exposed by recent retreat contain slopes sufficient for debris-flow initiation.

Field reconnaissance revealed evidence of mobilization of sediment in rills and gullies associated with debris flows. On the images, all debris-flow channels were traced to meltwater gullies with crenulated walls, possibly indicative of discrete mass wasting or slumping into the stream channel in association with localized overland flow or in-channel entrainment of sediment. Following the 2006 debris flow in Pyramid Creek, the average and median changes in gully width were 0.7 m and 0.0 m, respectively. Although other events after the November 2006 storm and prior to 2008 may have eroded gully walls, the extraordinary magnitude of the 2006 event implies that changes in gully width can be attributed to debris-flow activity in 2006 with reasonable certainty. Image analysis of channels impacted by debris flows in 2001, 2003, and 2005 also suggests that slumping from steep-sided gully walls facilitated debris-flow initiation. Differencing of the 2000 and 2008 LiDAR DEMs for Tahoma Creek shows that erosion was distributed over a 500-m-wide proglacial area and stagnant ice zone described by Walder and Driedger (1994) downslope of the terminus of South Tahoma Glacier, where glacier retreat over the period appears negligible (Figs. 1 and 5). This swath of erosion converges to a single channel 700 m downslope of the terminus; below this channel head, deposition predominates for nearly 200 m; next, erosion is concentrated in a discrete channel for more than 250 m; the erosional zone then appears to cover a 50–150-m-wide valley floor for the next 1300 m; erosion is then apparently limited to the valley margins, whereas deposition occurs along the valley center. Differences indicative of erosion near the terminus of Tahoma Glacier appear largely coincident with areas of glacier retreat in 1994–2008, and differences in areas at the head of the proglacial streams, where gullying might be expected, indicate negligible erosion. The DEM differencing evidence is consistent with debris-flow initiation over a widely distributed network of rills and gullies, although bulking by erosion of the bed and banks apparently contributed significant volume downstream of the initiation zone (Fig. 5).

Slopes are shown in Figure 2B. Average simulated meltwater flows are shown in Figure 6A, and the results of an example meltwater flow simulation are shown in Figure 6B. Standard deviations of simulated meltwater flows are shown in Figure 6C. Statistics of the distributions of distance from the nearest retreat zone, average simulated meltwater flow, standard deviation of simulated meltwater flow, and slope for 1,422,607 valid tuples (i.e., points downstream of glaciers) are shown in Table 4. The full distributions of the model variables and the values at the debris-flow initiation sites show that the latter generally cluster within limited ranges of each parameter (Fig. 4) and form a distinct cluster in the four-dimensional space defined by the model variables (Fig. 7). Identification of susceptible areas in terms of their proximity to these points in the four-dimensional data space with Equation (1) unsurprisingly highlights the areas adjacent to and downstream of glacier termini (Fig. 8). Given the formulation of the model, all sites of known debris-flow initiation are of course identified as high-susceptibility areas, as are the heads of most gullies (Fig. 3). With respect to gullies, the model does not appear able to distinguish particularly well between gullies that initiated debris flows and gullies that did not, nor does the model distinguish between areas with greater and lesser gully head concentrations.

The locations of the initiation sites within the chosen parameter space do clearly identify one outlier among them: The initiation site of the West Fork White River debris flow near Winthrop Glacier has a much lower slope than all other initiation sites. This difference suggests an initiation mechanism different than that of the other debris flows.

In the effort to quantitatively determine the hazards represented by debris flows and their associated risks where human life and property are at stake, it is necessary to first make a quantitative, spatially distributed determination of relative susceptibilities to debris-flow initiation, erosion, and deposition, of which we address only the susceptibility to initiation here. Subsequent studies can then determine likely extents of erosion and deposition; potentials of storms of different intensities and durations to cause debris-flow initiation in susceptible areas; probabilities of such storms striking susceptible areas when conditions (e.g., antecedent snowfall, rainfall, and vegetation) favor debris-flow initiation; and probabilities of debris flows destroying human life and infrastructure. The ideal goal of debris-flow susceptibility mapping is to include all areas that actually represent greater susceptibility within the area mapped as high susceptibility and for the mapped high-susceptibility area to include no areas actually associated with low susceptibility. For Mount Rainier, where recent debris flows have all initiated in areas uncovered by glacial recession within the past century or even, in many cases, within the past 20 yr, we would also like to know whether the uncovering of those areas has significantly increased the probability of debris-flow initiation on the mountain. The modeling approach presented here has the potential to address both of these issues, but the current formulation of the data model is fundamentally unable to answer the latter question and appears to fall short of the ideal with respect to the former goal as well. The current model does, however, illustrate the potential power of the data-modeling method and, with our detailed investigations of initiation zones, improve our understanding of the mechanisms involved in debris-flow initiation.

A data model ideally identifies areas that are most similar to known sites based on their most important characteristics. In this study, we wanted to choose those characteristics that made the sites most prone to debris-flow initiation. While data modeling specifically avoids formulation of a physics-based process model, the method requires accurate determination of the variables required by such a process model. Our results suggest that the current model does not employ the correct set of variables. While the model does appear to correctly identify high-susceptibility areas, it also appears to misrepresent lower-susceptibility areas: non-debris-flow gullies are not well differentiated from debris-flow gullies, and areas of lower gully density are not well differentiated from areas of higher density (Figs. 3 and 8).

Despite any shortcomings, the current model and field data provide insight into relevant processes and suggest a way forward. Field observations and mapping indicate that debris flows begin as overland water flows, typically from glacier margins, that entrain sediment from beds and walls of gullies and perhaps even from upslope rills (Figs. 3 and 5). For the model, we modeled meltwater flows as though they were emanating from locations of glacier retreat, but it seems more likely that the actual overland water flows associated with debris-flow initiation are sourced by rainfall on, or melt from, the entire upslope glacier area, as well as other impervious surfaces such as bedrock. This conclusion is consistent with two results: (1) as noted above, the model implies similar debris-flow susceptibility in areas with and without gullies and for gullies with and without debris flows; (2) we found that gullies with debris flows have greater upslope contributing areas than gullies without debris flows (Fig. 3). The latter fact indicates that the gullies themselves are maintained by fluvial processes, whereas erosion within, or erosive expansion of, the gully network leads to the initiation of debris flows. Whether the gullies are originally formed by fluvial or debris-flow processes remains an open question.

Whereas debris flows initiating from shallow rapid landslides typically require slopes greater than 20° (40° is typical of landslide slopes in the Oregon Coast Range; Robison et al., 1999), debris-flow initiation sites on Mount Rainier had slopes ranging from 2°, for a gully originating at a proglacial lake below a fragment of Winthrop Glacier, to 38°, according to measurements from the decimated 10-m DEM. Note that gully slopes measured from the 1-m DEM averaged 39° and covered the range 3° to 65° (Table 3). This range of slopes indicates that, while debris flows generally favor steep slopes, other factors, such as overland water flow, are more determinate of the likelihood of debris-flow initiation at a point. The measured changes in elevation downstream of South Tahoma Glacier and gully widths at all initiation sites support the inferences that the mapped debris flows are initiated by combinations of bed and bank erosion and bank failure during high overland flow events, and these debris-flow sources are distributed over wide areas rather than concentrated within discrete landslides.

An improved data model would likely result from a different set of variables. Slope appears somewhat diagnostic and should therefore remain in the model. Water discharge appears to play a large role and should appear in the model in some form, as either total upslope contributing area or impervious or ice-covered upslope contributing area. Some measure of the gullies themselves would likely enhance the model. Even if gully dimensions are not determinate, and presence of gullies is not sufficient for debris-flow initiation, they do appear to be necessary. The decimated 10-m DEM is unlikely to capture gully features because the gullies are typically ∼10 m wide (Fig. 3). Gullies might be included by calculating topographic curvature (∇²z). This and other variables should also be calculated on the 1-m DEM, and calculation of contributing area should employ routing through closed depressions (Tucker et al., 2001). Given that all mapped debris flows initiated close to glacier margins in areas recently uncovered by glacier retreat, some measure of glacier proximity might be appropriate. Rather than minimum distance to an area of recent glacier retreat, minimum upslope distance to an existing glacier might be more appropriate.

The current method does not allow optimization to maximize correct identification of high-susceptibility sites and to minimize incorrect identification of low-susceptibility sites. For the current type of formulation using distance in the four-dimensional model space to points representing known debris-flow initiation sites, we would ideally create a model with a subset of the mapped initiation sites and verify or optimize it with the remaining subset. Such a partitioning of the data would allow evaluation of the ratio of the incremental fraction of points correctly identified among the data used for optimization to the incremental fraction of area identified as susceptible as a function of, say, a threshold distance from the data used in model creation in the four-parameter space (e.g., Miller and Burnett, 2007). Or, perhaps the ground-based data could be used to create isosurfaces of site concentration in the parameter space, or distance could be calculated from the geometric center of the known initiation sites rather than the minimum distance to a single site in the parameter space.

Debris flows on Mount Rainier during the period 2001–2006 began when and where unusual amounts of overland flow of water led to distributed erosion associated with a proglacial gully network. An event in 2001 was triggered by redirection of summer meltwater, but most events, occurring in 2003, 2005, and 2006, were triggered by intense rainfall at high elevations with little or no antecedent snow cover. The debris flows originated at 17 different gully heads on predominantly steep topography in close proximity to, and historically exposed by, retreating glaciers. The locations of these gully heads were determined from field mapping and a high-resolution DEM (with 1-m grid cells) derived from airborne laser altimetry (LiDAR) collected in 2007 and 2008. Gully locations were confirmed by low intensities of infrared laser returns from gully bottoms, where low return intensities are consistent with the presence of water on the surface. We identified many gullies not associated with debris flows, and these are distinguished from debris-flow gullies by the greater upslope contributing area of the latter (Fig. 3). The DEM, laser intensity map, and orthophotos (from 2006) were used to revise glacier outlines previously mapped for 1994 and show that the total area of Mount Rainier's glaciers has decreased since that time. The debris flows were effectively initiated by overland flow of water, where that flow led to widely distributed erosion of beds and banks of gullies. We cannot determine the relative contributions of progressive scour and mass failure, and scalloped gully walls offer direct evidence of only the latter mechanism.

As a demonstration of the method, we developed a data model in which the values of four variables, topographic slope, mean and standard deviation of simulated meltwater flow from areas of glacier retreat, and minimum distance from areas of glacier retreat, were calculated at all points downstream of the glaciers, and relative debris-flow initiation susceptibility was calculated as a decreasing function of the minimum distance between each point and known sites of debris-flow initiation in the four-dimensional model space. This data modeling method may prove useful in a revised form, but the current model does not appear to distinguish between debris-flow and non-debris-flow gullies or between areas of varying areal density of gullies. Although debris flows initiated in areas uncovered by glacier recession, we cannot determine whether such areas have greater susceptibility due to material properties (e.g., fresh unconsolidated debris) or overland flow hydraulics associated with glacial outflow.

Paul Kennard (Regional Geomorphologist, U.S. National Park Service) and the staff of Mount Rainier National Park provided generous support and guidance of fieldwork efforts. The LiDAR data were made possible by the Puget Sound LiDAR Consortium and Mount Rainier National Park. Watershed Sciences, Inc., acquired and processed the airborne LiDAR data and provided expert assistance. Sky Coyote (Sky Coyote Science Software, www.skycoyote.com) wrote the software programs constituting the data model and performed the calculations for that model. Thomas Sisson (U.S. Geological Survey, Cascades Volcano Observatory) recommended the methodology for creating the 2008 glacier outlines. Carolyn Driedger (U.S. Geological Survey, Cascades Volcano Observatory) generously provided unpublished information regarding debris flows on Mount Rainier. Sarah L. Lewis (College of Earth, Ocean, and Atmospheric Sciences, Oregon State University) provided valuable feedback and skilled editing of the manuscript. This work was supported by grants from the National Science Foundation's Geomorphology and Land Use Dynamics Program (EAR-0756825 and EAR-0844017).