Application of Numerical Transient Fluid Dynamics to Snow Avalanche Flow. Part I. Development of Computer Program Avalnch

Abstract

A two-dimensional, transient fluid-dynamics computer code has been modified for specific application to the avalanche-runout problem. This code, called AVALNCH, permits the separation of path geometry effects from such flow factors as friction and viscosity. The longitudinal profile of the avalanche path is divided into cells, 10 to 20 m long, each of which can be assigned specific values for slope gradient, surface friction, and internal kinematic viscosity. The program gives average avalanche speed cell-by-cell down the path and the location and depth of avalanche debris. Internal kinematic viscosity and surface friction were modeled on an avalanche path of simple geometry and were found to be of about equal significance in predicting runout distance. Additionally, surface friction is represented by an exponentially increasing function as speed decreases in the runout zone, in order to model observed avalanche terminal-motion characteristics. Program AVALNCH is reduced to a basic algorithm that is efficient to run, and contains the essential mechanics to model avalanche flow accurately. The most pressing need is for more physical data to permit the matching of program output to observed results under a variety of conditions.

Résumé

On a modifié un programme informatique bi-dimensionnel pour 1’écoulement dynamique à surface libre d’un fluide, en vue de l’adapter au problème du parcours des avalanches. Ce programme, appelé AVALNCH, permet de séparer les effets de la géométrie du couloir et ceux de facteurs d’écoulement comme le frottement et la viscosité. Le profil longitudinal du couloir est divisé en cellules de 10 à 20 m de long, à chacune d’elle étant assignées des valeurs particulières pour la pente, le frottement superficiel et la viscosité cinématique interne. Le programme donne la vitesse moyenne de l’avalanche, cellule par cellule le long du couloir, l’emplacement et l’épaisseur du culot de l’avalanche. Les coefficients internes de viscosité cinématique et de frottement superficiel ont été essayés sur un couloir de géométrie simple et on a trouvé qu’ils étaient d’à peu près la même importance pour la prévision des distances d’arrêt. De surcroît, le frottement doit être représenté par une fonction qui croît exponentiellement quand la vitesse décroît dans la zone de dépôt, en vue de reproduire les caractéristiques du mouvement terminal de l’ avalanche. Le programme AVALNCH se réduit à un algorithme de base qui est efficace pour le mouvement et contient les mécanismes essentiels pour simuler l’écoulement des avalanches avec précision. Le besoin le plus urgent est de disposer de plus de données physiques pour permettre un calage du programme sur les résultats réellement observés dans des conditions variées.

Zusammenfassung

Entwicklung des Rechenprogramms AVALNCH. Ein Rechenprogramm für die zweidimensionale Dynamik nichtstationärer Strömungen wurde für die spezielle Anwendung auf das Problem des Lawinenabgangs eingerichtet. Dieses AVALNCH genannte Programm gestattet die Trennung der Einflüsse Bahngeometrie von Fliessparametern wie Reibung und Viskosität. Das Längsprofil der Lawinenbahn wird in Abschnitte von 10 bis 20 m Länge unterteilt, deren jedem spezielle Werte für Hangneigung, Oberflächenreibung und innere kinematische Viskosität zugeordnet werden können. Das Programm liefert abschnittsweise die mittlere Lawinengeschwindigkeit längs der Bahn sowie Lage und Mächtigkeit des Lawinenmaterials. Die innere kinematische Viskosität und die Oberflächenreibung wurden einer Lawinenbahn mit einfacher Geometrie angepasst; sie erwiesen sich als etwa gleich bedeutungsvoll für die Berechnung der Reichweite. Zusätzlich wird die Oberflächenreibung durch eine Exponentialfunktion erfasst, die mit abnehmender Geschwindigkeit im Endbereich der Lawine zunimmt, um den beobachteten Merkmalen der Bewegung im Endbereich gerecht zu werden. Das AVALNCH-Programm wird auf einen Grundalgorithmus reduziert, der den Ablauf erfasst und die wesentlichen Mechanismen zur genauen Beschreibung des Lawinenabgangs enthält. Es besteht dringender Bedarf nach weiteren physikalischen Daten, um die Rechenergebnisse mit beobachteten Werten unter verschiedenen Bedingungen vergleichen zu können.

Introduction

