A regional spatio-temporal analysis of large magnitude 1 snow avalanches using tree rings 2

Snow avalanches affect transportation corridors and settlements worldwide. In many mountainous 10 regions, robust records of avalanche frequency and magnitude are sparse or non-existent. However, 11 dendrochronological methods can be used to fill this gap and infer historic avalanche patterns. In this study, 12 we developed a tree-ring based avalanche chronology for large magnitude avalanche events using 13 dendrochronological techniques for a portion of the northern United States Rocky Mountains. We used a 14 strategic sampling design to examine avalanche activity through time and across nested spatial scales (i.e. 15 from individual paths, four distinct sub-regions, and the region). We analysed 673 total samples from 647 16 suitable trees collected from 12 avalanche paths, from which 2,134 growth disturbances were identified over 17 years 1636 to 2017 Common Era (C.E.). Using existing indexing approaches, we developed a regional 18 avalanche activity index to discriminate avalanche events from noise in the tree-ring record. Large magnitude 19 avalanches common across the region occurred in 30 individual years and exhibited a median return interval 20 of approximately three years (mean = 5.21 years). The median large magnitude avalanche return interval (321 8 years) and the total number of avalanche years (12-18) vary throughout the four sub-regions, suggesting 22 the important influence of local terrain and weather factors. We tested subsampling routines for regional 23 representation, finding that sampling eight random paths out of a total of 12 avalanche paths in the region 24 captures up to 83% of the regional chronology, whereas four paths capture only 43% to 73%. The greatest 25 value probability of detection for any given path in our dataset is 40% suggesting that sampling a single path 26 would capture no more than 40% of the regional avalanche activity. Results emphasize the importance of 27 sample size, scale, and spatial extent when attempting to derive a regional large magnitude avalanche event 28 chronology from tree-ring records. 29 30 1 https://doi.org/10.5194/nhess-2020-253 Preprint. Discussion started: 20 August 2020 c © Author(s) 2020. CC BY 4.0 License.

