Automated snow avalanche release area delineation in data sparse, 1 remote, and forested regions 2

. Potential avalanche release area (PRA) modelling is critical for generating automated avalanche terrain maps which 9 provide low-cost large scale spatial representations of snow avalanche hazard for both infrastructure planning and recreational 10 applications. Current methods are not applicable in mountainous terrain where high-resolution (≤ 5 m) elevation models are 11 unavailable and do not include an efficient method to account for avalanche release in forested terrain. This research focuses 12 on expanding an existing PRA model to better incorporate forested terrain using satellite imagery and presents a novel approach 13 for validating the model using local expertise, thereby broadening its application to numerous mountain ranges worldwide. 14 The study area of this research is a remote portion of the Columbia Mountains in southeastern British Columbia, Canada which 15 has no pre-existing high-resolution spatial data sets. Our research documents an open source workflow to generate high-16 resolution DEM and forest land cover data sets using optical satellite data processing. We validate the PRA model by collecting 17 a polygon dataset of observed potential release areas from local guides, using a method which accounts for the uncertainty of 18 human recollection and variability of avalanche release. The validation dataset allows us to perform a quantitative analysis of 19 the PRA model accuracy and optimize the PRA model input parameters to the snowpack and terrain characteristics of our 20 study area. Compared to the original PRA model our implementation of forested terrain and local optimization improved the 21 percentage of validation polygons accurately modelled by 11.7 percentage points and reduced the number of validation 22 polygons that were underestimated by 14.8 percentage points. Our methods demonstrate substantial improvement in the 23 performance of the PRA model in forested terrain and provide means to generate the requisite input datasets and validation 24 data to apply and evaluate the PRA model in vastly more mountainous