The ability to predict the runout distance for a snow avalanche is an important factor in the establishment of zoning restrictions in mountainous regions. An existing method for runout prediction by Voellmy (1955) is based upon the equations for open-channel equilibrium flow. Using this method, an avalanche slope is modeled by a sequence of straight-line segments to which the equations are applied after the selection of values for equivalent surface friction and fluid turbulence. These parameters are adjusted so that they have their intended influence on the flow, and account for geometric and flow anomalies that are characteristic of the particular section of the slope under evaluation. Accurate application of this method is generally dependent upon the user selecting appropriate values for surface friction and fluid turbulence. Indeed, a variation of almost an order of magnitude in these parameters is possible, depending on the application (Schaerer, [1975]). The method is particularly difficult to apply when snow runs out onto a shallow slope or a flat terrain. Efforts to verify and refine the theory experimentally help to establish a more definite basis for its application (Salm, 1968).

Here we report the application of a numerical fluid-dynamics approach to the problem as we recognize that flow over continually changing terrain has a transient behavior. The basic flow parameters can be varied at each incremental distance along the path, and geometry variations are treated separately from the assignment of flow parameter values. The numerical algorithm is based upon an Eulerian finite differencing of the time-dependent two-dimensional Navier-Stokes equations, to provide pressure and velocity distributions. These distributions are iteratively updated to satisfy compatibility requirements, which in turn require flow advancement consistent with impressed forces and constraints. A fluid-dynamics code, designated SOLA, recently reported by Hirt and others (1975), has the essential features which are needed to model avalanche flow. These features include the representation of a free-surface geometry, the incorporation of non-zero fluid viscosity, and the modeling of a changing lower-boundary profile. Problems of the long-duration modeling of a snow avalanche, and the economics of the modeling scheme have been addressed in the development of a special-purpose code entitled AVALNCH. This code incorporates essential features of the SOLA code, combined with major modifications tailored to the specific requirements of the avalanche-runout problem.

In Part I of this report we include descriptions and discussions of the basic modifications made to the SOLA code in order to apply it to avalanche-runout prediction. The test cases which were used to establish values for the friction and viscosity parameters are outlined. Finally, we include a description of data preparation for program AVALNCH using information from topographic maps. In Part II of this report (Lang and Martinelli, 1979) some further applications of program AVALNCH are discussed, and the results of these applications are summarized.

Description of the AVALNCH computer code

The two-dimensional Navier-Stokes equations are the fluid-dynamics formulations which govern the problem. These equations, which are solved by finite-difference techniques, are:

(1)

where u and v are velocity components, g_x and g_y are gravity components, p =p/p, is the pressure/density ratio, and v is the kinematic viscosity. Additionally,

,

and

.

Each cell of the superimposed finite-difference grid (through which flow advances) is restricted by the incompressibility condition,

(2)

so that total mass is conserved, at least to within an error specified in the numerical algorithm. In the usual finite-difference approach the flow domain is partitioned into an array of uniform flow-domain cells surrounded by a single layer of boundary cells. For modeling the avalanche geometry, the cross-section of the avalanching snow-pack at the center line is confined to the flow domain.

The procedure adopted within the code is that, starting with an initial, stationary distribution of snow under gravitational loading, a pressure and velocity distribution are computed using Equation (1). The velocity and pressure thus computed will not balance to ensure zero divergence of mass, so that small adjustments using Equation (2) and linearized approximations of Equation (1) are imposed iteratively on the cells which contain snow. In this calculation, movement of mass between cells and into new cells may be necessary in order to achieve continuity. The initial calculation of velocity and pressure using Equation (1) is the beginning of a cycle (CYCLE ^* ) of calculations, and each adjustment of parameters in a sweep of the cells is an iteration (ITER). Clearly, the cost in this type of iterative analysis is related to the number of cells within the model. To standardize the runout problem and minimize cost, the array of flow-domain cells can be described as a “ribbon” along the avalanche path up to IBAR = 200 cells in length, and JBAR = 2 cells high. This array is then surrounded by a layer of boundary cells for a maximum array dimension of IMAX = 202 and JMAX = 4. The length of each cell is determined by dividing the avalanche path distance into 200 or fewer uniform lengths. Cell height is set equal to the maximum height of the snow-pack which releases in the starting zone. In many cases, average slab height will be less than cell height depending upon the distribution of the snow in the starting zone. Once the snow is in motion, subsequent flow height predicted by the program is the height of the core material, and not that of a possible snow-dust cloud. Most field observations of the flow height are to the top of the dust cloud. With these constraints, cell dimensions may range from 10 to 30 m in length and up to three or more meters in height. Thus, actual avalanche flow is confined to a height of one flow-domain cell, with a second cell allowed in order to account for any possible pile-up of the flow on an adverse slope. This occurs, for example, when an avalanche crosses a valley and starts to climb the opposite side. The profile layout that has been described is illustrated by the Ironton Park avalanche path (Ironton Quadrangle, Colorado) shown in Figure 1. The procedure is to read elevations and horizontal distances from a topographic map, plot the profile showing major slope variations and obstacles, and then lay-off cell lengths along the profile. In the case of Ironton Park, the path is represented by 110 cells, each 10 m in length. From this plot, elevation change (in meters) between successive cells is scaled and is one of the input arrays to the program. A second input array is the height of snow-pack in each cell that releases when flow initiates. For example, in the case of the Ironton Park path, snow-pack is specified in cells 1 through 9, with zero height in all remaining cells. The time to read the topographic map, set up the profile, lay-off the cells, and scale the elevation change between successive cells is on the order of one to three hours depending somewhat upon the complexities of the path.