Snow avalanches are hazardous to human safety and infrastructure (Mock et al., 2016;Schweizer, 2003) as 33 well as an important landscape disturbance affecting mountain ecosystems (Bebi et al., 2009). In the United 34 States an average of 27 people die in avalanche accidents each winter (CAIC, 2020). Avalanches, especially 35 large magnitude events, also affect transportation corridors and settlements throughout the world. For 36 example, avalanches impact numerous roadways and railroad corridors in the western United States 37 (Armstrong, 1981;Hendrikx et al., 2014;Reardon et al., 2008). Consequently, understanding general 38 avalanche processes and associated large magnitude avalanche return intervals is critical for local and 39 regional avalanche forecasters, transportation agencies, and land use planners. 40 Long-term, reliable, and consistent avalanche observation records are necessary for calculating avalanche 41 return intervals, which can be used in infrastructure planning and avalanche forecasting operations. However, 42 such records are often sparse or non-existent in many mountainous regions, including areas with existing 43 transportation corridors. Thus, inferring avalanche frequency requires the use of dendrochronological 44 methods to document damaging events or geomorphic response within individual trees at individual path to 45 regional scales. Even in regions with historical records, tree-ring dating methods can be used to extend or 46 validate uncertain historical avalanche records, which has led to the broad implementation of these methods 47 in mountainous regions throughout the world (e.g. Corona  Numerous studies reconstructed avalanche chronologies in the United States using tree-ring methods 50 (Burrows and Burrows, 1976;Butler et al., 1987;Carrara, 1979;Hebertson and Jenkins, 2003;Potter, 1969; 51 Rayback, 1998). Butler (Table A1). 58

Framework and objectives 59
Tree-ring avalanche research is resource and time intensive. Like other scientific fields, it is not feasible to 60 completely sample the variable of interest with infinite detail due to logistical and financial constraints 61 (Skøien and Blöschl, 2006). Thus, a strategic spatial sampling method is necessary. Here, we strategically 62 sampled 12 avalanche paths in four distinct sub-regions of the U.S. northern Rocky Mountains of northwest 63 Montana to examine spatial differences at a regional scale. The sampling strategy is based on the concept of 64 scale triplet, which defines the spacing, extent, and support of our sampling scheme (Blöschl and Sivapalan, 65     (Mock and Birkeland, 2000), but it can exhibit characteristics of both continental or coastal climates. 121 The elevation of avalanche paths within the study sites range from approximately 1100 m to 2700 m and the 122 starting zones of these paths are distributed among all aspects (Table 1). 123 We eliminated or minimized influence from exogenous disturbance factors such as logging and wildfire by 124 referencing wildfire maps extending back to the mid-20 th century. We selected sites undisturbed by wildfire 125 since this time except for Lost Johnny Creek, which was purposeful as this area burned most recently in 2003. 126 We also minimized the influence of logging by selecting sites not previously logged. Using historical logging 127 parcel spatial data, we determined logging in some sites was limited to very small parcels adjacent to the 128 farthest extent of the runout zones. 129 The historical observational record in this area is limited. In this study region, the Flathead Avalanche Center 130 (FAC), a regional U.S. Forest Service backcountry avalanche center, records all avalanches observed and 131 Transportation and railroad company records, National Park Service ranger logs, and popular media archives. 138 In this sub-region, avalanche mitigation is conducted on an infrequent and inconsistent basis in emergency 139 situations, which is typically only once a year, if at all. Thus, the record approximates a natural avalanche 140 record. We compared the reconstructed avalanche chronology of the JFS sub-region to the historical record 141 of large magnitude years for qualitative purposes. A quantitative comparison would not be reflective of the 142 true reliability of tree-ring methods because of the incomplete historical record. 143

Sample Collection and Processing 144
Our sampling strategy targeted an even number of samples collected from both lateral trimlines at varying 145 elevations and trees located in the main lower track and runout zone of the selected avalanche paths. This 146 adequately captured trees that were destroyed and transported, as well as those that remained in place. The 147 definition of large magnitude avalanche in this study refers to avalanches of approximately size D3 or greater 148 (Greene et al., 2016) and may not run the full length of the avalanche path. We sampled spatial extents within 149 each avalanche path representative of runout extents ≥ size D3 avalanches. We also used recent (within 150 previous 10 years) observed large magnitude avalanche activity in these paths to constrain our sampling. 151 Sample size for avalanche reconstruction using tree-ring data requires careful consideration. Butler and 152 Sawyer (2008) suggested that a few damaged trees may be sufficient for avalanche chronologies, but larger 153 target sample sizes increase the probability of detecting avalanche events (Corona et  increases in the probability of extending chronologies with sample size greater than 40. Thus, given the large 156 spatial footprint (~3500 km 2 ) of this study and feasibility of such a large sample size, we sampled between 157 26-109 samples per avalanche path resulting in 655 trees (Table 2). Eight trees were unsuitable for analysis  158   leaving us with 673 total samples from 647 trees. Of the 673 total samples, we collected 614 cross sections  159 and 59 cores. Shed 10.7 (S10.7) path was the focus of previous work (Reardon et al., 2008), and the 160 dendrochronological record extends up to 2005 (n=109 trees). Little Granite Path (LGP) was collected in the 161 summer of 2009 (n=109 trees). We sampled the remaining 10 paths (437 of the 655 total trees) in the summer 162 of 2017. 163 We collected three types of samples: (1) cross sections from dead trees, (2) cross sections from the dead 164 leaders of avalanche-damaged but still living trees, and (3) cores from living trees. We used predominantly 165 cross-sections in this study for a more robust analysis as events can potentially be missed or incorrectly 166 identified in cores. We emphasized the selection of trees with obvious external scars and considered location, 167 size, and potential age of tree samples. A limitation of all avalanche dendrochronology studies is that large 168 magnitude events cause extensive damage and high tree mortality, thereby reducing subsequent potential 169 tree-ring records. 170 We sampled stem cross-sections at the location of an external scar or just above the root buttress from downed 171 or standing and dead trees, and from stems of trees topped by avalanche damage. We extracted tree-ring core 172 samples from live trees with obvious scarring or flagging along the avalanche path margins and runout zone 173 using a 5 mm diameter increment borer. We collected a minimum of two and up to four core samples per tree 174 (two in the uphill-downhill direction and two perpendicular to the slope). We photographed each sample at 175 each location and recorded species, Global Positioning System (GPS) coordinates (accuracy 1-3 m), amount 176 of scarring on the cambium of the tree, relative location of the tree in the path, and upslope direction (Peitzsch 177 et al., 2019). We also recorded location characteristics that identified the tree to be in-place vs. transported 178 from its original growth position (i.e. presence or absence of roots attached to the ground or the distance from 179 an obvious excavated area where the tree was uprooted). 180 To prevent radial cracking and further rot, we dried and stabilized the cross sections with a canvas backing. 181 We sanded samples using a progressively finer grit of sandpaper to expose the anatomy of each growth ring, 182 and used the visual skeleton plot method to account for missing and false-rings and for accurate calendar 183 year dating (Stokes and Smiley, 1996). We assessed cross-dating calendar-year accuracy of each sample and 184 statistically verified against measured samples taken from trees within the gallery forest outside the avalanche 185 path, and from preexisting regional chronologies (Table A2)

Avalanche Event Identification 189
We analyzed samples for signs of traumatic impact events (hereafter "responses") likely caused by snow 190 avalanches. We adapted a classification system from previous dendrogeomorphological studies to 191 qualitatively rank the trauma severity and tree growth response from avalanche impacts using numerical 192 scores ranked 1 through 5 (Reardon et al., 2008). This classification scheme identified more prominent 193 avalanche damage responses with higher quality scores, and allowed us to remain consistent with previous 194 work (Corona et al., 2012;Favillier et al., 2018) (Table 2). To compare our ability to capture 195 avalanche/trauma events using cores versus those captured using cross-sections, we sampled a subset (n=40) 196 of the cross-sections by analyzing four 5 mm wide rectangles to mimic a core sample from an increment 197 borer. The four subsamples on each cross section were made perpendicular to one another (i.e. 90 degrees) 198 based on the first sample taken from the uphill direction of each stem to replicate common field sampling 199 methods. We then summarized results from the four subsamples for each tree by taking the highest response 200 score for each growth year. Finally, we compared the number, quality response category, and calendar year 201 of the avalanche/trauma events derived from the core subsamples to those identified from the full cross 202 sections. 203 204

C1
• Clear impact scar associated with well-defined reaction wood, growth suppression or major traumatic resin duct development. • Or, the strong presence of some combination of these major anatomical markers of trauma and growth response recorded in multiple years of growth and occurring at a year that multiple samples from other trees at the site record similar trauma and scaring. • C1 events are also assigned to the death date of trees killed by observed avalanche mortality at the collection site; the presence of earlywood indicates an early spring, or late avalanche season, event killed the tree.

C2
• Scar or small scar recorded in the first ten years of tree growth without associated reaction wood, growth suppression or traumatic resin ducts. • Or, obvious reaction wood, growth suppression or significant presence of traumatic resin ducts that occur abruptly after normal growth that lasts for 3 or more years.

C3
• The presence of reaction wood, growth suppression, or traumatic resin ducts recorded in less than 3 successive growth years.

C4
• Poorly defined reaction wood, growth suppression or minimal presence of traumatic resin ducts lasting 1-2 years. • Or, a C3 class event occurring in the first 10 years of tree growth where the cause of damage could result from various biological and environmental conditions.

C5
• Very poorly defined reaction wood, growth suppression, or minimal presence of traumatic resin ducts isolated in one growth year. • Or, a C4 class event occurring in the first 10 years of tree growth where the cause of damage could result from various biological and environmental conditions. 207

Chronology and Return Period Calculation 208
To generate avalanche event chronologies and estimate return periods for each path and for the entire study 209 site, we utilized R statistical software and the package slideRun, an extension of the burnR library for forest 210 fire history data (Malevich et al., 2018). We calculated the age of each tree sampled and the number of 211 responses per year in each avalanche path and computed descriptive statistics for the entire dataset. Estimates 212 of avalanche path return intervals should be viewed as maximum return interval values due to the successive 213 loss of samples and decreasing sample number back through time. 214 We used a multi-step process to reconstruct avalanche chronologies on three different spatial scales: 215 individual paths, four sub-regions, and the entire region. We also calculated a regional avalanche activity 216 index (RAAI) (Figure 2). The process involved first calculating the ratio of trees exhibiting growth 217 disturbance (GD) over the number of samples alive at year t to provide the index It (Shroder, 1978): 218 where R is the number of trees recording a GD at year t with At representing the number of trees alive in our 219 samples at year t. 220

For each path, sub-region, region
Filtering Events based on classification of events (1-5)

For each sub-region and region
Patterns/Trends of RAAI POD for an avalanche year POD using any one path  POD year (Eq. 4) The probability of detecting an avalanche year in the regional chronology by sampling any given individual path. POD path (Eq. 5) The probability of detecting the full regional chronology using any one given path's chronology. (2) where the sum of trees with scars or injuries (C1 -C5) were multiplied by a factor of 7, 5, 3, 1 and 1 236 respectively (Kogelnig-Mayer et al., 2011). 237 Next, we classified Wit into high, medium, and low confidence events using the thresholds detailed in discriminating the avalanche response from noise. We included all events with medium to high confidence 240 in the next analysis. We then estimated the number of avalanche years, descriptive statistics for return 241 intervals (RI), and the annual probability (1/RI) for each path, sub-region, and region. We use these RI 242 derived after filtering events for confidence as the intervals throughout the study. We then compared return 243 intervals for all individual paths and sub-regions using analysis of variance (ANOVA) and Tukey's Honest 244 Significant Difference (HSD) (Ott and Longnecker, 2016). In the final step of RI analysis, we subset the 245 period of record for each path from 1967-2017 to compare RI from this condensed time series to the full 246 period of record for each path. 247 Next, we compared the number of avalanche years and return periods identified in the full regional 248 chronology to subsets of the region to determine the number of paths required to replicate a full 12-path 249 regional chronology. We assessed the full chronology against a subsampling of 11 total paths by sequentially 250 removing the three paths with the greatest sample size. We then randomly sampled two paths from each sub-251 region for a total subsample of eight paths, followed by generating a subsample of four paths by choosing 252 the path in each sub-region with the greatest sample size. Finally, we selected a random sample of one path 253 from each sub-region to compare against a total of four single path subsamples. 254

Regional Avalanche Activity Index and Probability of Detection 255
Next, we used the It statistic from each path to calculate a regional avalanche activity index (RAAI) for the 256 sub-regions and overall region (Germain et al., 2009). The RAAI for each year across the sub-regions and 257 region provides a more comprehensive assessment of avalanche activity within the spatial extent. For each 258 year t, we calculated RAAI: 259 where I is the index factor as per Eq. (1) for a given avalanche path for year t and P is the number of paths 260 that could potentially record an avalanche for year t. For the calculation of the overall RAAI, we required 261 each path to retain a minimum sample size of ≥ 10 trees with a minimum number of three paths for year t, 262 and a minimum of one path from each sub-region. We performed a sensitivity test to establish the minimum 263 number of paths necessary to calculate an RAAI value for any given year. 264 We also calculated the probability of detecting an avalanche year identified in the regional chronology as if 265 any given individual path was sampled. The probability of detection for a given year (PODyear) is defined as: 266 where a is the number of individual avalanche paths that identify any given avalanche year in the regional 267 chronology and b is the total number of avalanche paths (n=12). We calculated PODyear for every year in the 268 regional avalanche chronology. We then compared the PODyear of individual paths to the number of active 269 avalanche paths as defined in Eq. (3). 270 We also calculated the probability of detection for each path for the period of record (PODpath): 271 where c is the number of years identified in any given path that is included in the regional chronology and d 272 is the number of years in the regional chronology that are not identified in the chronology for the given path. 273 Finally, we examined trends in the RAAI through time using the non-parametric modified Mann-Kendall test 274 for trend (Mann, 1945;Hamed and Rao, 1998). We parsed the dataset into four periods to allow comparability 275 due to the loss of evidence and a decreasing sample size going back in time: the entire period of record, 1933 276 to 2017, 1950 to 2017, and 1990 to 2017. We selected these time periods based on the years with greatest 277 responses and peak RAAI values (1933, 1950, and 1990). We excluded intervals after 1990 to retain a 278 minimum of ~30-year record. 279

Geomorphological characteristics 280
Using a 10 m digital elevation model (DEM), we calculated a number of geomorphological characteristics 281 for each path, including mean elevation (m, full path and starting zone), elevation range (m), eastness 282 (sin( )) and northness (cos( )) (radians), slope (degrees, full path and starting zone), curvature 283 (index (0-1), profile and planform), roughness (index, full path and starting zone), perimeter (km 2 ), area 284 (km 2 ), length (m), and vertical distance from starting zone to runout zone (m). We also calculated the mean 285 of these characteristics for all paths in the region. The geomorphological characteristics allowed for a 286 determination of the representativeness of the region as a whole (i.e. are the paths similar across the region?) 287 as well as a comparison of the return interval for each path relative to these characteristics. Finally, we 288 estimated the potential relationship between path length, starting zone slope angle, the number of avalanche 289 years, and median return interval for each individual path using the Pearson correlation coefficient. 290

Sub-region Chronologies 360
When the paths were aggregated into sub-regions (three paths per sub-region) the median return periods for 361 each sub-region were similar and all less than 10 years (Figure 6(e) and Table 4). The number of avalanche 362 years for all of the sub-regions ranges from 12-18 with the greatest number of identified years in the JFS sub-363 region and the fewest in the WF sub-region. The JFS sub-region has the shortest median return interval 364 followed by the Swan, WF, and GTSR sub-regions. The number of avalanche years for each aggregated sub-365 region is greater than the number of avalanche years for any individual path within each sub-region except 366 for the JFS sub-region where 18 avalanche years were identified but Shed 10-7 totaled 20 avalanche years 367 (

Regional Chronology and RAAI 387
We identified 30 avalanche years in the overall region and a median return interval of 3 years (Table 5). The 388 number of samples increases through time to a peak during 2005 and as expected the number of GD also 389 increases through time (Figure 8(a)). The Wit index also increases, particularly from year 2000 onward with 390 the largest spikes in 2014 and 2017 (Figure 8(b)). The regional assessment of avalanche years identified 391 fewer years (n=30) than the simple aggregation of all unique avalanche years identified in the individual 392 paths (n=49) ( Table A2).

401
When we included all paths but S10.7 (one of two paths with the greatest sample size), we captured 80% of 402 all avalanche years and added one new year to the chronology (Table 7). When we removed LGP (the other 403 path with the greatest sample size), we still captured all of the years in the regional chronology but introduced 404 four new years into the chronology for a total of 34 years. A random sample of eight (two from each sub-405 region) of the 12 avalanche paths captured 83% of the years in the chronology and identified two new 406 avalanche years. Finally, when using only one path from each sub-region with the largest samples size (Shed 407 10-7, 54-3, LJA, and RMA), we captured 73% of the avalanche years identified in the full regional 408 chronology. When using a random sample of one path from each sub-region (1163, LGP, LJC, RMB), we 409 captured only 43% of the years included in the regional chronology of all 12 paths. The RAAI is insensitive 410

418
Paths Region (All Paths) All but S10.7 To assess potential long-term trends in regional avalanche activity we implemented the modified Mann-420

All but LGP
Kendall test since the chronology exhibits weak serial autocorrelation. The full period of record of RAAI 421 The probability of detection for the avalanche years (PODyear) identified in the regional chronology ranged 425 from 8 to 58% when we examined individual paths (Table 7). The year with the highest POD was 2014. The 426 mean POD for all years was 21%. When we examined avalanche paths that exhibited at least one scar during 427 avalanche years identified in the regional chronologies, the POD is generally greater. 428 429  67  1970  17  50  1971  25  50  1972  33  83  1974  33  75  1976  17  50  1982  42  92  1990  42  83  1993  17  50  1997  8  92  1998  17  50  1999  17  58  2002  25  75  2003  17  33  2004  17  75  2009  17  33  2011  8  33  2012  17  42  2014  58  58  2017  17  25  Mean  21  52  432 Finally, the probability of capturing all of the avalanche years identified in the regional chronology by each 433 individual path ranges from 7% to 40% (Table 8). The greatest PODpath value from any given path is S10.7 434 (POD = 40%) in the JFS sub-region followed by RMC in the Whitefish sub-region (POD = 37%). In general, 435 the paths within the Whitefish sub-region capture the regional chronology most consistently. 436 437

Discussion 439
The processing and analysis of 673 samples spanning a large spatial extent allowed us to create a robust 440 regional large magnitude avalanche chronology reconstructed using dendrochronological methods. Cross-441 sections provided a more robust and complete GD and avalanche chronology compared to a subsample 442 generated from cores alone. Due to the reduced information value of working only with cores, Favillier et al. 443 (2017) included a discriminatory step in their methods to distinguish avalanche signals in the tree-ring record 444 from exogenous factors, such as abnormal climate signals or response to insect disturbance. By using cross 445 sections to develop our avalanche chronologies, we were able to view the entire ring growth and potential 446 disturbance around the circumference of the tree as opposed to the limited view provided by cores. This 447 allowed us to place GD signals in context to both climate and insect disturbance without the need for this 448 processing step. 449 We targeted sample collection in areas in the runout zones and along the trim line where large magnitude 450 avalanches occurred in recent years. However, at several sites we also collected samples into the bottom of 451 the track (S10.7, Shed 7, and 1163) rather than just the runout zone. Thus, some additional noise in the final 452 chronology for those specific paths could be due to more frequent small magnitude avalanches. Though the 453 oldest individual trees extended as far back as the mid-17 th century, the application of the double thresholds 454 processing steps restricted individual path avalanche chronology lengths since the minimum GD threshold 455 requirements were not met. It is difficult to place confidence in these older recorded events due to the 456 decreasing evidence back in time inherent in avalanche path tree-ring studies. Therefore, we chose to examine 457 more recent time periods dictated by the avalanche years identified through the double threshold methods. 458 All of the paths in the study are capable of producing large magnitude avalanches with path lengths greater 459 than 100 m (typical length for avalanche destructive size 2, D2), and all but RMC have a typical path length 460 of close to or greater than 1000 m (for avalanche destructive size 3, D3) (Greene et al., 2016). As Corona et 461 al. (2012) noted, the avalanche event must be large enough to create an impact on the tree, and size D2 or 462 greater will be evident from the tree-ring record (Reardon et al., 2008). However, the successive damage and 463 removal of trees from events size D2 or greater also impacts the future potential to record subsequent events 464 of similar magnitude. In other words, if a large magnitude avalanche removes a large swath of trees in one 465 year, then there are fewer trees available to record a slightly smaller magnitude avalanche in subsequent 466 years. Therefore, dendrochronology methods inherently underestimate avalanche events by up to 60% 467 (Corona et al., 2012), and our results suggest these methods captured about 10-50% of the available historical 468 record for JFS canyon. 469

Regional Sampling Strategy 470
By examining three different spatial scales (individual path, sub-region, and region) we produced a large 471 magnitude avalanche chronology for the region captured in a small subset of the total number of paths across 472 the large region. Accordingly, this sampling strategy may also alleviate the issue of recording large magnitude 473 avalanches within a region in the successive years following a major destructive avalanche event that 474 removed large number of trees within specific paths but not others. Overall, a regional sampling strategy 475 enables us to capture large magnitude avalanche events over a broad spatial extent that is useful for regional 476 avalanche forecasting operations and future climate association analysis. This strategy also allows us to 477 understand large magnitude avalanche activity at scales smaller than the regional scale. 478

Chronologies for Individual Paths and Sub-Regions 479
We applied the Wit threshold specifically to weight higher quality responses. The number of identified 480 avalanche years does not change for any individual avalanche path when we applied the Wit process. This 481 lack of change suggests that many of the responses in our samples were ranked as high-quality (i.e. C1-C2). 482 The high quality of responses can be attributed to the use of cross sections which allowed for a more complete 483 depiction and assessment of the tree-ring signal (Carrara, 1979). 484 We developed avalanche chronologies for 12 individual avalanche paths. The path with the greatest number 485 of identified avalanche years, S10.7, contains two major starting zones that are both steeper (35 and 39 486 degrees) than Shed 7, which also contains two separate starting zones. Reardon et al. (2008) collected a 487 substantial number of samples at higher elevations in the S10.7 avalanche path. However, the location data 488 for these samples were not available. Many of those samples were the living stumps that captured smaller 489 annual events. This is likely the root of the difference for S10.7 and the reason this path contains the largest 490 numbers of avalanche years in this analysis. 491 The range of return intervals across all paths (2 -28.5 years) is similar to those reported for 12 avalanche 492 paths across a smaller spatial extent in the Chic Choc Mountains of Quebec, Canada (2 -22.8) (Germain et 493 al., 2009). Although the authors in that study used a different avalanche response index, their study still 494 suggests considerable variation in avalanche frequency across avalanche paths within a region. 495 JGO contains the maximum return interval for any path in the study, and the return intervals are significantly 496 different than numerous other paths. A lack of recording data after one large avalanche event could easily 497 skew this value. To understand if this value is accurate, we would have to sample adjacent tracks to determine 498 if the return intervals are similar or not. Therefore, we cannot fully explain the large maximum return interval 499 for this path. However, one potential explanation is that this path is the only one located east of the 500 Continental Divide where the snowpack is often much shallower, particularly at lower elevations (Selkowitz 501 et al., 2002), thus inhibiting frequent large magnitude events from impacting the sampled runout zone. The 502 fetch upwind of this avalanche path is characterized by steep, rocky terrain harboring scoured slopes. This 503 topography limits the amount of snow available for transport to the JGO starting zone which may also 504 influence the load and stress placed upon this starting zone and subsequent large magnitude avalanches. 505 The return intervals for LJC in the Swan sub-region were the greatest in this sub-region and this is likely due 506 to wildfire activity in this path in 2003. LJC was heavily burned, and this created a steep slope with few trees 507 that was once moderately to heavily forested. Substantial anchoring and snowfall interception likely created 508 an avalanche path without many large magnitude avalanches for decades since slope forestation plays a 509 substantial role in runout distance and avalanche frequency in forested areas (Teich et al., 2012). In addition, 510 wildfires in 1910 burned a majority of the JFS sub-region as well, and the higher frequency of avalanche 511 years recorded between 1910 and 1940 in S10.7 suggests wildfire impacts may also be a contributor to the 512 high frequency of avalanche events in that location (Reardon et al., 2008). In addition, the fire in LJC may 513 also have removed evidence of previous avalanche activity. 514 Our results also suggest that return interval increases as path length increases, though the sample size for this 515 correlation analysis on individual paths is small (n=12). This result is likely because only very large 516 magnitude avalanches affect the far extent of the runout of the paths. This finding differs from a group of 517 avalanche paths in Rogers Pass, British Columbia, Canada, where path length was not significantly correlated 518 with avalanche frequency (Smith and McClung, 1997). However, that study used all observed avalanches, 519 including artillery-initiated avalanches, as opposed to a tree-ring reconstructed dataset. 520 The greatest number of identified avalanche years is in the JFS sub-region. The avalanche paths in this sub-521 region are all south or southeasterly facing whereas the other sub-regions span a greater range of aspects. 522 This narrow range of aspect may cause a bias toward overrepresentation of those aspects compared to the 523 inclusion of other aspects in other sub-regions. 524 The differences between individual avalanche paths as well as sub-regions are likely due to localized terrain 525 and weather/climate factors and the interaction of the two (Chesley-Preston, 2010). For example, Birkeland 526 (2001) demonstrated significant variability of slope stability across a small mountain range dependent upon 527 terrain and weather. Slope stability and subsequent large magnitude avalanching are likely to be highly 528 heterogeneous across not only the sub-region, but across a large region. This is also consistent with findings 529 by Schweizer et al. (2003) that suggest substantial differences in stability between sub-regions despite the 530 presence of widespread weak layers. Finally, climate drives weather, but is not a first order effect on 531 avalanche occurrence in any one given avalanche path. In this study, we derived a regional avalanche 532 chronology to provide a spatial scale that aligns more with the spatial scale of climate drivers than any one 533 individual path. These are relationships that should be examined in future work. 534

Regional Chronologies and RAAI 535
The regional chronology we developed through the use of tree-ring analysis on collections made across 12 536 avalanche paths suggests, unsurprisingly, that the inclusion of more avalanche paths across a large spatial 537 extent produces a more robust identification of major avalanche winters. When we aggregate all 12 paths 538 together and apply thresholds to discriminate the signal from the noise, we identified 30 avalanche years 539 throughout the region. This type of analysis allows us to place each avalanche year in context of the region, 540 or full extent of the scale triplet, rather than simply collating all major avalanche winters identified in each 541 individual path or sub-region. However, we also account for the support and spacing by including adjacent 542 avalanche paths within a sub-region and multiple sub-regions throughout the region. This sampling strategy 543 combined with the large sample size collected throughout the region allowed for a robust assessment of 544 regional avalanche chronology derived from tree-ring records. 545 We tested the sensitivity of the term regional by removing specific and random paths. Our results suggest 546 that removing paths from this structure, and subsequently compromising the sampling strategy, introduces 547 noise. By reducing the sample size, we reduce the ability of the thresholds to filter out noise, thereby 548 increasing the actual number of avalanche years in the region. However, the sample size reduction also 549 reduces the number of identified avalanche winters common to the full 12 path regional record (Table 6). 550 Our results emphasize the importance of sampling more paths spread throughout the region of interest as well 551 as a large dataset across the spatial extent. 552 Avalanche path selection is clearly important when trying to assess avalanche frequency (de Bouchard 553 d'Aubeterre et al., 2019). The importance of path selection is supported by our results suggesting that S10.7 554 is more influential than any other path in our study (Table 6) and is also illustrated by the large number of 555 avalanche years detected in S10.7 due to increased sampling in the track. However, by selecting multiple 556 paths representative of the range of geomorphic and potentially influential local weather controls throughout 557 the region we are able to provide a reasonable assessment of regional avalanche activity in areas without 558 historical records. By quantifying the sensitivity of the number of avalanche paths within a given region, we 559 illustrate that sampling a greater number of avalanche paths dramatically increases the probability of 560 identifying more avalanche years as well as increases the ability to reconstruct major widespread avalanche 561 events. However, as previously noted, dendrochronological techniques tend to underestimate avalanche 562 frequency, which implies that caution should be used when interpreting a regional avalanche signal using 563 this technique, particularly as sample numbers and qualities (e.g. cores vs. cross sections) decline. 564 Interestingly, the difference in median return interval throughout the "region" using 12 paths compared to 565 using four or eight paths changes only slightly suggesting that fewer paths are still able to represent the major 566 avalanche return intervals across a region. However, choosing fewer paths appears to introduce more noise 567 and therefore fewer years identified than a regional chronology with more avalanche paths. 568 The RAAI provides a measure of avalanche activity scaled to the number of active avalanche paths across 569 the region through time. The years with the greatest RAAI value coincide with substantial activity provided 570 in the historical record as well as previous dendrochronological studies from the JFS sub-region (Butler and 571 Malanson, 1985a, b; Reardon et al., 2008). The winter of 1932-33 was characterized by heavy snowfall and 572 persistent cold temperatures leading to extensive avalanche activity that destroyed roadway infrastructure in 573 the JFS sub-region, 1950 saw a nearly month-long closure of U.S. Highway 2 due to avalanche activity, and 574 in 2002, an avalanche caused a train derailment. While these are all confined to the JFS sub-region, with the 575 exception of 2002, they are also years shared by at least two other sub-regions. 576 We examined the probability of detecting an avalanche year throughout the region by sampling any one given 577 path. In seven of 30 years, the PODyear is only 8% and in all but three years the PODyear is less than 40%. The 578 low POD values are distributed throughout the time series, suggesting decreasing sample size back in time 579 or the number of active avalanche paths is not an influential factor. The POD is likely reflective of the spatial 580 variability of large magnitude avalanche occurrence across a region. It also aligns with the observational 581 findings of Mears (1992). During a major storm in 1986 throughout much of the western United States that 582 deposited 30-60 cm of snow water equivalent, Mears (1992) reported that in the area around Gothic, 583 Colorado, less than 40% of avalanche paths produced avalanches and less than 10% produced avalanches 584 approaching the 100-year return interval. This finding also supports the wide variability of avalanche years 585 between sub-regions recorded in our tree-ring record. Additionally, some of the avalanches in a given cycle 586 may not be large enough to be reflected in the tree ring record. Therefore, low values of PODpath when 587 considering only one avalanche path and identifying only one common year of large magnitude avalanche 588 activity (1982) amongst the sub-regions through dendrochronology is not surprising. Paths with at least one 589 GD (i.e. without applying thresholds) during avalanche years identified in the regional chronology exhibit a 590 greater PODpath, but this greater PODpath comes at the expense of introducing more noise if we were to simply 591 use one scar per path to define an avalanche event. 592 Our results also suggest that our sampling design using scale triplet increases the probability of detecting 593 avalanche activity across an entire region. We note that we are only able to scale our probability calculations 594 to our dataset with a limited historical observational record. However, our results illustrate the importance of 595 sampling more paths if the goal is to reconstruct a regional chronology. In our dataset, the greatest value of 596 PODpath is 40%, suggesting that if by chance, we sampled this path, we would have captured the regional 597 avalanche activity 40% of the time. 598 The trends in RAAI over the entire period of record are likely influenced by the decreasing number of samples 599 available to record an event further back in time. Despite the RAAI accounting for the number of avalanche 600 paths (minimum of n = 3), the small sample size from the late 19 th century precludes us from suggesting there 601 is a true increase in regional avalanche activity from 1867 to 2017. This is also supported by the absence of 602 positive or negative trend from 1950 to 2017 and 1990 to 2017. 603

Limitations 604
Overall, our results suggest that sampling one path, or multiple paths in one sub-region, is insufficient to 605 extrapolate avalanche activity beyond those paths. Multiple paths nested within sub-regions are necessary to 606 glean information regarding avalanche activity throughout those sub-regions as well as the overall region. 607 Our study is still limited by the underrepresentation inherent in dendrochronological techniques for 608 identifying all avalanche events. While we analyzed 673 samples over the extent of the region, some of the 609 paths in our study had relatively small sample sizes per individual path as compared to recent suggestions 610 (Corona et al., 2012). This may have influenced the number of avalanche years identified and subsequent 611 return intervals per individual path. However, we attempted to limit the influence of sample size by using 612 full cross-sections from trees, robust and critical identification of responses in the tree-rings, and appropriate 613 established threshold techniques. 614 We also recognize that sampling more avalanche paths in our region would certainly provide a more robust 615 regional avalanche chronology, but time, cost, and resource constraints required an optimized strategy. 616 Finally, our study would undoubtedly have benefited from a longer and more accurate historical record for 617 comparisons and verification of the tree-ring record in all of the sub-regions. Overall, our study illustrates 618 the importance of considering spatial scale and extent when designing, and making inferences from, regional 619 avalanche studies using tree-ring records. 620

Conclusions 621
We developed a large magnitude avalanche chronology using dendrochronological techniques for a region 622 in the northern U.S. Rocky Mountains. Implementing a strategic sampling design allowed us to examine 623 avalanche activity through time in single avalanche paths, four sub-regions, and throughout the region. By 624 analyzing 673 samples from 12 avalanche paths, we identified 30 years with large magnitude events across 625 the region and a median return interval of ~3 years (from 1866-2017). Large magnitude avalanche return 626 interval and number of avalanche years vary throughout the sub-regions, suggesting the importance of local 627 terrain and weather factors. Our work emphasizes the importance of sample size, scale, and spatial extent 628 when attempting to derive a regional large magnitude avalanche chronology from tree-ring records. In our 629 dataset, the greatest value of PODpath is 40%, suggesting that if we sampled only this path, we would have 630 captured the regional avalanche activity 40% of the time. This clearly demonstrates that a single path cannot 631 provide a reliable regional avalanche chronology. Specifically, our results emphasize the importance of 1) 632 sampling more paths spread throughout the region of interest; 2) collecting a large number of cross-sections 633 relative to cores; and, 3) generating a large dataset that scales to the appropriate spatial extent. Future work 634 should include conducting a similar study with a number of paths in the same sub-regions for verification, or 635 in an area with a more robust regional historical record for verification. 636 637 6. Appendix A 638 Table A1: List of previous avalanche-dendrochronological work with more than one avalanche path in study -to 639 place our regional work in context with other regional/multiple path studies. Number of samples, paths, growth 640 disturbances (GD), and spatial extent (linear distance between most distant avalanche paths in study area) are 641 included. For spatial extent, NA is reported in studies where spatial extent is not reported or could not be inferred 642 from maps in the published work. Where spatial extent is not reported directly in previous work, it is estimated 643 by using maps from the published work and satellite imagery. We included only the initial studies using a dataset 644 with more than one avalanche path. For example, if a study used the same dataset again in subsequent work, we 645 did not include it.    EP responsible for study conception and design, data collection, analysis, writing. JH contributed to 663 development of study design, methods, editing, and writing. DS responsible for data collection, tree-ring 664 processing and analysis and writing. GP, KB, and DF contributed to study design, editing and writing. 665

Disclaimer 666
Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by 667 the U.S. Government. 668