Introduction
Snow avalanches are a significant natural hazard for traffic and settlement infrastructure as well as for individuals who travel in snow covered mountainous regions.Roads, railroads, utilities, and permanent structures that are located in areas with potential avalanche hazard can be destroyed by large avalanche impacts or blocked for extended periods during winter storm events causing financial losses and potential for injury or death from individuals being buried in the debris.In economically developed countries, the majority of avalanche fatalities occur during recreational activities (i.e., backcountry skiing, snowmobile riding, mountaineering) where individuals voluntarily expose themselves to avalanche hazard (Boyd et al., 2009), and accidental avalanches are mostly triggered by the party that is caught (Schweizer and Lütschg, 2001;Techel et al., 2016).
In North America and Europe, an average of approximately 140 people are killed in avalanches each year (Jamieson et al., 2010;Techel et al., 2016;Colorado Avalanche Information Center, 2020).
To mitigate avalanche hazard, locations with potential for avalanche release need to be identified so elements at risk can attempt to minimize their exposure.This can be achieved by avoiding avalanche prone areas, minimizing exposure time, or implementing avalanche control methods (McClung and Schaerer, 2006).Avalanche hazard mapping is a time honored practice for determining the spatial distribution of snow avalanche hazards (Margreth and Funk, 1999;Rudolf-Miklau et al., 2015).Traditional manual hazard mapping combines multiple methods such as terrain inspection, numerical simulations, avalanche event databases and expert experience to evaluate avalanche hazard exposure and spatial extent making it both labor and cost intensive.This highly detailed approach is the gold standard for determining avalanche zoning for permanent infrastructure, but the costs make it unsuitable for mapping large areas of mountainous terrain (Rudolf-Miklau et al., 2015, Bühler et al., 2018, 2022).
To overcome this challenge, automated GIS and remote sensing based methods have been developed to expedite the mapping process and produce avalanche terrain indication maps based on digital elevation model (DEM) and land cover data (Maggioni and Gruber, 2003;Gruber and Haefner, 1995).The foundation of automated avalanche terrain mapping is potential avalanche release area (PRA) modelling, which estimates the location of potential hazards based on the local terrain characteristics (Bühler et al., 2013(Bühler et al., , 2018;;Veitinger et al., 2016).PRA models can be applied to define the spatial extent of release areas in dynamic avalanche simulations, which estimate the runout distance, velocity, and flow height of avalanche debris (Christen et al., 2010), or as a standalone spatial layer to assist with hazard identification.Their ability to operate at the mountain range scale with limited human input dramatically reduces cost and time to develop spatial data sets which can assist infrastructure planners and recreationists in making more informed decisions about their avalanche hazard exposure (Bühler et al., 2018).
The development of large-scale avalanche hazard indication maps in Switzerland has led to them being applied as a tool to help backcountry recreationists visualize terrain hazards and incorporate them into their trip planning process (Harvey et al., 2018).
The current state of the art methods for PRA modelling have been developed and validated in regions with widely available high-resolution DEM and forest cover data as well as long term records of avalanche observations (Bühler et al., 2018;Veitinger et al., 2016).However, the majority of mountainous regions in the world do not have freely available high-resolution DEM or forest cover data yet, and long term spatially accurate records of avalanche release are very rare.This seriously limits the application and local validation of PRA models around the world.
An additional limitation of existing high-resolution PRA models is that they do not account for the interaction between forest characteristics and avalanche release.For example, both the Bühler et al. (2018) and Veitinger et al. (2016) PRA models allow for forested areas to be excluded from PRA calculations based on a forest mask layer, but they do not explicitly capture forest avalanche interaction.This reduces the applicability of these PRA models in mountain ranges where a significant portion of the avalanche terrain is forest covered, such as in western North America.
To address these challenges and make PRA models applicable more broadly, the objective of this research is to develop a costeffective workflow for generating the required input datasets for the Bühler et al. (2018) PRA model using satellite data and open-source remote sensing methods.In addition, we present a relatively simple method for adapting the current PRA model to work in forested terrain.In the absence of long-term avalanche observations, we develop a novel approach for utilizing the expertise and terrain knowledge of local mountain guides to validate the PRA model output and optimize the input parameter for the unique terrain and snowpack characteristics of our study area.These three developments-the use of satellite data, the adaptation of the model to work in forested terrain, and the validation with local terrain expertise-together open new opportunities for applying state of the art avalanche terrain modelling in regions with limited existing datasets and resources.

Potential Avalanche Release Area Modelling
Early versions of GIS based avalanche terrain models (e.g., Gruber and Haefner, 1995;Maggioni and Gruber, 2003) struggled to outperform simple slope based avalanche release area estimates (Voellmy, 1955) due to the inability of low resolution DEMs (20-30 m) to detect small scale terrain features.Current PRA modelling methods evolved over the course of a decade and benefit from developments in high-resolution DEM production and remote sensing (e.g., Bühler et al., 2013Bühler et al., , 2018Bühler et al., , 2022;;Veitinger et al., 2016;Kumar et al., 2019).To define avalanche release areas the algorithms use different combinations of DEM derivatives (i.e., slope angle, terrain ruggedness, curvature, and aspect), which are calculated using focal functions of raster pixels.Bühler et al (2013) found that 5 m resolution is the optimal tradeoff between processing efficiency and smallscale feature identification for PRA modelling.With DEM resolution of 5 m, a common nine cell focal neighborhood (3x3) is 225 m 2 , which is well below the median slab size for human triggered avalanches of 4,000 m 2 (Schweizer and Lütschg, 2001).
Hence, high resolution input data is essential to capture sub-release area scale terrain characteristics which are critical for accurate potential release area modelling.
The development of these algorithms depends on a robust validation dataset of observed avalanche events to determine the optimal input parameter settings for the target study area.By comparing the extent of the PRA model output to the location of avalanche observations the overall accuracy of the PRA model can be evaluated and comparisons can be made between different combinations of input parameters.Such datasets can be created through recording of manual observations or generated by applying satellite mapping (Lato et al., 2012;Bühler et al., 2019;Hafner et al., 2021).The most comprehensive known avalanche release area validation dataset currently available is curated by the WSL Institute for Snow and Avalanche Research SLF in Davos, Switzerland with experienced staff manually mapping avalanche outlines throughout the winter in the surrounding mountain areas.This avalanche observation catalog began in 1970, and as of 2016 it included 5785 mapped avalanches (Bühler et al., 2018).This dataset is now expanded including data from satellite avalanche mapping (Bühler et al., 2019) as well as airplane (Bühler et al., 2009;Korzeniowska et al., 2017) and drone surveys (Bühler et al., 2017).
Using a subset of this validation data, Bühler et al. (2018) compared their PRA algorithm performance against another PRA model (Veitinger et al., 2016) and a simple slope based release area estimation method from Voellmy (1955).The validation study showed that the Bühler et al. (2018) and Veitinger et al. (2016) PRA models had lower probability of detection compared to the slope based model, by 3.5% and 2% respectively, but also had lower probability of false detection, by 13.3% and 8.4% respectively.This means that the slope only model detected a greater percentage of observed release areas compared to the PRA models but also overpredicted release areas at a higher rate than the PRA models.Due to the relatively smaller decrease in probability of detection (3.5% and 2%) compared to the decrease in probability of false detection (13.3% and 8.4%) the more sophisticated PRA models are considered more skillful than the slope based model.The relatively small improvements over the slope based release area estimate illustrates the fact PRA modelling is a field of marginal gains, but when applied over large areas marginal improvements can have a large impact on the output extent of PRA models.The results also show slightly better performance for the Bühler et al. (2018) PRA model over the Veitinger et al. (2016) PRA model, with lower probability of false detection and higher measures of Pierce skill score and Heidke skill score.
An additional advantage of the Bühler et al. (2018) PRA model is the ability to convert the raster based PRA model output to polygon features using object based image analysis.Converting the PRA model output to polygon features enables the PRA model to be paired with dynamic avalanche simulation software (Christen et al., 2010;Bühler et al., 2018)  The Bühler et al. (2018) PRA model has been applied in multiple regions worldwide, including Chile, Alaska, Afghanistan, and India.However, the input parameters have not been independently tested and optimized using local validation data.Therefore, it is unknown whether the input parameters optimized for Davos, Switzerland are appropriate for mountain regions with different topographic and snowpack characteristics.Our research aims to address this knowledge gap by applying an updated version of the Bühler et al. (2018) PRA model to the Columbia Mountains of southeast British Columbia, Canada and seeks to optimize the input parameters for the study area based on locally available validation data.

Avalanches in Forested Terrain
In addition to DEM derived terrain variables some PRA algorithms use forest coverage to define PRA based on the assumption that avalanche release is less common in areas with tall and dense vegetation.The snowpack in forested areas is generally more stable due to the anchoring effect of trees, forest canopy snow interception, the disruption of the continuity of weak layers due to snow drop from canopy, and altered snow surface radiation and temperature conditions.However, it is still possible for avalanches to release in forested areas, especially in areas with steep slope angles, low tree density, or in openings within forested areas (Bebi et al., 2009).Small and medium avalanches generally do not have enough impact force to damage trees or tree stands, and forests tend to reduce their runout potential by detraining snow from the flowing avalanche (Feistl et al., 2014).Larger avalanches can break or uproot trees and cause massive destruction to the forest ecosystem (Feistl et al., 2015;Bebi et al., 2009).The location of avalanche release areas in relation to the forest plays a large role in whether trees will impede avalanche flow or be destroyed and possibly entrained (Teich et al., 2012).
The ability to account for forest characteristics in avalanche terrain modelling is largely based on locally available data sets.
Laser scanning or LiDAR data provide high-resolution digital surface model (DSM) and digital terrain model (DTM) datasets to define the forest character, including canopy height, location and size of forest gaps, and basal area (Brožová et al., 2020;Dash et al., 2016).Vegetation height models derived from DSM and DTM data can be used to identify forests with protective function and input as forest masks in PRA models (Bebi et al., 2021;Bühler et al., 2018Bühler et al., , 2022;;Waser et al., 2015).Similar to their application for DEM production, the high accuracy of these data sets comes at a high cost.
Alternative lower cost methods for estimating forest characteristics include traditional field based sample plots and radar or optical remote sensing instruments (Hyyppä et al., 2000;Waser et al., 2015).The most accessible of these alternative methods is satellite based optical imagery, which can be used to create a forest land cover classification, to determine the extent of the forested area (Bühler et al., 2013), and can be combined with field plot observations of specific forest characteristics to create a predictive model based on the spectral and textural characteristics of the imagery (Dash et al., 2016).
Prior research has attempted to incorporate forest characteristics with PRA modelling (Sharp et al., 2018), but low resolution DEM and forest data combined with a limited validation data set make it challenging to evaluate the overall performance of the model.However, the principle of adjusting the potential for avalanche release based on forest character aligns with analytical and theoretical understanding of avalanche release in forested terrain (Bebi et al., 2009;McClung, 2001).This research aims to expand existing methods for capturing forest avalanche interaction in PRA models using satellite remote sensing methods that are cost-effective and efficient for processing large scale avalanche terrain models.

Methods
Applying the potential avalanche release area (PRA) model to the study area required three main analysis steps (Figure 1).First, developing a pipeline for producing high-resolution DEM and forest classification data from satellite imagery.Second, adapting the existing PRA model to better capture forested terrain and processing many versions of the PRA model using a predefined range of input parameters.Third, developing new methods to validate the PRA model using polygons collected from local experts in order to optimize the input parameters for our study area.Steps two and three required many iterations (Figure 1, Step 3c) to test different baseline input parameters and evaluate performance using our grid search validation procedure.The datasets and code required for replication of our DEM processing, forest classification, and PRA validation are available in our Open Science Framework (OSF) repository (Sykes et al., 2021).

Study area
The study area for this research is the tenure area of CMH Galena, a mechanized skiing operation that operates in the Selkirk Mountains of British Columbia, Canada, approximately 100 km southeast of Revelstoke (Figure 2).The tenure covers 1162 km 2 , ranges from 450-3,050 m in elevation and is composed of roughly 60% forested terrain.The Selkirk Mountains have a transitional snow climate with a maritime influence where persistent avalanche problem types are common.The most common persistent weak layers associated with these avalanche problems are surface hoar and faceted crystals associated with a crust (Hägeli and McClung, 2003;Haegeli and McClung, 2007;Shandro and Haegeli, 2018).The best existing DEM and land cover datasets for the study area are the Canadian Digital Elevation Model (CDEM) with a resolution of 18m and the 2015 National Land Cover Dataset (NLCD) with a resolution of 30 m.The resolution of both these datasets is too coarse for high-resolution PRA modelling.

DEM generation
Based on our desire to develop a cost-effective and reproducible approach for applying PRA models across large areas, we chose to purchase raw satellite imagery and use open source photogrammetry software to produce our own DEM.At the time we purchased the imagery our estimate was that producing our own DEM would be roughly 2-10x less expensive than alternative methods to acquire a 5 m DEM based on price quotes from multiple commercial suppliers.However, the cost savings of producing a DEM using raw imagery come at a tradeoff of requiring significant technical knowhow to process the stereo imagery.One downside of this approach is that the vegetation cover inhibits the ability to create a bare ground DEM (known as a digital terrain model; DTM) and we end up with a digital surface model (DSM) that represents the reflective surface at the top of the vegetation.While a DSM is not the ideal representation of terrain in forested areas (Brožová et al., 2020), the high cost of LiDAR, the only remote sensing method that can produce a DTM in vegetation covered terrain, currently prevents its widespread use.
Producing a 5 m DEM requires satellite imagery with a spatial resolution of at least 1.5 m.After comparing the products from various providers (Pleiades 1, Worldview 1-4, GeoEye 1, SPOT 6/7, and KOMPSAT 2-3) we purchased SPOT 6/7 imagery based on our requirements of DEM resolution, study area size, and cost.The listed price for tasking new imagery collection for 1.5 m resolution SPOT 6/7 tri-stereo imagery at the time of acquisition was USD $12.65 per km 2 for a minimum study area of 500 km 2 , which does not account for any academic or other discounts available through imagery suppliers.The SPOT 6 tristereo satellite images were captured on August 19th, 2019 with 1.5% cloud cover and no visible atmospheric distortions (wildfire smoke, haze) in the images.Tri-stereo imagery captures forward, nadir, and backward looking images in a single pass and provides three stereo image perspectives which increases DEM accuracy in steep terrain and minimizes sensor shading.For a more detailed description of our DEM processing interested readers should reference the supplementary material 'DSM production in mountainous, forested terrain using SPOT 6 tri-stereo imagery with Ames Stereo Pipeline'.
To improve and assess the accuracy of our DEM we collected a set of 66 ground control points (GCP) distributed across our study area using a Trimble Geo7x handheld differential global navigation satellite system (DGNSS) unit connected to an Hstar base station network, from August 24-27th 2019.We collected GCP in locations with high contrast such as edges of snowfields, water body inlets, bridges, and land cover transitions (e.g., boundary of talus slope and vegetation) to make the locations accurately identifiable in the satellite imagery.The timing of our image collection (August 19th, 2019) and GCP data collection (August 24-27th, 2019) meant that there were minimal changes in the natural features we used as reference points (i.e., snowfields, water bodies).
To process the imagery, we used a combination of open source software tools from Geospatial Data Abstraction Software Library (GDAL), QGIS, and the Ames Stereo Pipeline (ASP) version 2.6.2 (Beyer et al., 2018;GDAL, 2021;QGIS, 2021).
Several steps of preprocessing were necessary to optimize our images prior to stereophotogrammetry, including bundle adjustment and orthorectification (Shean et al., 2016).The ASP stereo tool was developed for imagery containing bare rock and glacial landscapes.Differences in image texture in forested terrain are challenging for the default settings of ASP to produce accurate pixel matches.To address this issue, we extensively tested different stereo correlation algorithms and stereo processing settings to optimize performance for forested mountainous terrain.Our best results were achieved using the smooth semi-global matching (MGM) stereo correlation algorithm (Facciolo et al., 2015), which resulted in fewer DEM holes in forested terrain and terrain with suboptimal lighting conditions.Optimizing the settings of the ASP stereo tool produced accurate pixel matches in forested terrain and was only limited by artifacts in the original imagery (cloud, cloud shadow, poor lighting conditions).
Our stereo processing workflow generated 6 separate DSMs from the SPOT 6 tri-stereo imagery by taking all possible combinations of left and right stereo images.The goal of this method was to reduce DSM holes in steep or poorly lit terrain by taking advantage of the multiple view angles provided by the tri-stereo imagery.Before combining the individual DSMs to produce the final DSM mosaic, we removed pixels with a triangulation error greater than the resolution of the input images (1.5 m) to ensure robust elevation estimates (Figure 3a).Overall, we see low normalized median absolute deviation (NMAD) values across the DSM mosaic (Figure 3b), with a median NMAD of 0.32 m.To improve the alignment of the final DSM mosaic to our GCP, we used the ASP point cloud alignment tool to co-register the output DSM to the GCP (Shean et al., 2016).To evaluate the accuracy of our DSM we used 15 internal checkpoints which were not used as part of our GCP dataset (Höhle and Höhle, 2009).
Localized cloud cover and poor lighting on steep north facing terrain caused several holes in our SPOT6 DSM mosaic, which account for approximately 1% of the total DEM area (11.7 km 2 ).We filled these holes by down sampling the existing Canadian DEM to 5m, aligning the CDEM to our SPOT6 DSM mosaic using the point cloud alignment tool in ASP, and then blending the two DEM datasets together.To avoid smoothing the entire SPOT6 DSM we progressively blended the datasets across a 60 m buffer from holes in the SPOT 6 DSM.
The methods described here were only tested on a single set of SPOT 6 tri-stereo images, but the performance in forested terrain was vastly improved compared to the default ASP settings.For more detailed information on the ASP workflow or the computer resources used to calculate the DSM please see the supplementary material or contact the authors.

Forest classification
The existing PRA model of Bühler et al. (2018) uses a binary forest mask based on photogrammetric vegetation height model classification to mask release areas in forested terrain.We tested several approaches to generate a binary forest mask for our study area.Since our SPOT 6 imagery was limited by poor lighting conditions on steep north facing terrain due to early morning sun angle, we substituted Planet Labs' RapidEye imagery, collected on July 14th, 2018 (Planet Team, 2017).An advantage of the RapidEye imagery is that it includes a red edge band which provides additional spectral resolution to differentiate between forests and other types of vegetation (Dash et al., 2016).
The overall accuracy of the classifier is critical for providing a distinction between forested land cover and other types of vegetation, such as shrubs and herbaceous plants.For avalanche release area modelling this distinction is important because shrubs and herbaceous plants are buried or pressed down beneath the winter snowpack and therefore have minimal effect on the potential for avalanche release.Trees with rigid trunks that resist being laid over by the winter snowpack and canopy heights greater than the snowpack depth (approximately 2-3 m) have an anchoring effect on the snowpack which is essential to capture accurately in order to account for their effect on avalanche release.By iteratively fine tuning the training dataset we were able to control how the classier identified forested terrain and opted to select a model that primarily captured densely forested areas and omitted areas with isolated smaller trees surrounded by shrubs and herbaceous plants.
To perform the classification, we used a random forest algorithm on the blue, green, red, red edge, and near infrared image bands utilizing the python libraries Numpy, GDAL, Rasterio, and SciKit Learn (GDAL, 2021;Gillies et al., 2013;Harris et al., 2020;Pedregosa et al., 2011).To improve the classification accuracy, we included the normalized difference red edge index (NDRI), normalized difference vegetation index (NDVI), and normalized difference water index (NDWI) as additional bands for the random forest classifier.Training data were created by manually drawing polygons around individual land cover types (forest, water, bare ground, snow and ice, shrub, moss and lichen) based on RapidEye, SPOT6, and Google Earth imagery from our study area.Our training dataset is composed of 253 individual polygons (12.0 km 2 ), with 73 polygons of forested terrain (3.6 km 2 ).For further details on the analysis methods used for the forest classification interested readers are referred to our OSF directory where the data and code are available for review (Sykes et al., 2021).

Integration of forest information into PRA model
Our development of additional PRA model functions to improve performance in forested terrain was guided by two principles; 1. Minimize additional complexity when running the PRA model compared to the original version.2. Utilize remote sensing datasets that are available in data sparse locations and do not require extensive field validation.

Forest density
To estimate forest density, we used a focal function to calculate the total number of forested pixels within a five-cell neighborhood (625 m 2 ).The function simply summed up the total number of forested pixels and did not account for the location of the forested pixels within the five-cell neighborhood.This step resulted in a forest sum raster with values ranging from 0 to 25, with 0 meaning no forested cells and 25 meaning all cells within the five-cell neighborhood are forested (Figure 4, step 3).We included this step to capture the fuzzy transition between forested and non-forested snowpack characteristics.In areas adjacent to forested terrain the snowpack can be altered by forest cover (i.e., wind dynamics, radiation balance, canopy snowfall interception) despite not being directly covered by the forest canopy (Bebi et al., 2009).This method also helps to identify glades or meadows within the forest canopy by creating a fuzzy buffer around small non-forested islands within densely forested terrain.The size of the neighborhood function (625 m 2 ) is representative of small human triggered avalanches that have the potential to bury or injure a person, especially if they are carried into a forested area (Schweizer and Lütschg, 2001).We then reclassified each forest sum cell into an ordinal variable with the forest density categories open (0 cells), sparse (1-10 cells), moderate (11-20 cells), dense (21-24 cells), and very dense (25 cells) (Figure 4, step 4).We chose this uneven classification scheme to bias the application of the forest slope scalar parameter towards increasing the slope angle minimum more strongly in densely or very densely forested areas (i.e., cells with 21 to 25 neighboring forested cells).Since areas with more surrounding forested pixels likely represent more mature forests, this approach captures the fact that more mature forests have a greater potential impact on avalanche release.The resulting forest density layer provides a foundation to control how forested cells are included in the PRA model.

Forest slope scalar
As an additional control on how the PRA model is applied in forested terrain, we introduced a forest slope scalar parameter to increase the slope angle minimum based on the forest density value.Applying this parameter assumes that steeper slopes are necessary for avalanche release in forested terrain, which is supported by prior research (Campbell and Gould, 2013;Schneebeli and Bebi, 2004).The rate of slope angle increase is controlled by the forest slope scalar parameter (0.0-2.0), which is applied as an exponent to the forest density value (0-4) and added to the slope angle minimum value (e.g., 30°).For example, a slope angle minimum of 30° and a forest slope scalar value of 1 would result in the following slope angle minimums for forested terrain: open (0) 30°, sparse (1) 31°, moderate (2) 32°, dense (3) 33°, very dense (4) 34°.Whereas a slope angle minimum of 30° and a forest slope scalar value of 2 would result in the following slope angle minimums: open (0) 30°, sparse (1) 31°, moderate (2) 34°, dense (3) 39°, very dense (4) 46° (Figure 5).Altering the slope angle minimum input parameter changes the starting position of the forest slope scalar function but does not impact the rate of increase for each forest density value.

Parameter tuning and validation
To develop a meaningful validation dataset in the absence of long term records of avalanche events, we collaborated with two CMH Galena guides, who each have decades of experience in the study area, to develop a novel method that takes advantage of their local expertise to optimize the PRA model for our study area.For technical details on our statistical calculations and processing workflow, our validation processing script and data necessary to reproduce our results are available in our OSF repository (Sykes et al., 2021).

Validation data collection
CMH Galena primarily operates on approximately 300 defined ski runs within their tenure.The runs range in size from 0.2-19.0km 2 and their locations have been mapped with polygons that outline the typical skiing terrain (Figure 2).The frequency of how often these runs are used varies significantly depending on terrain characteristics, weather conditions for flying, and snowpack conditions.To validate the PRA model, the two collaborating guides selected five runs (highlighted in Figure 2) based on their familiarity with the terrain, their representativeness of the terrain characteristics relative to the entire study area, and the balance of forested and alpine avalanche terrain contained in the runs.
The process of collecting validation polygons from the CMH guides was carried out on a custom designed website.The website platform enabled us to develop and present meaningful reference layers (e.g., satellite imagery, topo maps, terrain data, GPS tracks, heat maps) and provide the guides with multiple perspectives of the study area to assist with drawing the validation polygons.Both guides drew release area polygons for the five validation runs individually before creating a final consensus set of polygons in collaboration.Through the process of developing the validation data collection workflow we found that mapping the precise location of start zones based on personal recollection without being in the terrain at the time is extremely difficult.ThereforeTherefore, we developed an alternative method that would explicitly accommodate this uncertainty.Instead of forcing the participating guides to explicitly outline all avalanche release areas, we asked them to draw validation polygons around terrain features with similar characteristics (i.e.slope angle, forest density, ruggedness) and specify for each polygon what proportion represent potential release areas (0%, 25%, 50%, 75%, 100%) (Figure 6).Polygons of obvious probable release areas or non-release areas where guides had high confidence about their spatial extent were labeled with 100% and 0% respectively.Areas with scattered probable release areas, such as open forests with glades, where the identification of each probable release area would be cumbersome and unreliable, were marked as larger polygons and labelled with the estimated spatial proportion of the probable release areas (25%, 50% or 75%).Outliers, such as infrequent release areas with very low slope angles that require specific snowpack structures and weak layer types (e.g., surface hoar), were not included in the validation dataset in order to avoid biasing the validation dataset toward rare events that are not representative of typical conditions in the study area. .Our fuzzy approach to mapping probable release areas has several advantages.Foremost, accommodating uncertainty in the spatial extent of release areas is a requirement when relying on human memory to generate the validation data as specifying probable release areas with higher precision from memory is simply unrealistic.This method also accounts for the variability in release area extent that results from the dynamic nature of snowpack and weather conditions.The workflow also minimizes the effects of local errors in the reference layers that we provided the guides with to record their validation polygons.Specific examples of uncertainty caused by reference layers are variations in satellite imagery lighting due to sun angle and artifacts of the DSM generation process, such as over steepened slope angle values caused by transitions from forested to non-forested terrain.Our final validation dataset consists of 167 polygons across five runs with a total area of 8.42 km 2 , with sample sizes of 100% = 91, 75% = 23, 50% = 23, 25% = 18, 0% = 7, run polygons = 5.In locations where the polygons overlapped, we retained the highest proportion value of the overlapping polygons.The overlapping region was also clipped from the total area of the lower probability polygon.Locations within the run polygons that were not explicitly mapped by the guides were assumed not to be release areas.However, our validation approach differentiates between these implied and the explicit 0% validation polygons because we have more confidence in the latter.

PRA model grid search
In contrast to the raster based validation approach of Bühler et al. (2018), our validation dataset requires analysis on the scale of individual polygons.Since we do not know the explicit locations of the release areas in polygons with release area proportions of 25%, 50%, or 75%, we cannot directly compare the PRA model output to the validation polygons on a pixelby-pixel basis.Instead, we have to compare the total area within each polygon that is considered a PRA by the model to the proportion provided by the local guides.To calculate the PRA error we subtracted the proportion of the polygon that the PRA model determines as a release area (0-100%) from the release area proportion determined by the guides (0%, 25%, 50%, 75%, 100%) for each validation polygon.For example, if the PRA model output predicted that a polygon contained 60% PRA and the guides designated that polygon as containing 50% release area, then the PRA error would be -10%.This PRA error value is the basis of our grid search process and can range from -100% to 100% depending on whether the PRA model overpredicted or underpredicted the guides estimated release area proportion.
To properly reflect the validation data collection process in our analysis we also need to consider the hierarchical structure of assessment polygons collected from the local guides.The highest value validation data are the 100% and 0% polygons because they provide explicit spatial extents for PRA locations.These polygons are from locations the participating guides are most familiar with and have the highest level of confidence in.We therefore placed more emphasis on PRA model performance in these areas when selecting the optimal inputs.The validation polygons with the greatest uncertainty are the run polygons.They were not explicitly drawn by the guides and the absence of PRA within these polygons was implicit and not explicitly specified.
Hence, the accuracy of these polygons was weighted least in selecting the optimal PRA input parameters.To select optimal input parameters for the PRA model we performed a grid search as described by Bühler et al. (2018) using the following values: slope angle minimum (default 30°, range 20°-40°), slope angle maximum (default 60°, range 45°-65°), ruggedness window (default 9, range 3-15), ruggedness maximum (default 6.0, range 0.5-10.0),curvature maximum (default 6.0, range 0.5-10.0),forest density (default NA, range 0-4), forest slope scalar (default NA, range 0.0-2.0)(Table 1).It is computationally not feasible to test all possible combinations of input parameters, therefore we used a set of default parameters from Bühler et al. (2018) as a baseline and iterated over each parameter to analyze the impact on the accuracy of the model.
Based on validation using the guide polygons we systematically updated the default parameters to optimize the PRA model accuracy for our study area (Figure 1, Step 3c).The input parameters slope angle minimum, slope angle maximum, ruggedness window, ruggedness maximum, and curvature maximum are derived from the DEM (Figure 7 a-c).The forest density input parameter is derived from the forest mask (Figure 7 d).Selecting the optimal set of input parameters did not rely on any single statistic.Each PRA model iteration was compared using the mean absolute error (MAE), mean bias error (MBE), proportion of accurate polygons, and proportion of underestimated and overestimated errors.MAE values can range from 0 to 100, with lower values indicating a more accurate model.MBE values can range from -100 to 100, with 0 indicating a balance between positive and negative errors.Polygons were considered accurately predicted if the PRA error was within ± 12.5%, meaning that the area of the PRA model output and guide estimate were within a 25% range of each other which is equivalent to 1 step in the guides rating scale (0%, 25%, 50%, 75%, 100%).Underestimated and overestimated polygons are defined as having a PRA error greater than ± 12.5%, and polygons with a PRA error greater than ± 25% were considered severely overestimated or underestimated.
The accuracy statistics for each grid search iteration were calculated on the basis of the total number of polygons (n = 167).
We elected not to weight the statistics based on polygon size because the highest value validation polygons (0% and 100%) are generally the smallest.Selecting the optimal input parameters for our PRA model required evaluating performance across all these statistics and taking the structure of our validation dataset into account.
When selecting the optimal set of input parameters we erred on the side of a model that overestimates the extent of potential avalanche release areas, which is indicated by a negative MBE.We consider this an appropriate approach because the guides' polygons reflect only the avalanche conditions that they have experienced and recall.Despite their multiple decades of experience, the guides have not witnessed all potential combinations of snowpack conditions, which could cause avalanche release in uncommon areas.In contrast, the PRA model is a terrain based tool which aims to identify locations in the study area which have the potential for avalanche release independent of snowpack conditions.

Results and discussion
Since the context of the input data, parameter settings, and output from the original model are vital for evaluating the performance of our updated version of the PRA model, we combine the results and discussion into a single section.After presenting and commenting on the results, we conclude this section with an evaluation of some likely sources of error for our updated PRA model and share our thoughts on the limitations of a purely satellite remote sensing based method for capturing forest character in the PRA model.

Data preparation pipeline
The data preparation pipeline produced a 5 m resolution satellite DSM and forested land cover data set as input for the PRA model.Using 15 internal check points (ICP), the DSM accuracy can be described with a median vertical error of -0.43 m and normalized median absolute deviation (NMAD) of 4.72 m (Table 2).These accuracy metrics indicate good performance of the stereo DSM method, especially considering the rugged mountainous terrain across our study area and close proximity of steep slopes to some of the ICP.Compared to the best available existing DEM for our study area (18 m resolution CDEM), the SPOT 6 DSM provides vastly improved small scale terrain feature identification (Figure 8).This indicates that the red edge band was by far the most important imagery band to delineate forested pixels.
Creating the forested land cover classification using the same satellite imagery as the stereo DSM processing would be the most efficient workflow for producing the necessary input data sets for PRA modelling, because it uses the least possible input data and thereby minimizes data acquisitions costs and effort.However, in our study, we elected to utilize Rapid Eye imagery as an alternative due to better overall lighting conditions and improvements in accuracy, primarily due to the red-edge spectral band.The overall accuracy of our classifier and the feature importance of the red edge band highlight the strength of RapidEye imagery for forest classification modelling.
Our processing pipeline provides a cost-effective approach for creating high-resolution DEM and forested land cover classification data in remote and data sparse regions.Compared to alternative methods, such as LiDAR and commercial satellite stereo DEM products, purchasing raw satellite stereo imagery to produce a high-resolution DEM provides significant cost savings, control over the DEM generation settings, and produces a DEM product with sufficient accuracy (Kramm and Hoffmeister, 2019;Shean et al., 2016).The primary limitations are the inability to resolve bare ground terrain features, susceptibility to DEM holes due to cloud cover and lighting conditions, and degree of technical knowhow and computer processing resources required to convert the raw imagery to a DEM product.Despite these limitations, the processing pipeline enhances accessibility for high-resolution PRA modelling in remote regions.
An alternative approach, which has the advantage of decreasing the technical skills required to produce a stereo DEM while still having significant cost saving benefits over LiDAR, is to purchase an off the shelf stereo DEM from a commercial satellite imagery provider.Costs vary greatly depending on resolution, location, and whether archival imagery is available for a given study area.In our case existing DEMs or stereo imagery were not available in our study area, so the added costs of new image acquisition and processing made producing our own DEM more advantageous.Those interested in applying these methods to their own area should carefully evaluate costs of acquiring a 5 m DEM to assess the feasibility of high resolution PRA modelling.This trend strongly contrasts the other polygon types (Figure 109, panels b and c), which have higher percentages of accurate polygons for slope angle minimum values > 26°.For 0% and 100% polygons the percentage of accurate polygons declines steeply above 26° accompanied by an increase in severely underestimated polygons.The MAE and MBE statistics follow a similar trend, with relatively uniform values until 27° followed by steeply increasing error rates and positive bias for the remaining grid search inputs.The run polygons (Figure 10, panel c) have the highest accuracy with slope angle minimum greater than 31°.However, the percentage of severely overestimated polygons decreases drastically at 27°.Below 28°, the MAE and MBE have steeply increasing error rates and negative biases, respectively.Above 28° the curves flatten out and trend towards 0 for both MAE and MBE.
Our choice of a 27° slope angle minimum strikes a balance between PRA model performance for each polygon type with a priority towards optimizing performance on the 0% and 100% polygons, which are the most spatially explicit and have the highest degree of certainty.Setting the slope angle minimum lower than 27° would result in too strong of a bias towards minimizing underestimated errors which is not appropriate given our target of a frequent avalanche scenario.This is illustrated by a decrease in overestimated and severely overestimated polygons at a slope angle minimum value of 27° for the 25%, 50%, 75% polygon dataset as well as the run polygons (Figure 10, panels b and c).

Forest density and forest slope scalar
Determining the optimal value for forest density was the most straightforward of the three parameters we optimized because the percentage of accurate polygons, lowest MAE, and lowest proportion of underestimated polygons all occur at a density value of very dense (4) (Figure 9, panel f).Setting forest density to very dense ( 4) means that the PRA model is not restricted by any forest mask and the forest slope scalar is applied across the full range of forest density values.
Out of the three parameters we optimized, forest slope scalar has the least variation in percentage of accurate polygons, MAE, and MBE across the range of values tested in the grid search (Figure 9, panel g).This indicates that the PRA model performance is less sensitive to changes in forest slope scalar compared to slope angle minimum and forest density.However, setting this parameter to 1.25 helps to create a more balanced model by decreasing the number of overestimated polygons, which is illustrated in the upward trend of the MBE value.
Similar to slope angle minimum, we see a decrease in the percentage of severely overestimated polygons for the 25%, 50%, and 75% and run polygon datasets for higher values of forest slope scalar (Figure 11, panels b and c).This is a trade off with a slight decrease in the percentage of accurate polygons and increase of percentage of underestimated polygons for the 0% and 100% polygons (Figure 11, panel a).This is reflected in the 0% and 100% polygon MBE value of -0.81 at 1.25, which is relatively high compared to the other polygon types.Given our target of a frequent avalanche scenario this trade off is justified to create a balanced PRA model and account for the influence of forested terrain on avalanche release.

PRA model output and comparison
The final PRA model captures 57.5% (96 of 167) of the consensus validation polygon data set accurately, which we define as the PRA model predicted area is within ± 12.5% of the area specified by the guides for each validation polygon (Table 3).The remainder of the validation polygons were either underestimated 10.2% (17 of 167) or overestimated 32.3% (54 of 167), compared to the guides' consensus estimates of release area proportion.The MAE value is 18.2, which is a measure of the average error across all polygons.The MBE value is -10.9, which indicates that the PRA model errors are negatively biased towards overestimating release areas.This interpretation of the MBE value aligns with the skewed distribution of underestimated and overestimated polygons.Overall, we see improvements in the MAE, percent of accurate polygons, and percent of underestimated polygons using the locally optimized input parameters.The MAE for the present model is 18.2 compared to 21.4 for the 'Bühler 2018no forest mask' version, demonstrating a slight improvement in overall model error (Table 3).The present model improves the percent of accurate polygons by 11.7 percentage points over the 'Bühler 2018 -no forest mask' PRA model, which is a substantial improvement given the marginal gains observed in prior PRA model comparisons (Bühler et al. 2018).Similarly, the reduction of 14.8 percentage points for underestimated polygons between the present model and the 'Bühler 2018 -no forest mask' demonstrates the improved performance of the grid search optimization.These improvements can be attributed to optimizing the slope angle minimum and forest slope scalar input parameters using the local validation data.
The trade off of the optimized input parameters for the present model is a bias towards overestimation, which is indicated by the MBE of -10.9 compared to -3.7 for the 'Bühler 2018 -no forest mask'.This is also shown by the slight increase of 3.2 percentage points in overestimated polygons from the 'Bühler 2018 -no forest mask' to the present model.Producing a more negatively biased PRA model is in line with our mindset of creating a PRA model that errs on the side of overestimating observed release areas.In our opinion, the benefits of improved percentage of accurate polygons and underestimated polygons outweighs the downside of a slight increase in overestimated polygons.
The present model has a substantially lower slope angle minimum of 27° compared to the default value of 30° from Bühler et al. (2018), which results in a notable increase in the overall area of the PRA output due to expansion into lower angle terrain (Figure 12).The fact that the validation data led us to a substantial decrease in slope angle minimum is likely due to differences in the terrain and snowpack characteristics in our study area compared to the region of Davos in Switzerland where the model was initially validated.The avalanche character in our study area is prone to persistent avalanche problem types with the most common weak layers being either surface hoar or faceted crystals associated with a crust (Hägeli and McClung, 2003;Haegeli and McClung, 2007;Shandro and Haegeli, 2018).As a weak layer, surface hoar can release at lower slope angles and has increased potential to propagate across terrain features compared to other weak layer types (McClung and Schaerer, 2006).
Despite our aim of excluding outlier release areas with extremely low slope angles that are only capable of producing avalanches under very specific snowpack conditions from the validation dataset in order to target a more frequent avalanche scenario, the widespread influence of surface hoar as a weak layer in our study area still contributes to an overall lower minimum slope angle threshold.The fact that our validation data set and grid search approach produced a PRA model that also aligns with our theoretical understanding of the snowpack properties in our study area is an encouraging result.However, in terrain within the study area that is not prone to surface hoar development, such as alpine terrain with a high degree of wind and sun exposure, our PRA model is likely to overestimate PRA extent.

Potential sources of PRA model errors
Based on discussions with our collaborating guides and exploring spatial patterns of discrepancies between our validation data set and PRA model output, we have highlighted two likely sources of error in our PRA model.First is the limitation of using a relatively simple remote sensing based approach to account for forested release areas in the PRA model, which does not explicitly capture forest characteristics that are known to have a strong bearing on the interaction of avalanches and forest, such as crown cover, stem density, and gap size (Bebi et al., 2009;Teich et al., 2012).Second is the inherent uncertainty of relying on human experience to generate validation data, which can be subject to individual biases and faulty recollection.
Overall, we believe that the forest characterization is responsible for a larger portion of the PRA model error and is the most fruitful direction for future research to try and address.This section provides examples of these sources of error and discusses how we have attempted to minimize their impact on the PRA model accuracy.

Forest characteristics
To shed light on potential sources of PRA model errors we applied two different approaches that consider different spatial scales.First, we visualized the spatial patterns in the PRA errors for each validation run and consulted the local guides to provide their insight.Second, we extracted the terrain characteristics of the entire set of validation polygons and compared the distributions of the terrain characteristics based on the PRA error value.Both approaches yielded similar insight, which highlight the challenge of capturing forested avalanche release areas accurately using an approach based purely on satellite imagery.
Visualizing the patterns of PRA model error by validation run reveals concentrated clusters of higher PRA error on specific runs or subregions within runs (Figure 13).The 'Lunatic Fringe' run has by far the highest proportion of overestimated polygons out of the five validation runs, with 22 out of the 42 validation polygons being overestimated (Figure 13, panel a).
Based on information provided by the local guides, this run is characterized by a steep continuous face with several welldefined large avalanche paths dissecting mostly forested terrain.The forest is very dense and impassable for a guided group at the upper elevations of this run.In contrast, the 'Red Baron' run, which is located directly across the valley from 'Lunatic Fringe', contains lower slope angle terrain with a large proportion of mature forest (Figure 13, panel b).The forest has greater canopy height with widely spaced gaps between the individual trees.The forest canopy between each tree extends horizontally enough that the land cover classification is unable to detect many of the gaps on the forest floor.This run contains 7 out of 8 of the severely underestimated validation polygons, with the other polygon located in a forested area with similar characteristics on the 'Bandito' run.This further illustrates the guides' interpretation that the severely underestimated polygons on 'Red Baron' have high forest density and the limitation of our forest slope scalar approach for accounting for forested terrain with highly variable characteristics.It is important to note that the sample size of 'severely underestimated' polygons is small with only eight polygons.For context, the distribution 'severely overestimated' (n = 37) polygons also include high percentages of forest cover and forest density, which can be partially attributed to the dense and tightly spaced forested terrain on 'Lunatic Fringe'.
The PRA errors on 'Lunatic Fringe' and 'Red Baron' demonstrate the limitations of our approach in capturing the real world forest characteristics.Further improving the performance of the PRA model in forested terrain would require more detailed data sets such as LiDAR or a field based forest inventory which could capture additional forest characteristics such as stem spacing (Waser et al., 2015;Hyyppä et al., 2000;Dash et al., 2016), which are beyond the scope of this research.A notable publicly available source of LiDAR vegetation height measurements which could be used to interpolate forest height or overall biomass and potentially improve the performance of PRA models in forested terrain is the NASA ICESat-2, which collects LiDAR point measurements across the globe.The benefit of our method is to create cost-effective and high-resolution avalanche terrain maps based exclusively on remotely sensed data which can be applied in any location, regardless of remoteness or accessibility.For this purpose, our approach allows forested terrain to be captured in the PRA model on a basic level and broadens the range of avalanche terrain that the model can be applied to.

Uncertainty in validation data
One of the key differences in relying on local expertise for model validation is the necessity to incorporate uncertainty.There are two distinct types of uncertainty that are relevant for validating the PRA model: 1) Uncertainty in the accuracy of the observations, recollection, and experience of our collaborating guides, and 2) uncertainty in the reference datasets we provided them with to transfer their knowledge into spatial datasets and precision of polygon drawing.
In the case of guide observations, the primary sources of uncertainty in determining the location of avalanche release areas are the variability of avalanche conditions, how often the terrain is observed throughout the season, the guides recollection of avalanche events, and the potential for altered snowpack structure due to frequent guiding.These limitations are inherent to relying on human recollection as a source of validation data.However, our approach for capturing validation polygons from local experts accommodates these limitations by allowing for fuzzy boundaries in drawing polygons, collecting validation data from multiple guides independently, and intentionally minimizing the specificity that we ask the guides to label the release area proportions (0%, 25%, 50%, 75%, 100%).
The process for collecting validation data from our collaborating guides evolved through frequent back and forth discussions.
When applying the validation polygons to select optimal input parameters for our study area we accounted for the nature of the data collection by placing more emphasis on the performance of 0% and 100% polygons, which have the highest level of certainty for the guides and are the most spatially explicit.We also preferred input parameters that resulted in a PRA model that is biased toward overestimating release areas in order to account for the potential that the guides have not witnessed all possible combinations of snowpack and weather conditions in our study area, despite their extensive experience.
An example of how the guides' experience can influence our validation data set can be seen in the right half of the 'Rendezvous' ski run, where there are many severely overestimated validation polygons (Figure 13, panel c).According to our DEM, the slope angles in this area are predominantly in the low to mid thirties, which are within the range observed for human triggered avalanches (Schweizer and Lütschg, 2001).However, the guiding operation frequently uses this piece of terrain and intentionally manages the snowpack using skier traffic to minimize the potential for weak layers to form and persist on the surface (e.g., surface hoar).Frequent guiding use and intentional maintenance of weak layers can create a modified snowpack structure (Haegeli and Atkins, 2016) and has the potential to impact the guides' perception of release area potential.In areas where the guide's experience is largely based on modified snowpack structures there is a high potential for the PRA model to overestimate avalanche release compared to the validation data set.
While our workflow for collecting validation data from local guides was customized for our use case, these methods could be adapted to other professional communities such as avalanche forecasters or ski patrol.We recommended considering the following key principles for developing meaningful PRA validation datasets: 1. Identify a manageable size area to create the validation dataset that is representative of the terrain and snowpack conditions in the larger study area you want to apply the PRA model.
2. Solicit feedback from collaborators to identify sources of uncertainty in their ability to translate their local experience into polygons that can be compared to the PRA model output.
3. Incorporate that uncertainty into the validation process by allowing them to use fuzzy boundaries to identify potential release areas.
4. Take the structure of the validation data into account when performing statistical comparisons to the PRA model output.
This process can be time consuming and iterative, but it is critical to ensure shared understanding of the validation data between researchers and collaborators.In the absence of long-term observations of avalanche events in most mountainous regions the development of methods to extract local knowledge from human experts is critical to the application and validation of PRA models.