Fig. 1. Longitudinal profile of Ironton Park aualanche path.

Two physical parameters that must be input are the surface friction fand the kinematic viscosity v. Surface friction is input either on a basis of one value per cell (third array of input data), or as a single value over the entire path. A physical basis exists for the input of surface friction as a variable along the path. For example, in coastal regions, dry snow in the starting zone may override dry snow until at some lower elevation the base snow becomes wet. Another case may involve snow running-out onto different base conditions such as rocky ground asphalt, etc. Thus, in the AVALNCH code, provision is made to input variable friction along the path. Currently, viscosity is treated as a constant over the entire path because data on any possible variation in v is lacking. However, if a user wishes v to be input as a variable, the coding changes to AVALNCH are straightforward.

The initial difficulty in applying AVALNCH to particular avalanche runs is in selecting appropriate values for f and v. By running a number of cases of recorded avalanches a range on v and f is established. This is discussed in detail in Part II (Lang and Martinelli, 1979).

The surface-friction models that are contained in the SOLA code are not adequate for avalanche modeling. In general, surface friction is controlled by specifying the velocity component parallel to the slope in the lower boundary cell u₁ to be a multiple of the corresponding velocity component in the flow-domain cell u₂ (Fig. 2). In AVALNCH the relationship is

(3)

which is equivalent to a surface velocity given by

For f= 0, u₁ = u₂ and the surface is slip-free. For f = 1.0, u₁ = —u₂ and u_s = 0 which is a no-slip condition. For intermediate values of f, the boundary is a partial-ship surface. Cases of f > 1.0 are valid and can be interpreted as a penetration of the surface roughness r into the flow. The surface where u _s is equivalently zero is defined by

(4)

as illustrated in Figure 2. However, there is an artificial aspect to this interpretation in that the quantity of fluid flowing is not reduced correspondingly for r > 0. Thus, the analogy which perhaps most closely represents the case f > 1.0 is that of sparsely-spaced rock protrusions or trees that slow the flow significantly, but do not trap large fractions of the flowing material.

Fig. 2. Effect of friction on boundary flow.

A condition which was observed in early computer runs was a tendency for the flow to reach an almost steady-state shallow flow not unlike equilibrium flow in open channels, and to runout for considerable distances. Since this phenomenon is not observed in most avalanches, a surface-friction law which is velocity dependent is incorporated in the code. At low flow speeds, the friction coefficient is increased exponentially with decrease in speed. The relationship established between f and velocity u is

(5)

where f ₀ is the initial value selected for f. This surface-friction-velocity relationship is written into the program so as to operate only after the avalanche has gone through its release, during which timef = f ₀. Although the quantitative basis for Equàtion (5) is unknown, this “fast-stop” option is a better representation for most avalanche flows than formulations which exclude it. Both friction and internal viscosity ought to be modeled as functions of both past and current states of flow, but in the absence of physical measurements to quantify this effect, recourse is made to the approximation described above.

In addition to the three arrays described previously which constitute part of the input to AVALNCH, another eight parameter values complete the input. These are:

• IMAX = number of flow-domain cells along path,

• JMAX = number of flow-domain cells normal to the slope,

