Introduction
Avalanche-warning programs in the United States and Europe are currently operated without the benefit of models. Forecasts are made by the judicious use of both experience and a basic knowledge of the physical processes affecting snow stability. This technique is flexible because it allows a synthesis of the important snow-structure and snow-texture variables together with pertinent weather data, avalanche response, and antecedent conditions. Heavy reliance is placed on the avalanche activity which obtains when the forecast is made, because this is one of the few definitive ways of determining stability: it also drastically reduces the number of false alarms. The main weaknesses of such programs are their heavy reliance on timely and accurate avalanche data, their high degree of subjectivity, the limited number of qualified personnel, and the immense task of remembering antecedent conditions for large, remote mountain areas.
Future use of statistical and simulation process models for guidance will eventually alleviate some of these problems. Current research efforts arc directed toward this goal, but progress has been slow due to the complexity of the problem. Reference PerlaPerla (1970), Reference Bois and ObledBois and Obled (1973), Reference Judson and EricksonJudson and Erickson (1973), Bois and others ([1975]), Salway (unpublished), Reference BovisFöhn and others (1977), and Reference BovisBovis (1977), made encouraging progress using techniques ranging from the use of scatter diagrams through discriminant function analysis to time-series modeling. All investigators had difficulty with important variables, including snow structure, snow texture, and wind. Snow-cover variables have always been difficult to analyze and are seldom available on a daily basis. Wind was not selected as a variable in most stepwise multiple regression techniques because it failed to discriminate between avalanche days and non-avalanche days. Frequently, snow available for transport is removed long before the wind drops its threshold speed. Use of a snow-transport model solves this problem. Because of difficulties associated with statistical treatment of snow and avalanche variables, we tried simulation modeling with a process-oriented model. The basic premise is that snow stability in avalanche starting zones is related to snow structure, texture, and loading rates, after minimum depth requirements are satisfied.
Modeling Technique
Our simulation process model uses the approach developed by Leaf and Brink (1973[a], 1973[b]) to forecast streamflow from Rocky Mountain alpine and subalpine zones. It uses the primary subroutines in the streamflow model. The streamflow model simulates: (1) winter snow accumulation, (2) the energy balance, (3) snow-pack conditions in terms of its energy level and free-water requirements, and (4) resultant snowmelt in time and space on all combinations of elevation, aspect, and slope. The snow-pack is assumed to behave as a dynamic heat reservoir, so input elements are expressed in units of heat. Precipitation events are added algebraically as caloric heat to develop the heat reservoir or snow-pack. Midpack temperatures are computed using unsteady heat-flow theory. Snowmelt is released only when the pack becomes isothermal and its free-water holding capacity is saturated.
The original streamflow model was modified and extended to simulate layer depth, age, densification, temperature gradient, temperature-gradient metamorphism, melt-crust formation, redistribution of snow by wind, loading rates, and avalanche danger. Six-hour indices of avalanche danger are based on snow structure, texture, and weighted average load values over the previous 72 h. Load values in the previous 24 h receive full weight; lesser weights are assigned to older loads. The model handles wind (loading) from any direction.
The model simulates snow-pack stability on one or more avalanche units selected by aspect and terrain to represent avalanche paths characteristic of a given area. The avalanche unit (Fig. 1) consists of an up-wind fetch, a down-wind deposition area (starting zone), and a down-wind zone representing the avalanche track. An idealized ground profile represents the major physical features of each unit. Specific numerical descriptors specify aspect, slope, and length of the fetch and deposition areas. The fetch is treated in three parts or cells to account for variations in the amount of snow trapped by vegetation and terrain irregularities along the fetch. The starting-zone area has a simplified triangular configuration with a maximum accumulation height H which varies according to the slope angle ϕ. Twenty-nine cells, each one meter long, are located in this area. Two cells maintained in the track give the unit a total of 34 cells. The 34-cell total is caused by limitations in the computer memory, but programming techniques allow an extension of the starting zone beyond 29 m with little impact on the simulation.
Fig. 1. The idealized avalanche unit consists of an up-wind fetch, a down-wind deposition area (starting zone), and a down-wind zone representing snow conditions in the avalanche track. Thirty-four cells, representing various segments of the unit, are depicted. Cell 34 represents snow conditions in a level study plot.
Input Data
Weather, snow, and avalanche data are from the Forest Service avalanche station at Berthoud Pass, Colorado. Temperature, precipitation, and snow depth were measured in a sheltered clearing at 3450 m. Windspeed and direction come from instruments on a nearby ridge at 3620 m. Avalanche data are from 12 paths situated within 1000 m (horizontal) of the weather sites. The 6 and 24 h inputs are:
Subroutines
General
Figure 2 depicts the general flow of the model and displays the primary subroutines discussed. They deal with the seasonal snow-cover characteristic of the unit. Special subroutines for a permanent snow cover were developed for use on other units, but are beyond the scope of this presentation.
It should be noted that the 6 h avalanche danger ratings derived from the load analysis are based on both 6 and 24 h subroutines. For this reason, the avalanche danger rating subroutine will be discussed last.
Layers
To simulate the development of a layered snow cover, 6 h precipitation and snow- transport events are deposited in layers of water-equivalent with a minimum thickness of 5 mm. Falling snow and windblown snow must be treated separately, because programming Fig.2 techniques do not allow simultaneous treatment of both events. Snow-fall is deposited at the end of each 6 h period. If the water-equivalent equals or exceeds 5 mm, a layer is completed and a new layer will begin with the next event. When less than 5 mm falls, the layer fraction or sub-layer is stored until additional water-equivalent completes the layer. In the case of snow transport, layers are decremented and incremented on the fetch, and are incremented (occasionally decremented) in the starting zone. Layer age is stored for later use in computing avalanche danger, and in the densification subroutine which simulates snow depth. A maximum of 100 layers can form in every avalanche unit cell.
Fig. 2. Generalflow of the simulation model with primary inputs and process steps.
4.3 Snow transport
(a) Wind direction
Wind direction is introduced when the speed exceeds the threshold value of 7.3 m/s in the model. Wind direction is specified to (1) have no effect on the unit, (2) transport snow out of the unit, or (3) transport snow within the unit. In the first case, wind is not permitted to move snow in the unit due to obstructions in flow such as an arête or forest. In the second case, snow is transported out of the unit by cross-winds. There are two options in the third case. First, snow may be removed from one part of the fetch (cell 1), deposited on another (cells 2 and 3), removed from cell 3, and deposited in the starting zone. Second, snow can be removed from the starting zone and deposited on fetch cell 3.
(b) Theory
Theory for this part of the model is based on the work of Reference Budd and RubinBudd (1966), Reference Budd and RubinBudd and others (1966), Reference SchmidtSchmidt (1972), Reference Tabler and SchmidtTabler and Schmidt (1973), and Reference TablerTabler (1975). Basic equations are presented but not derived.
The volume of the snow transported from the fetch (cells 1 through 3 in Fig. 1) is:
where Q is the volume in units of water-equivalent per unit width perpendicular to the wind, Pr is the water-equivalent on the fetch available for relocation, in mm, Rµ is the transport distance before an average sized particle evaporates, in m, and Rc is the fetch length, in m.
Water-equivalent available for relocation (Pr) from the fetch is given by the expression:
where P is the water-equivalent of snow on the fetch, and Ps is the sum of the water-equivalent of snow trapped by vegetation, plus that which is stored by terrain irregularities, plus that which melts and evaporates. The amount of snow stored by vegetation and terrain irregularities is specified by field estimates.
There is a complex relationship between environmental conditions and the transport distance Rµ which is the distance an average-sized snow particle travels before it completely sublimates. It depends on wind-speed and sublimation processes. According to Reference Tabler and SchmidtTabler and Schmidt (1973):
Equation (3) is expressed in terms of 10 m wind-speed U10 (in m/s), temperature T (in °C), relative humidity RH (in per cent), and time t (in hours) since the last snow-fall. Transported snow is available for deposition when 
and when more snow is available for transport than required by the up-wind storage given in Equation (2).
The rate of snow transport F, in g/m s, given by Reference Budd and RubinBudd and others (1966) is:
where
Equation (4) is programmed into the model with two modifications. First the to 10 m windspeed is modified by a simple trigonometric function so as to approximate changes in wind-speed along the fetch due to slope. Then, to restrict the flux when surface melt crusts are present (F limits Q):
where α is an erodibility index based on crust thickness. Surface crusts reduce transport rates sharply until such time as the crusts themselves are eroded by strong winds.
Snow deposition (wind-transported snow)
Deposition areas in the avalanche unit include cells 2 and 3 on the fetch and cells 4 through 32 in the starting zone. Fetch cell 2 is 100m long and fetch cell 3 is 10 m long. For modeling purposes all cells are one meter wide and starting zone cells are one square meter in area, and have a finite volume obtained by specifying the maximum drift height above ground. Additionally, cells 1 to 32 are given an initial holding capacity based on field estimates. No transport takes place from a cell that contains less snow than its holding capacity.
To approximate the deposition process in avalanche terrain, it is necessary to make gross assumptions based on a general knowledge of trapping efficiency. In general, trapping efficiency is a function of the angle subtended by the slopes of the erosion and deposition surfaces, the wind-speed over the area, and the amount of snow already in traps. These factors vary with each transport event. For our purposes, the following efficiency estimates were made:
On the relatively shallow-angled fetch surfaces, snow is allowed to accumulate in uniform layers, whereas loading in the steep starting zone is specified to occur in vertical, wedgeshaped layers. Snow transported from fetch cell 1 is available for deposition (after sublimation) on cell 2, and to a lesser extent on cell 3. Snow from cell 2 does not reach the starting zone cells until it is completely filled. When this occurs, the entire fetch is treated as cell 3.
Snow is deposited into starting zone cells using a 6 : 3 : 1 ratio pattern (after Berg and Caine, unpublished). The total amount available for deposition Q is deposited beginning with 60% in cell 4, 30% in cell 5, and 10% in cell 6. If cell 4 will not hold 60% of Q, it is filled, Q is reduced, and the 6:3:1 ratios are distributed into cells 5, 6, and 7 (see Fig. 2 for more detail). The sequence continues in this manner throughout the season. Discontinuities which occur in filled cells due to settlement and evapotranspiration are refilled by further deposition. Maximum cell heights can be exceeded briefly if snow falls on filled cells when wind-speeds are below the threshold value. The excess is sheared off, added to the total amount Q, and deposited with the next transport event.
Radiation balance
This subroutine is used without modification as it appears in the streamflow model described by Leaf and Brink (1973[a]). Incoming solar radiation is adjusted for slope and aspect according to tables published by Reference Frank and LeeFrank and Lee (1966). The long-wave components are derived using the average air temperature in the Stefan-Boltzmann function. Reflectivity of the snow-pack is changed according to the energy balance, precipitation, and time (Reference ArmyU.S. Army, 1956).
Snow temperatures
The total heat energy input to the snow cover is indexed by the air temperature measured about one meter above the snow surface. For all practical purposes, studies have shown that temperature in this boundary layer approximates the integrated effects of wind, temperature and humidity of an air mass, and the radiation balance. The heat-flow equation solved by this method is given by Reference QuickQuick (1967) and Reference Riley, Riley, Chadwick and EgglestonRiley and others (1969):
where ks is the thermal conductivity, cs the specific heat capacity, ρs the snow density, Ts the snow temperature, Z the depth below the snow surface, and t is the time. Equation (6) can be rewritten as:
where KƲ is the thermal diffusivity.
Thermal diffusivity is assumed to vary with density as (Reference SchwerdtfegerSchwerdtfeger, 1963):
where ki is the thermal conductivity of ice, ρi the density of ice, and ci the specific heat capacity of ice.
Density is assumed to be constant across the pack for snow-temperature simulation, although it may vary with time. For our purposes, a weighted-average layer density, which is computed from the model, is used.
A finite-difference solution following the procedure outlined by Reference SmithSmith (1965), and Reference Richmyer and MortonRichmyer and Morton (1967) is obtained for Equation (7) using average air temperature, a preset ground temperature, and previous snow temperature at each node to simulate present snow temperatures. The number of simulated temperature nodes is dynamically proportional to snow depth with a spacing on the order of 15 cm.
Temperature-gradient metamorphism
Development of temperature-gradient snow is simulated using the expression developed by Reference Giddings and LaChapelleGiddings and LaChapelle (1962):
where t is the formation time in days, D the final crystal diameter in mm, IIa the site specific atmospheric pressure in mbar, II0 the sea-Ievel pressure in mbar, and dTs/dz the snow temperature gradient in deg cm-1.
The rate of formation is obtained by specifying initial and final crystal diameters Dand solving for t. For modeling purposes we initialize D at 150 µm for wind-transported snow and 600 µm for newly fallen snow. Growth is incremented according to Equation (9) using snow temperature and snow temperature gradients provided by Equation (7). When crystal diameters reach an arbitrary size of 0.7 mm they are designated as temperature-gradient (TG) snow. Thickness of TG layers is monitored every 24 h. A high ratio of this layer thickness TG to the total snow depth HS increases avalanche danger as will be explained later. Routines for equi-temperature metamorphism were not developed.
Layer densification
The model uses five density functions to simulate snow depth and layer densification throughout the year. Density-age relationships were derived from long-term field data from Berthoud Pass, Colorado and Alta, Utah. The function chosen depends on time of year, snow structure, and load rate. The functions are:
where x is the number of days. Function 4 is also used to account for rapid densification when young, low-density layers are loaded by heavy deposition.
Avalanche danger
(a) Methodology
The development of avalanche danger in the model closely parallels techniques used by operational forecasters except that data on current avalanche activity and the reaction of snow to destructive tests are not used. While this information is important in operational programs,. we found no practical method of incorporating these inputs into the model due to problems of quantification and reliability on a short-term basis (i.e. every 6 h).
When current avalanche data are unavailable the basic approach is to (1) determine if snow depth and stratigraphy satisfy initial conditions for avalanching,(2) evaluate existing snow stability on the basis of structure and load, and (3) estimate the effects of either a change in the present load or in the existing snow strength. Information for such determinations is available from the present model (i.e. snow depth, basic stratigraphy, and load in starting zone cells 4-32). Initial conditions for avalanching in this unit are:
(b) Theory
Perla and Judson (unpublished) developed a theory for assigning weights to load increments with time. Given a continuous record of load as a function of time L(t), the stability index S should behave as:
where M(t) is a positive monotonically decreasing function of time, and M(o) = 1. Equation (10) states that the stability index S(t) evaluated at time t consists of a weighted sum of the load in the past interval o ≼τ≼ t.
Discrete data points were used to approximate Equation (10) as:
where 1 ≼ a1 > a2 a3 >.…
The coefficients a1, a2, a3,… were determined from existing precipitation and avalanche occurrence data using a discriminant-function analysis on data collected at several sites on the Forest Service’s West-wide network (Reference JudsonJudson, 1970). Results of this statistical work showed that loading events more than 72 h earlier had a negligible contribution to avalanching on any given day. The data agree with work done by Reference ObledObled (1970) and Reference Yefimov and KozikYefimov and Kozik (1974). Actual load weighting factors used in the model are: 1 for the first 24 h, 0.6 for 25-28 h, and 0.2 for 49-72 h. In subroutine load analysis, the program scans load cells 4-32 and selects the two cells containing the greatest weighted load over the previous 72 h. A mean value for the two cells is then scaled, and this value gives the load factor A for the avalanche unit every 6 h.
(c) Basal instability
The final step in developing avalanche danger for the model involves linking the load factor A with structural weaknesses in the pack caused by the presence of temperature-gradient crystals. Examination of a large number of snow profiles taken in avalanche starting zones provided a means for incorporating this phenomenon into the model. The profiles were collected by project personnel and by members of the San Juan Avalanche Project (Reference Armstrong and IvesArmstrong and Ives, 1976). The ratio of the depth f TG to the depth of the snow-pack HS was chosen to index this important snow-pack feature. The ratio is called the basal instability factor BIF. BIF values range from 1 where TG/HS = 0.1 to a maximum of 2.1 following an S-shaped curve fitted by the expression:
where TG/HS > 0.1. The model optionally limits BIF to 1 once snow depth exceeds 2 m.
(d) Avalanche-danger ratings
The avalanche-danger rating ADR is the product BIF A depicted in Figure 3 on a 24 h basis from December 1973 through February 1974. ADR is set at zero until initial conditions with respect to snow depth, age, and density arc satisfied. At this point a background ADR of 0.1 is assigned to account for minimal danger which persists through the winter. Six-hour ADRs are computed using the current six-hour load factor A and the previous day’s BIF because layer thickness TG comes from a 24 h subroutine.
Fig. 3. Comparison of computed versus observed avalanche danger, December-February for winter 1973-74.
Observed Avalanche Danger (Observed Avalanches)
Twelve east-facing avalanche paths are in this first unit. Five paths are uncontrolled, four receive light ski control, and the remaining three receive moderate to heavy control. Excellent weather and avalanche data arc available from 1962-79. Observed avalanche danger in Figure 3 is a function of the number of avalanches, sizes, types, and avalanche-path constants. The path constant is a product computed as the inverse of the path frequency and its control intensity. For example, the smallest path constant of 0.2 was assigned to the path with the highest frequency and the heaviest control. The remaining path constants range from 0.3 to 1.2. Maximum possible daily avalanche danger for the subunit is 63.6, and the greatest observed one-day value was 15.6 for the period. This formulation of the predictand allows use of natural and artificial avalanches, and accounts for the severity of avalanching by combining several aspects of the phenomena.
Results
Baseline model runs were made for the eight winters, 1970–77.model predicted avalanche potential on 86% of the 175 avalanche days during the initial test period. For our purposes, avalanche potential is defined as any computed avalanche-danger rating greater than the background ADR of 0.1. The model indicated avalanche potential 50% of the time on non-avalanche days, for an overall error rale of 45%. A sensitivity analysis is currently under way to improve model performance and test the validity of the process steps in this prototype model. For the time being, initial results appear promising, and avalanche danger for the unit will be simulated on a real-time basis during the winter 1979-80. Model output will be available every 6 h in Colorado’s avalanche-warning center for direct comparison with conventional techniques.
Future Work
Considerable work remains in basic research and analysis of present model capability. Additional basic research is needed on threshold windspeeds,trapping efficiency,depositional mechanisms,and development of equi-temperaturc algorithms. Analysis of present model capability has just begun,using a fraction of the data available for determining snow stability. Some intriguing options include assessment of the load factor A over more than two cells,introduction of subroutine AVALANCHE for ADR adjustments,use of additional snow-cover variables,and adjusting numerous parameters present in the prototype model.
Upon completion of the analysis for this avalanche unit,additional units will be added to represent avalanche danger for an entire region. Further analysis will entail a comparison of simulated responses of the physical processes affecting snow stability with avalanche data from all representative units. If this work produces satisfactory results the model will be extended to other avalanche areas in the west.
Acknowledgement
The authors express sincere thanks to R. A. Schmidt of the Fort Collins Mountain Snow and Avalanche Project for his interest and technical support in the development of the snow transport and deposition sections of this paper.