Limitations
The primary limitations of this research are direct consequences of our aim to minimize the cost of input data production and create a flexible workflow to apply and validate the PRA model in remote and data sparse regions.Using a DSM as input for a PRA model has not been thoroughly tested, and the inability to detect bare ground features within forest canopy likely causes localized errors in the PRA model.Recently, a comparison of high-resolution DSM and DTM models for avalanche runout modelling demonstrated some of the limitations of a DSM for dynamic avalanche simulation (Brožová et al., 2020).We were unable to test the accuracy of the SPOT6 DSM compared to a DTM due to the lack of alternative high-resolution data in our study area.
Relying exclusively on optical satellite imagery to account for forest avalanche interaction provides limited detail on meaningful forest characteristics.Explicit modelling of stem density, gap size, or crown cover could improve the PRA model's ability to capture forest avalanche interaction (Dash et al., 2016;Wallner et al., 2015).However, our focus is on minimizing field data collection to create a workflow that is applicable in remote areas.
Finally, the experience of local experts is not an ideal source to generate validation data compared to long term observation records.Observations from individual experts are prone to biases in their experience and potential for faulty recollection.We attempted to minimize these effects on our dataset by collaborating closely with the guides to develop a system for recording their observations that allows for uncertainty and is based on independent observations of multiple guides.

