Snow Avalanche Frequency Estimation (SAFE): 32 years of monitoring remote avalanche depositional zones in high mountains of Afghanistan

. Snow avalanches are the predominant hazards in winter in high elevation mountains. They cause damage 9 to both humans and assets but cannot be accurately predicted. Here we show how remote sensing can accurately 10 inventory large avalanches depositional zones every year in a large basin using a 32-y snow index derived from 11 Landsat satellite archives. This Snow Avalanche Frequency Estimation (SAFE) built in an open-access Google 12 Engine script maps snow hazard frequency and targets vulnerable areas in remote regions of Afghanistan, one of 13 the most data-limited areas worldwide. SAFE correctly detected of the actual depositional zones of avalanches 14 identified on Google Earth and in the field (Probability of Detection 0.77 and Positive Predictive Value 0.96). A 15 total of 810,000 large depositional zones of avalanches occurred since 1990 within an area of 28,500 km 2 with a 16 mean frequency of 0.88 avalanches km²y-1 , damaging villages and blocking roads and streams. Snow avalanche 17 frequency did not significantly change with time, but a northeast shift of these hazards was evident. SAFE is the 18 first robust model that can be used worldwide and is especially capable of filling data voids on snow avalanche 19 impacts in inaccessible regions.


Introduction 22
Snow avalanches are among the fastest, up to 61 m s -1 , and therefore most dangerous natural hazards in mountain 23 areas (Louge et al., 2012). Casualties associated with avalanches are numerous; in 2021 alone, 37 fatalities occurred 24 in the US (Colorado Avalanche Information Center, 2021) and 127 in Europe (European Avalanches Warning 25 Services, 2021), but avalanche monitoring is not consistent across the globe. Most remote mountain regions and 26 communities are not systematically monitored for avalanche occurrence. Avalanche surveys amongst remote 27 villages are sparse because regions are uninhabited; however, avalanches can block connecting roads every year 28 since avalanche volumes range from hundreds to several tens of thousand cubic meters (Gubler, 1987). Where 29 weather stations exist, avalanches can be predicted based on snow depth and other weather parameters (Greene et 30 al., 2016). However, the global weather monitoring of mountainous areas is scattered and very sparse in developing 31 nations. 32 detection. However, the acquisition of frequent radar images is too recent to use this technique to detect historical 66 avalanches. In addition to optical, radar or Lidar data, other studies used Digital Elevation Models (DEMs) and 67 topographic parameters to determine the influence of terrain on avalanches in Switzerland (Maggioni and Gruber, 68 2009). Other studies incorporated other parameters such as morphology and vegetation to define potential 69 avalanche zones and ran the Avalanche Flow and Run-out Algorithm to automatically detect potentially affected 70 regions by avalanches (Barbolini et al., 2011). Moreover, the combination of snow measurements (depth) and high 71 resolution DEMs have proved useful in snow hazard detection (Bühler et al., 2018a(Bühler et al., , 2022. Lidar is being used in 72 the same regard with a higher level of precision. Lidar sensors measure snow depth before and after events at 73 submeter resolutions (Prokop, 2008;Deems et al., 2013;Prokop et al., 2013;Hammond et al., 2018). However, 74 this technology remains very expensive, and the spatial coverage is limited. Therefore, Lidar data are not suitable 75 for avalanche detection at the basin scale. In addition to optical, radar or Lidar data, other studies used Digital 76 Elevation Models and topographic parameters to determine the terrain influence on avalanches in Switzerland 77 (Maggioni and Gruber, 2009). Some other studies added other parameters such as morphology and vegetation to 78 define potential avalanche zones and ran Avalanche Flow and Run-out Algorithm to automatically detect regions 79 that are potentially affected by avalanches (Barbolini et al., 2011). Moreover, the combination of snow 80 measurements (depth) and high resolution DEMs was shown to be efficient for snow hazard detection (Bühler et 81 al., 2018a). 82 83 Optical data are the most available data in terms of spatial and temporal resolution as well as historical archives. 84 Thus, we used optical data to detect avalanches on a long-term basis. Landsat-5, 7 and 8 products were used as 85 their resolution (30 m, 900 m²) is sufficient to detect small avalanches . Most of these 86 data are available at a global scale. Optical sensors can detect areas covered or not covered by snow and this 87 approach has been used in multiple studies during the past decade. Manual approaches or indices have been used 88 in such studies. For example, Landsat-8 Panchromatic images (15 m) in combination with radar images were used 89 to detect avalanches in Norway (Eckerstorfer et al., 2014). Such combinations were also recently used in west 90 Greenland to map a large number of avalanches after an unprecedented snow event (Abermann et al., 2019). To 91 our knowledge, only one recent study automated the detection of avalanches using remote sensing products and 92 an open-access scripting approach (Smith et al., 2020). This study downloaded avalanches annually for a given 93 region of interest using available Landast-8 images and computed NDSI for each image. NDSI differentiated so-94 called 'supraglacial debris' from snow cover, for the date of interest. However, this approach only covers high 95 elevation areas while our study aims to detect avalanches proximate to local communities at lower elevations 96 (typically valleys). Manual and visual approaches, despite the time consuming process, can also be applied to 97 detect avalanches using high resolution images (e.g., SPOT-6), mid-resolution (e.g., Sentinel-2A 1989 1991 1993 1995 1997 1999 2001 2003 2005 2007 2009 2011 2013 2015 2017 2019 2021

Observed years
Landsat-5 Landsat-7 Landsat-8 As the aim of this study is to detect and map the annual occurrence of depositional zones during the past 32 years 171 within the study area, the monitoring approach must be reasonable and transferable from year to year. Based on 172 frequent field observations and literature , the authors noticed that depositional zones 173 can be detected using the contrast between snow cover and bare cover, but the timing is perhaps the most important 174 consideration. Indeed, the script is based on the assumption that snow packages exist in lowlands, especially along 175 rivers and streams, as late as May through mid-July. At this time of the year, the terrestrial snow cover has largely 176 melted and only snow packages triggered by avalanches remain. The location of those snow packages is also very 177 critical (i.e., along riverbanks). These zones are indeed detectable by delineating the depositional zones of the 178 avalanches (not their release or transition zones); in most cases these were located on river or stream banks as 179 observed in the field because the hillslopes always route snow avalanches in this direction. We cannot differentiate 180 between dry, wet, or powder snow because the process detects the remaining snow packages as avalanches in the 181 late season (spring and summer), not in winter, nor can we delineate multiple deposits within the same depositional 182 feature, only the combined deposit zones. In winter, we were not able to differentiate contrasts between snow cover 183 and avalanches, thus our focus is on the late season. 184 185

Google Engine interface and code availability 186
The concept of detecting the 'remaining snow avalanches deposits in the late season' was written in Java Script 187 using the Google Engine platform. The script SAFE is available at 188 https://code.earthengine.google.com/2a3bc50702a76d6989119a4a6c07af80. We selected Google Engine for its 189 relative simplicity of use and open access code, which is available to all stakeholders involved in hazard and 190 vulnerability assessments. Additionally, internet connections in remote areas, such as within the Amu Panj basin, 191 are limited and powerful computers required to run scripts and process big data are sparse. Our script can be run 192 by anyone in a reasonable amount of time, even with a low internet capacity. As an example, yearly depositional 193 zones of avalanches in our study area were downloaded and mapped from Badakhshan (SAFE was processed from 194 Khorog, University of Central Asia campus, in Tajikistan A 200 m riparian buffer was used as a mask to clip the Landsat images. Because our area of coverage encompasses 211 very different elevations, the date of snow melt is not uniform throughout the basin. Therefore, distinguishing 212 between the depositional zone and bare land requires different times depending on elevation. To accomplish this, 213 we calculated the average elevation of the snowline for the last 20 years using MODIS products. To distinguish 214 the different melt timing between highlands from lower areas, we selected the summer snowline (June-July-215 August; JJA). The average elevation of the JJA snowline was 4420 m during the past 20 years. Two masks were 216 therefore produced: one with a river buffer in lowlands and another for highlands. Those masks are only relevant 217 if the user carefully selects the date of interest. For lowlands (below 4420 m), our time window was 15 May to 15 218 June, indicating that the script downloads and compiles all available Landsat images acquired in this range and 219 detects the deposit zones efficiently because during that period the terrestrial snow cover has already melted and 220 the deposits are easily recognised. For higher elevations (above 4420 m), snow cover melted later; dates to 221 accurately distinguish the remaining snow packages ranged from 15 June to 15 July. After many tests, it was 222 confirmed that these date ranges reproduced the desired snow conditions during the entire 32-y period. In the 223 script, users can modify these dates (lines 24 and 112) to conform to local conditions. 224 225

Snow index reclassification 226
After the construction of the mask, SAFE proceeds as outlined in Figure 3. NDSI is selected to detect snow of 227 deposit zones in the script for its transferability from one Landsat generation to another. NDSI computes a ratio 228 between VIS and SWIR bands of Landsat satellites with negative NDSI representing non-snow cover and positive 229 values indicating snow coverage (Equation 1). Three cover types were distinguished to detect depositional zones 230 of avalanches at the correct time: (1) bare lands; (2) water bodies; and (3) snow. The values in Table 1 were 231 established after multiple tests before obtaining sufficient precision to distinguish deposit zones from other land 232 covers. On each mosaic (composite of the available images during the period of interest), a cloud mask is applied 233 using Landsat QA bands in the script to remove clouds from the scene. 234 Equation 1 235 236 designates snow avalanche deposit zones is selected in the script. From the selected reclassification, the script 241 removes the standalone pixels because their classification might not be precise or representative of actual cover. 242 Next, the selected 'avalanche pixels' are verified into the script avoiding manual vectorization after the 243 downloading process. The vectorization procedure of depositional zones of avalanches is justified by the analysis 244 steps and post-processing after downloading data. Depositional zones of avalanche statistics, elevations, and 245 surface areas are extracted from vector files. Finally, annual avalanche deposit zone shapefiles are exported into 246 the Google Drive user's account. 247 248

Figure 3. Flow chart for Snow Avalanche Frequency Estimation (SAFE) using Landsat archives in Google 249
Engine. 250 251

Depositional zones of avalanche surface area classification 252
Once the data are downloaded and imported into the GIS environment, statistical analysis commences. Every year, 253 the number and areas of deposit zones are calculated to quantify their evolution. Moreover, the surface areas of 254 the depositional zones are classified. Although a generic surface area classification exists (Greene et al., 2016), 255 we decided to classify avalanches by surface areas based on local conditions. We segregated four discrete 256 categories of deposit zones: small (< 1000 m²); medium (1000-5000 m 2 ), large (5000-15,000 m 2 ), and very large 257 (15,000-100,000 m 2 ). Such a classification enabled us to assess the intensity and potential impact of those hazards 258 in specific locations. SAFE is not able to detect the avalanches at their time occurrence, and since these hazards 259 are detected weeks after initiation, their surface area is underestimated by SAFE due to melting. However, the 260 estimated surface areas in SAFE are still useful for classifying depositional zones of avalanches by surface area 261 since large snow deposits melt slower than small snow deposits. The small avalanches that occurred in winter will 262 appear as small deposit at the time of extraction in SAFE and the large events as large hazards since visible snow 263 deposits can be seen in late spring. 264 265 3. Results 266

Validation 267
The performance of SAFE in correctly detecting snow avalanche depositional zones required careful assessment. 268 To achieve this, we collected datasets that show actual locations ( The results suggest a good reliability of SAFE (Table 2). The overall POD is 0.77 which means that SAFE 294 identified a significant number of the depositional zones of avalanches that impacted valley bottoms. Moreover, it 295 seems that SAFE performs better in detecting true positive deposit zones (that occurred on the ground), as shown 296 by the high PPV scores (average: 0.96). SAFE almost never detected depositional zones of avalanches that did not 297 exist. However, SAFE might miss some deposit zones due to cloud cover on the Landsat images, especially in 298 Another source of error arises when SAFE cannot detect depositional zones due to a dark color on the snow surface 303 associated with surface debris or a debris flow on top of the deposit zones. NDSI may have identified those debris 304 layers as bare soil in the classification. Moreover, it should be understood by the users that another limitation is 305 that SAFE does not detect early winter avalanche deposits due to melting and snow coverage on and around the 306 snow deposit, which might affect the deposits frequency estimations. However, based on our findings, SAFE can 307 be considered as a conservative, yet robust and efficient tool to automatically identify snow avalanche depositional 308 zones in very remote areas and can be applied in any mountainous region. 309 310 3.2 SAFE outputs compared with outlined avalanches using SPOT-6 images 311 As a potential method of strengthening our testing of SAFE, outputs of our model were compared with a method 312 that applied a more precise and expensive remote sensing product in Switzerland in 2018 (Bühler et al., 2019b;313 Hafner and Bühler, 2018). The Swiss area encompassed 12,500 km² where more than 18,000 snow avalanches 314 were manually digitized using very high-resolution SPOT-6 images (in January 2018). While our dataset is quite 315 different from the Swiss data, the objective of this comparison was to assess how many snow avalanche deposits 316 SAFE could detect compared to the approach using SPOT-6 (Table 3). Figure 4 shows an illustration of this 317 comparison. It appears that the deposit zones detected by SAFE are in line with SPOT6 outlined avalanches. The 318 later however covers the entire avalanches while SAFE only detects, automatically, the deposit zones. 319   (Bühler et al., 2019b;Hafner and Bühler, 2018) 327 328 Importantly, not all avalanches manually digitized on SPOT-6 images were comparable to SAFE results. To make 329 this comparison more consistent, we clipped the outlined avalanches with the valley bottom mask used in SAFE. 330 Following this modification, the SPOT-6 digitization process identified 7574 avalanches deposits in valley bottoms 331 compared with 9948 by SAFE. Overlapping these two datasets, we found that both approaches detected 2194 332 deposit zones in common. Much of this discrepancy is due to the timing of SAFE images, which examine deposits 333 that remain in late spring and early summer, whereas SPOT digitization covered only January. The larger number 334 of snow deposits detected by SAFE occur during late season snow avalanches that impact valleys. This suggests 335 that SAFE could not detect all January snow deposits because many of those already melted by the time of SAFE 336 detection (early April to late June in the Swiss case). In addition, optical image quality strongly depends on cloud 337 cover that may cause avalanches to be obstructed. For instance, we could not compare the 2019 SPOT-6 derived 338 dataset in eastern Switzerland  due to cloudy images at the end of winter and early spring 339 because these snow avalanches had already melted, implying that SAFE is more suitable for high mountain areas 340 (>4000 m) where snow deposits remain longer in valleys, thus inflicting greater damages and obstructions. Using 341 LANDSAT images, SAFE somewhat circumvents this problem of cloud cover by assessing many years of data 342 (in our case 32 years). However, SAFE does not distinguish individual events and considers overlapping snow 343 deposits as one, in contrast to SPOT-6 which distinguishes these as discrete events. This, in addition to the different 344 methods and spatial resolution difference between SAFE and SPOT, explains the somewhat low number of 345 overlapping snow deposits between SAFE and SPOT. Moreover, the SPOT digitization procedure found a total 346 avalanche area of 362,187,741 m² in January, while SAFE detected 494,454,599 m² of deposits at the end of the 347 avalanche season, including 223,907,868 m² in common. The area detected by SAFE is naturally larger than SPOT-348 6 since SAFE maps all detectable deposits at the end of the winter. Moreover, SAFE did not detect the small 349 avalanches of January that rapidly melted after they occurred. The polygons extracted by SAFE using Landsat 350 images are obviously coarser than those outlined with SPOT-6 images, which partly explains the low number of 351 overlapping snow deposit zones, but a much more comparable detected area (62%) between the two methods. 352 Much of the discrepancy is related to SAFE's inability to detect individual events and missing deposits that rapidly 353 melt (mostly from the early winter snow avalanches), as well as the very different resolution of these products. 354 355

Snow avalanche depositional zone frequency from 1990 to 2021 356
By compiling 32 years of satellite images (see Methodology), the frequency of avalanche deposits at a 900 m² 357 pixel scale was determined ( Figure 5 and 6a). SAFE inventories snow avalanche deposits that occurred within a 358 year and therefore identifies the most vulnerable areas, but it does not aim to forecast the timing of future 359 avalanches. During this period, some 810,000 depositional zones impacted valleys within the Amu Panj basin 360 (28,500 km²), i.e., approximately 28 depositional zones km -2 . Each year these avalanche deposits cover an average 361 of 1.23% of the basin area but surface area varies from year to year. Avalanche depositional zone surface area 362 ranged from 900 to 100,000 m 2 and is categorized into four classes: small (< 1000 m²); medium (1000-5000 m 2 ), 363 large (5000-15,000 m 2 ), and very large (15,000-100,000 m 2 ). The most frequent are medium-size surface area 364 deposit zones; 342,000 events during the past 32 years. Our approach also identifies very large snow avalanche 365 depositional zones that pose the greatest danger to local populations and infrastructure. We found no correlation 366 between altitude of depositional zones and their surface areas. Avalanches deposits in this region have an average 367 altitude of 3820 m and the lowest depositional zone occurred at 1755 m. 368 369 These spatial and temporal statistics allow for a geographic assessment of the avalanche deposits. In total, ten sub-370 catchments (ranging from 18 to 240 km 2 ) were impacted by more than one avalanche depositional zone km -2 y -1 , 371 with an average frequency of 0.26 deposit zones km -2 y -1 throughout the Panj Amu basin (Figure 7) Our remote sensing approach facilitates innovation in snow avalanche depositional zone monitoring: i.e., detecting 383 avalanches deposits outside of populated areas, especially along roads that are frequently blocked by avalanches 384 (Figure 9). More than 2000 roads in the basin (5.47% of the road network) were affected by avalanche deposits 385 every year. Additionally, more than 400 roads in Upper Badakhshan and Wakhan regions experienced more than 386 2 avalanche deposits y -1 km -1 of road (within a 1 km buffer). The average frequency along roads is 0.86 avalanche 387 deposits km -1 y -1 during the past 32 years, most of these in the medium-size surface area category. 388

Stream blocking and resultant flooding
Damages to infrastructure and blocking of roads by avalanche deposits are not the only consequences of these mountain 415 hazards. Because depositional zones typically reach rivers in this steep, incised terrain, the sudden and rapid arrival of several tons of snow can temporarily block rivers inducing short-term localized flooding. By cross-checking the map of the rivers in the Amu Panj basin with SAFE outputs, it appears that 26.2% of the river network is impacted by avalanche deposits, mainly in the high mountains. During the past 32 years, 12% of the streams have been blocked by at least 10 avalanches deposits km -1 representing a significant risk for villages and farms in floodplains. The accumulated snow mass impounds river water until 420 it can break through releasing a large discharge surge. Thus, depending on the surface area of the avalanche deposit with respect to the channel dimensions, damages to villages and farmlands may occur both upstream due to impounded water (hours to weeks) and to downstream following the sudden release of water.

Snow avalanche depositional zone trends during the past 32 years 425
This long-term monitoring of snow avalanche deposits facilitates the assessment of the evolution of these rapid mass movements. During the 32 years of avalanche depositional zone assessment, no significant temporal trends in impacted areas were detected (Figure 10). In addition, there was no significant trend of the surface are of snow deposits (p-value > 0.05).
Nevertheless, some years posed much greater risk than others. In the last 32 years, ten years have been more at risk with aboveaverage avalanche deposit coverage: 1990coverage: , 1991coverage: , 1992coverage: , 1993coverage: , 1994coverage: , 1995coverage: , 1996coverage: , 2003coverage: , 2005coverage: , 2007coverage: and 2012 2003 had many avalanche depositional zones that occupied almost 6% of the surface area of the entire basin. That year was locally noted as having heavy snowfall and farmers benefited from more snowmelt in the spring, leading to higher than average crop yields in 2003 (FAO, 2003;Guimbert, 2004). Notably, the higher risk years were also characterised by lower altitudes for avalanche deposits. There is a slight negative correlation (-0.55, Pearson test) between altitude and total annual avalanche area.
With larger avalanche areas, deposits reach closer to villages. For example, in 2003, the lowest avalanche depositional zone 435 occurred at an altitude of only 1871 m, very close to housing clusters and roads. It is therefore possible that communities below 2000 m are also impacted by snow avalanche deposits and in many mountain regions of the world this represents a significant proportion of the communities living proximate to these altitudes.

Temporal geographic shifts of snow avalanche deposit zones
Long-term monitoring also shows the evolution of the spatial distribution of snow avalanche depositional zones. The pattern 445 of snow avalanche deposits has changed with time and slightly shifted to the northeast portion of the basin; thus, more avalanches are now occurring in the northeast than in the southwest (Figure 11a and b). Nevertheless, snow coverage did not shift simultaneously according to our remote sensing analysis nor did the snowline evolve, but rather remained variable over the last 32 years. The geographic shift of avalanche depositional zones is therefore likely due snow depth evolution. Deeper snowpacks trigger snow avalanches. There are no available data on snow depth at such a scale. However, the slope of the 450 temperature was calculated and a Mann-Kendal test was applied for each pixel of the land surface temperature images (MOD11C3). Remotely sensed land surface temperature changed during the last 20 years (Figure 11 a), with a warmer band occurring through the central portion of the basin in December (p-value 0.03 with an increase of 0.88 C°y -1 ). This central portion is mainly mountainous, and this temperature pattern may have shifted the avalanches to the northern mountains of the area, while the south is characterized by lower mountains. Overall, avalanche depositional zone locations tend to follow the 455 spatial distribution of snow depth . This means that despite the high variability of the snow line and snow coverage, the distribution of snow avalanche deposits can significantly change over time in response to temperature changes and local communities must be prepared for shifting hazards.

Sensitivity analysis of SAFE
To better understand how SAFE works and assess its performance, a sensitivity analysis was conducted between the model parameters. The number and surface areas of avalanche depositional zones vary according to the buffer used, the dates of Landsat images, and finally the NDSI range during the snow classification. The sensitivity analysis was conducted for the year 2019 when SAFE was most robust in valleys where actual avalanche deposits were quite visible on Google Earth images 470 (POD: 0.84 and PPV: 0.94). First, we run SAFE with different buffer widths (25 m of difference between each buffer). There is a strong positive correlation (0.98) between the number of avalanche deposits detected by SAFE and the buffer width ( Figure   12a). The wider the buffer around the rivers, the more avalanche deposits SAFE will detect. On the other hand, for narrower buffers, the average surface area of avalanche depositional zones is smaller (positive correlation of 0.71). This is because a large buffer extends upslope where small snow patches reside, which are not avalanches deposits since they are located at the 475 top of hillslopes. This means that the user should not select a buffer that is too wide, rather the area should only include the riparian zone of rivers and streams where the snow avalanche deposits are located. As such, we used a value of 200 m for the entire region studied. The number and surface areas of avalanche depositional zones detected by SAFE depends on the NDSI range when classifying 485 snow. NDSI is used to differentiate between water bodies, bare lands, and snow. By varying the NDSI ranges of snow in the script, we notice a strong positive correlation with the number of avalanche deposits detected by SAFE. The closer the index is to 0, the more hazards SAFE finds. However, this correlation shows that the choice of NDSI range is important because we notice a threshold at 0.31 (Figure 12b). Avalanche depositional zones seem to be more numerous with an NDSI lower than 0.31 because the snow pixels are confused with water bodies. It is therefore essential for the user to select an NDSI higher than 490 0.31 to distinguish water bodies (rivers, flood areas or lakes) and snow. However, there is no correlation between the NDSI ranges and the average surface areas of avalanche deposits because NDSI cannot interpret pixels other than 'snow' above the 0.31 threshold. Finally, the date of interest is a key parameter in SAFE. The number of avalanche depositional zones detected by SAFE is highest at the end of winter due to the almost constant cloud cover since January, but also due to the inability to distinguish avalanche deposits from snow cover in winter (with Landsat images) (Figure 12c). May is a key month for SAFE 495 applications in high mountains: the snow coverage, which is thinner than the avalanche depositional zones, begins to melt and the number of avalanche deposits detected can then be assessed. It is therefore essential to select post-May images to detect avalanche deposits while taking care not to select post-July images as avalanches deposits melt in summer precluding detection.
Some avalanche deposits are also visible on successive images after the snow cover melts. SAFE was specifically designed to 500 detect avalanche depositional zones at their earliest stage after snow melts. Indeed, starting from May (when the avalanche deposits are not confused with snow coverage), snow avalanche deposits start to melt and the surface area will begin to be underestimated. For that reason, it is important to select late spring images for lowland avalanche depositional zones and early summer for highland deposits, not later. Cloud cover is another issue in avalanche deposit locations and surface area detection since cloud cover may partly or fully obstruct the avalanche deposits at the time of the image. This is another reason to select 505 images starting from late spring when regional cloud cover is lowest and even absent in early summer. If cloud cover is high even late spring, the users can still select later images, but there will be a risk that detected avalanche deposits will have started to melt. To summarize, we recommend the following three parameters in the SAFE script: buffer of 200 m to include only snow avalanche deposits; NDSI > 0.31 to distinguish water bodies from snow, and images from May to July to distinguish avalanche depositional zones from snow cover. 510

Excluding snow coverage
Interpreting the remaining snow packages as avalanche deposits can lead to some errors. Indeed, despite a precise masking operation (excluding summits and very high plateaus where snow persists), in some cases the use of NDSI might not properly segregate avalanche depositional zones from large areas of remaining snow. After assessing the surface area of true avalanche 515 deposits (the ones that SAFE correctly detected based on Google Earth images), it appeared that snow cover areas > 100,000 m² were not avalanche depositional zones but rather only snow cover; thus, these were removed. However, in highlands, even along riverbanks, some snow packages interpreted as avalanche depositional zones may be remaining snow cover. As such, the date range for highlands was selected as late as possible in the year. Thus, it is advised to keep the mask at the very bottom of valleys (maximum 200 m buffer along the river) to exclude high plateaus and potential snow-covered areas. 520

Water bodies in SAFE
A final limitation of using this remote sensing and NDSI approach for avalanche deposit detection is the possible confusion between some small water bodies and avalanche depositional zones. Indeed, in some cases certain river reaches (stream order > 4 in our study area) could be interpreted as snow because they were frozen and appeared as white pixels on Landsat archives. 525 The same issue can occur with ponds and lakes. This limitation was foreseen before processing the images in our study and we excluded these large water bodies from the region of interest (in the mask) by using available shapefiles. For example, Shiva Lake, one of the largest water bodies in Amu Panj basin (15 km²), was removed from the analysis. Another way to avoid the water pixel selection is to adapt the NDSI reclassification itself, depending on the study area. This is possible in the script lines 51-53 for low elevations and lines 139-141 for high elevations in the script. 530

SAFE outcomes compared to other snow avalanches detection studies
SAFE contributes to the literature on snow avalanche detection, but in a unique way using remote sensing. As noted, many studies and models exist using various products: radar, optical, and topographic. The strength of remote sensing is the automatic processing across wide scales and over long timeframes. SAFE uses the capabilities of remote sensing by processing more 535 than one image per year at the catchment scale. Moreover, the use of Landsat archives allows assessment over the last 32 years, which is not yet possible with recent radar data such as Sentinel-1. Most of the current avalanche detection models use freely available products, with acceptable if not good accuracy (Table 4). The accuracy of these studies using radar images ranges from 53 to 81% making this a relatively robust tool. One of the reasons why SAFE does not use radar images is the weight of the images (data storage), especially Sentinel-1, which is mostly above 1 Gb/image. These heavy images are not suitable for a 540 model like SAFE, which was specifically designed for remote study areas where internet connections may be very limited.
Other models also exist with optical images with high accuracy ranging from 71 to 93% (Table 4). In the optical domain, SAFE showed a POD of 77% over an area of 28,500 km². SAFE is therefore in the high range of models with optical, medium resolution (Landsat) images.  Bühler et al., 2018 95** DTM *55 avalanches were detected using S1 image out of 102 on the field. **POD 555

Conclusion
SAFE can be considered as a universal approach to assess snow avalanche depositional zones in spring and early summer where ground data are very limited, such as in the Afghan mountains. Here we showed the capability of long-term remote sensing data to robustly detect snow avalanche deposits that impact valley locations. While we have successively applied 560 SAFE to assess the frequency and impacts of avalanche deposits in valleys and lower hillslopes of Afghanistan, arguably one of the most data-limited regions worldwide, this model should perform even better in areas where snow data are available making it an important tool for avalanche vulnerability assessment worldwide. More than 30 years after the launch of Landsat-5, it is now possible to compile all data and assess the temporal as well as spatial evolution of such hazards. NDSI is a relevant index to detect avalanches when selecting the correct region and dates of interest -i.e., riverbanks during the late melt season. 565 The thickness of the depositional zones facilitates the detection of these avalanche deposits after the snow cover has melted on hillslopes in spring or early summer. Moreover, the application of SAFE in Afghanistan, compared to its application in Switzerland, showed that the script can be applied worldwide, especially in high mountains (above 4000 m) since deposit zones are still detectable in late spring at those elevations.

570
The automation of snow avalanche detection using remote sensing technologies at regional scales is still new and SAFE was designed to guide decision-makers, planners, and disaster risk practitioners. Indeed, such people can now know where the most at-risk areas are located based on these frequency maps. Such information informs the relative risk of building sites and land use decisions in such mountainous terrain with greater precision. The level of exposure of roads to avalanche depositional zones can also be estimated using these frequency maps and can inform road planners and managers regarding road location, 575 maintenance practices, and mitigation structures. Moreover, villages of high mountains such as in Afghanistan are highly dependent on road access and connections to provide necessary food, energy, medical supplies, and life-support, especially in winter. It is therefore critical for local decision makers to assess the frequency of road blockage by avalanche deposits. Thus, open-access and user-friendly tools such as SAFE are highly applicable to interests of local stakeholders even with medium to low power computers since SAFE uses Google servers. The tourism sector can also benefit from this snow avalanche deposit 580 inventory, especially the winter sports industry. Furthermore, this method can also be used to prioritize areas for more sophisticated and data-intensive avalanche risk analysis (Keylock et al., 1999). SAFE can be applied by any user throughout mountainous regions of the world as it is designed to be user-friendly, and frequent users can contribute to the robustness of the snow avalanche deposit archive, thus improving recommendations for policy makers.