• DELX = length of cell in meters,

• DELY = height of cells in meters,

• NU = kinematic viscosity (m² s^–1),

• FRK = surface friction value, if constant over path,

• TWFIN = time estimate for avalanche flow in seconds,

• CWPRT = number of cycles between extended printouts of flow profile.

These are listed in order to demonstrate the minimal amount of input data required to set up an avalanche run.

Output from AVALNCH includes two types of data depending upon user preference. The most complete output at the end of a cycle is a listing of the velocity, pressure, and height distributions of the flowing snow in each cell. The frequency of this printout is controlled by the parameter CWPRT, which designates the number of cycles between printouts of the complete flow distribution. At the end of cycles for which the complete printout is not specified, a single line of the following parameters is printed:

• CYCLE = current cycle number,

• ITER = number of iterations in current cycle,

• DELT = time increment of current cycle (in seconds),

• TIME = total time into flow (in seconds),

• FVOL = fluid volume (in m³),

• UMAX = maximum velocity parallel to the slope (in m s^–1),

• UEDG = leading-edge velocity (in m s^–1),

• LDEG = cell number occupied by the leading edge.

The snow profile can be plotted from the extended printout, as well as the velocity and pressure distribution in the moving mass. Among physical effects modeled by AVALNCH is the surface-wave phenomenon often observed in some avalanche flows. However, it is costly to output the extended printout repeatedly in a computer run, so that CWPRT is often set to a value larger than the number of CYCLES needed to complete a run. In this case an extended printout is output only on the final CYCLE of the run, and this defines the final distribution of snow debris.

The execution time taken by a CDC 6400 computer to complete an avalanche run varies between 0.5 and 3.0 min. Factors which influence the execution lime are the number of cells along the path, the magnitude of the friction and viscosity coefficients, and the geometry of the path.

We conclude this section by describing the physical conditions that an avalanche must exhibit in order to be modeled by this program. The Navier-Stokes equations are equations of motion for laminar fluid flow. If the internal viscosity is non-zero, then the flow is rotational. In laminar rotational flow, the particles of the fluid traverse curved paths or streamlines, and circulation is in evidence. In contrast, the particle motion in turbulent flow is completely random. From study of numerous avalanche films, our conclusion is that for the flowing motion avalanche (Perla, 1976), the flow is rotational and laminar. Values of Reynolds number substantiate this conclusion. In contrast, the motion of particles in an attendant snow-dust cloud may approach a turbulent condition because of reduced viscosity of the fluid. Results, which will be fully documented in a later report, show that the entire flow of the dense core material of an avalanche satisfies the Prandtl boundary-layer equations, which further substantiate the assumption of rotational laminar flow. Additional evidence supporting a laminar flow model is the necessity of a “fast-stop” option. Turbulent flow is characterized by dissipation which is proportional to velocity raised to some power, whereas the “fast-stop” requirement shows velocity dependence to a negative power. Thus, at low speeds a dissipation mechanism which does not show a behavior consistent with turbulent flow is operating. Whether this dissipation arises from increases in the surface friction or the internal viscosity is not known, and remains to be tested experimentally. In summary, our general conclusion is that most observable avalanches exhibit characteristics of rotational laminar flow, with turbulent or rigid-block cases being exceptional extremes. Details as to the type of dense core flow that produces the piston action for attendant dust-cloud flow of avalanches with a significant air-borne component are not known. However, based upon evidence from debris, the assumption of rotational laminar flow appears reasonable in many cases.

Ironton Park avalanche path

The Ironton Park path is located in the San Juan Mountains, Red Mountain Pass North; the south-east slope of Hayden Mountain. The path is of uniform width with runout onto a frozen lake bed. From topographic map readings, the profile of Figure 1 is constructed. In December 1958, a large, dry-snow avalanche crossed the road, noted in Figure 1, and left a debris pile approximately 1.5 m high on the road. No information is known on starting-zone fracture height, which must be considered exceptional in that an avalanche flowing to the road is a rare occurrence for this path.