Conclusions
This research aimed to increase the range of application for existing high-resolution PRA modelling by developing a costeffective workflow for generating the required input datasets, expanding current PRA modelling methods to include avalanche Using locally optimized input parameters, our updated PRA model has a higher overall accuracy and less underestimated release areas compared to the default parameters developed for Davos, Switzerland in Bühler et al. (2018).Our validation approach utilizes local expertise to collect avalanche release area polygons via a custom-built online mapping tool and applies spatial and statistical analysis to quantify the accuracy of the PRA model.We leveraged this unique validation data set to develop a new polygon based grid search approach to optimize the PRA model input parameters.Creating a validation method that allows for optimization of the PRA model in areas without a long standing avalanche observation dataset is essential to evaluate the PRA model performance in new locations.This method also provides the opportunity for comparison of optimal input parameters in different snow and avalanche climates.Future research applying the PRA model in maritime and continental snow climates would provide additional insight into how the input parameters can be optimized for a broader range of snowpack and avalanche conditions, which are not captured in the existing Davos or Galena study areas.
To include forested terrain in the PRA model we focused on creating a simple addition to the existing PRA model which does not require any additional input data and remains an optional extension of the existing PRA model framework.We also focused on maintaining the ability to create the input data sets via optical satellite remote sensing methods to minimize the overhead cost and effort to produce forest characteristic data.Our approach allows the PRA model to capture the interactions between forests and avalanche release by controlling the forest density where the PRA model is applied and altering the slope angle minimum threshold based on the local forest density.These two changes are simple yet effective methods to account for forest cover in PRA modelling.
Additional research focused on satellite imagery based modelling of forest characteristics (Dash et al., 2016;Hyyppä et al., 2000), such as stem density and gap size, could further improve the performance of PRA models in forested terrain.While the availability of high-resolution LiDAR, laser scanning, or field measured forest characteristics are essential for meaningfully validating the derivation of these forest characteristics datasets (Waser et al., 2015), this type of development and analysis was beyond the scope of this research.The forest regions in our study area are dominated by coniferous tree species, which limits our ability to generalize the effectiveness of the PRA model in coniferous or mixed forest ecosystems.Hence, we encourage other researchers to explore our approach in other forest types.
Despite the limitations and shortcomings of our approach, the present research improves the accessibility of high-resolution PRA modelling by combining an existing state of the art PRA model with open source software tools and lower cost input data and presenting a flexible validation method to assess accuracy of the model output based on local terrain expertise.These developments have the potential to enable a more widespread application of high-resolution avalanche terrain indication modelling worldwide.
to estimate runout distance, impact pressures, flow depth and velocity of the avalanche flow.This powerful combination of release area and runout modelling represent the state of the art of current avalanche terrain indication modelling practices and are a valuable resource for large scale avalanche hazard indication mapping.Therefore, this research seeks to improve and expand upon the existing Bühler et al. (2018) PRA model.

Figure 1 .
Figure 1.Workflow diagram illustrating the necessary input datasets and processing steps to apply and validate the PRA model in our remote, data sparse, and forested study area.The dashed lines in step 3c indicate the option to either refine the baseline input parameters of the PRA model and re-start from step 2a or select the final PRA model and move to step 3d.

Figure 2 .
Figure 2. Study area map showing the CMH Galena tenure, lodge location, ski run polygons, and the subset of runs used to validate the PRA model.Forest data created using Planet Labs imagery (Planet Team, 2017), inset map made with Natural Earth.

Figure 3 .
Figure 3. SPOT 6 DSM error estimates.Triangulation error for each set of stereo pairs (a) with pixels where error is greater than image resolution (1.5 m) removed from DSM. Normalized median absolute deviation (NMAD) for mosaic of 6 stereo pairs (b) with inset map showing slope scale detail.Internal checkpoints (green points) with height difference in meters between DSM surface and DGNSS measurement (negative values indicate the DSM surface height is lower than ICP height).
The forest density parameter controls what classes of forest are included in the PRA model, while the forest slope scalar adjusts the slope angle minimum threshold based on the forest density class for each pixel.Including these parameters takes advantage of the existing forest mask functions of the PRA model and only adds two input parameters when running the PRA model, both of which are optional and can be omitted to run the PRA model in the prior configuration from Bühler et al. (2018).