The Ironton Park path was the first case considered in a check of AVALNCH. It was used to establish empirically the magnitude of viscosity v and surface friction f ₀ in modeling early winter, light, dry snow-pack. Snow-pack distribution in the release zone is assumed to span nine cells starting with cell 1. Cells 1 through 5 are assumed to be at full height, and cells 6 through 9 are assigned heights that correspond to a linear tapering decrease up to the leading edge of the flow. Later studies show that total length of the snow-pack in the release zone has negligible direct effect on runout distance, so that any reasonable modeling of the distribution of release-zone snow-pack is satisfactory. In all computer runs, the distribution of snow-pack is a linear scaling of the distribution described above.

A summary of the values of v and f ₀ used for different runout distances is given in Table I. Most of the runs are for a maximum release height of 2.0 m, which is considered an upper limit. As to snow type, the avalanche occurred in late December, which suggests a moderate mid-winter snow. On the basis of an initial slab height of 2.0 m, computer runs that characterize avalanche termination in the vicinity of the road are those for which the values of v+f ₀ are near to unity. For example, for v = f ₀ = 0.5, the avalanche leading edge was 60 m beyond the road, and the debris on the road was nominally 1.9 m deep. For v = 0.4 and f ₀ = 0.7, the leading edge of the avalanche reached the road, and for v = 0.5 and f ₀ = 0.6, the leading edge stopped 20 m short of the road. If v and f ₀ are further adjusted by small amounts, a closer correspondence between observed and computed results can be obtained. We further note from Table I the near correspondence of results when v and f ₀ values are interchanged. That is, for v = 0.4 and f ₀= 0.7, runout distance was approximately the same as for v = 0.7 and f ₀ = 0.4. Similar results were obtained for all other non-equal combinations of v and f ₀ reported in Table I, this indicates that these two parameters influence runout distance by almost equal amounts. Thus, in the present case we take v = f ₀ = 0.5 as representative of 2.0 m flow, and adapt the practice of selecting v and f ₀ to be equal for simple flow cases. Specific cases where exception to the v = f ₀ criterion is warranted are considered in Part II of this report (Lang and Martinelli, 1979). An equally valid result from Table I is the case v = f ₀ = 0.4 with an initial slab height of 1.5 m. As we lack specific information on fracture height, we cannot judge which case is the more accurate.

Figure 3 shows the total avalanche travel-distances obtained when the height of the initial-release slab is varied from 0.5 to 2.0 m, and v and f ₀ are both held constant at 0.5. These results show a non-linear dependence of travel distance versus slab height, the characteristics of which depend upon the local geometry of the path. All runs summarized in Table 1 are with the “fast-stop” option incorporated in the modeling. Without this option the results would be significantly different, particularly for the shallow-flow cases. Because of the lack of data on shallow-flow response, which would permit calibration of the “fast-stop” option for this type avalanche, AVALNCH must be considered as being primarily useful and accurate for the representation of moderate to large avalanches.

Fig. 3. Avalanche travel-distance for the Ironton Park avalanche as a function if initial slab height, with viscosity and friction cofficients held constant.

Table I. Computed runout for the ironton park avalanche obtained for different values of viscosity and Surface-Friction Coefficients

Summary

A standard computer code which models transient fluid-dynamics effects in two dimensions has been modified for specific application to the snow-avalanche-runout problem. Representation of avalanche flow as a transient phenomenon is an improvement over techniques developed previously. Modifications are made to the code which minimize computer running time, but which still retain accuracy by an empirical sizing of the material coefficients to match known avalanche responses. An important feature is the separation of geometric effects from surface-friction and internal-viscosity effects.

Results from a number of runs on a typical avalanche path with simple geometry show that viscosity and surface friction have approximately equal influence on avalanche runout Values of v and f ₀ in the range 0.4 to 0.6 yield flow conditions that resemble observed avalanche runout characteristics. It is shown that the height of a released snow-slab has a significant and non-linear influence on the avalanche travel distance. Quantitative evaluation of parameter sensitivity on avalanche flow is included in Part II (Lang and Martinelli, 1979).

Further empirical fitting of flow parameters to known avalanche cases will result in a rational refinement of parameter ranges, and this will aid program users by reducing the amount of time required for them to familiarize themselves with the program. AVALNCH can be used by those not familiar with computer languages and will even permit them to work out their initial avalanche case in the span of a day or two. For the user familiar with the Fortran IV computer language, the program is flexible in the sense that, as data which will permit a refined definition of avalanche response are accumulated, the code can be updated to accommodate the changes. The most pressing need in the development of more accurate applications of AVALNCH to different runout situations is for more physical data to permit a matching of program output to actual results.