Figure 5 .
Figure 5. Forest slope scalar functions applied to a 30° minimum slope angle threshold.

Figure 6 .
Figure 6.Validation polygons from one run at CMH Galena.Polygons are color coded based on the release area proportion of each polygon.Forest data created using Planet Labs RapidEye imagery (Planet Team, 2017).

Figure 7 .
Figure 7. PRA model input parameters.Slope angle, curvature, and ruggedness derived from the DEM (a-c) and forest density derived from the forest mask (d).Forest data created using Planet Labs RapidEye imagery (Planet Team, 2017).

Figure 8 .
Figure 8.Comparison of existing 18 m resolution CDEM to 5 m resolution SPOT6 satellite stereo DSM, derived from our data preparation pipeline.Canadian DEM data from Natural Resources Canada.

Figure 9 .
Figure 9. Results of PRA model grid search.In each of the panels, the left Y-axis shows the percentage of polygons in different PRA error classes with colored bars (accurateyellow, underestimatedred, overestimatedblue).Black squares and triangles show the values of MAE and MBE for each grid search iteration with a grey dashed horizontal line to show the 0 threshold which correspond to the right Y-axis.The vertical back lines indicate the optimized parameter settings.

Figure 10 .
Figure 10.Frequent avalanche scenario PRA model grid search results for slope angle minimum with validation polygons split based on the type of polygon.

Figure 11 .
Figure 11.Frequent avalanche scenario PRA model grid search results for forest slope scalar with validation polygons split based on the type of polygon.

Figure 12 .
Figure 12.Comparison of present PRA model to 'Bühler 2018 -forest mask'.Present model PRA area is dark purple with light purple for forested areas.'Bühler 2018 -forest mask' is shown in orange for comparison.Inset map shows detailed PRA comparison on a local scale, extent shown by black dotted line.Forest data created using Planet Labs RapidEye imagery (Planet Team, 2017).

Figure 13 .
Figure 13.PRA model accuracy for each validation run, with the downslope direction at the bottom of each panel.The validation polygons are labelled with their release are proportion and color coded based on the PRA error for each individual polygon.For overlapping validation polygons we retain the highest release area proportion value and clip that area from the surrounding lower proportion polygon.

Figure 14 .
Figure 14.Analysis of PRA error based on forest density and percentage of forested area for the validation polygon data set.The plots show the distribution of forest density and forest cover percentage for validation polygons based on their PRA error.
forest interaction, and by creating a novel approach for validating the model based on the local expertise of avalanche practitioners for data sparse regions.The research produced an updated version of theBühler et al. (2018) PRA model which enables high-resolution avalanche terrain modelling in a vastly greater proportion of mountainous terrain than previously possible.This is thanks to the widespread availability of the necessary satellite remote sensing input data and local expertise required to validate and optimize the PRA model input parameters.The updated model also allows for inclusion of forested terrain with varying densities, contributing to a substantial improvement in the performance of the PRA model in our study area.The data preparation pipeline developed for this research is based on open source software and intended to be reproducible in areas without existing high-resolution DEM and forest cover data sets, which achieves our goal of making high-resolution PRA modelling more accessible in remote and data sparse areas.Producing a satellite stereo DSM based on raw imagery provides control over the DSM characteristics and minimizes the cost associated with acquiring this essential data set.Further testing of the DSM pipeline developed for this research is required, especially in forested terrain, and could provide a meaningful direction for future research.Despite the dramatic cost reduction of our workflow, high resolution satellite stereo imagery areis still relatively costly so readers interested in applying PRA models in their own area should carefully evaluate costs of acquiring the necessary input data.

Table 3 . PRA model comparison
Bühler et al. (2018)our parameter optimization demonstrates meaningful improvement, we compared the accuracy statistics of the model using the optimized parameters (Present model) to theBühler et al. (2018)defaults both with and without a forest mask (Table3).The 'Bühler 2018 -forest mask' PRA model does not identify release areas in any terrain identified as forested based on the land cover classification, whereas the 'no forest mask' version allows the PRA model to calculate release areas in all terrain.Since the 'forest mask' version naturally performs substantially worse in most accuracy statistics due to the large proportion of forested terrain in our study area, we will focus the comparison on the 'Bühler 2018 -no forest mask' model version.