DYNAMIC MODELLING OF E X T R E M E SPEED PROFILES OF D R Y FLOWING A V A L A N C H E S by Christopher Paul Borstad B.Sc. Colorado State University, 2002 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Applied Science in The Faculty of Graduate Studies (Civil Engineering) The University of British Columbia December 2005 © Christopher Paul Borstad 2005 11 Abstract A brief review of different approaches to avalanche dynamics modelling is considered, drawing on observations and knowledge of avalanche behaviour and analogies with other physical phenomena. The numerical model used for this analysis, DAN, is described with respect to the formulation of the equations of momentum and continuity. Avalanche path geometry is approximated as a rectangular channel with the possibility of variable width. The flow mass is discretized longitudinally into a variable number of constant-volume blocks. The uncertainty surrounding the dynamics of frontal ploughing of undisturbed snow, basal erosion, and mass entrainment is discussed. In order to avoid this uncertainty, which is implicit when modelling a slab avalanche beginning from rest in the starting zone, the calculations in DAN begin near the middle of the slope where the flow is expected to achieve maximum speed. Knowledge of the maximum probable speed for an extreme avalanche is used to supply an initial, non-zero speed to the flow near the middle of the path, and entrainment is assumed negligible as the flow decelerates. A sensitivity analysis explores theoretical considerations surrounding the modelling of a flow mass in multiple dimensions, given this new procedure as well as traditional from-rest approaches. The model is shown to be more sensitive to changes inflowwidth than flow length or depth. DAN calculations are compared against examples of recorded extreme avalanche data, with generally good agreement. The new modelling procedure was able to reproduce the observed sharp deceleration of avalanche frontal speed in runout zones. This is worthy of remark, because most decisions about land use planning and hazard assessment are made in the lower reaches of avalanche paths. iii Contents Abstract ii Contents iii List of Tables v List of Figures vi Acknowledgements viii 1 1 1 3 3 6 8 10 13 14 15 A Brief Outline of Avalanche Dynamics Models 1.1 Introduction 1.2 Physical Description of Dense Flowing Avalanches 1.2.1 Observations and Measurements of Dense Flowing Snow 1.2.2 Frontal Ploughing, Basal Erosion, and Mass Entrainment 1.3 Relevant Granular Experiments and Models 1.4 Snow-Specific Models 1.4.1 Friction Parameters 1.4.2 Internal Pressure 1.5 Discussion 2 Numerical M o d e l Description and Features 2.1 The New Modelling Procedure 17 21 3 Theoretical Considerations on a Smooth Parabolic Path 3.1 Changing Path Width 3.2 Internal Pressure Formulation 3.3 Shape of Initial Flow Profile 3.4 Dependence of Flow Depth on Runout 3.5 Dependence of Flow Length on Runout 3.6 Variable Number of Mass Flow Blocks 3.7 Variable Initial Speed Profile 3.8 Discussion 23 24 25 29 30 38 41 43 44 4 Model and Procedure Validation 4.1 Validation of DAN With Simulations Starting From Rest 4.1.1 Full-Simulation Discussion 4.2 Validation of the Proposed Method 4.2.1 The Aulta avalanche path and measured events . 4.2.2 The Ariefa-Samedan avalanche of January, 1951 4.3 Examples of DAN Calculations for Surveyed Avalanche Paths 4.4 Discussion . 46 47 53 54 56 58 62 64 . Contents iv 5 Conclusions 5.1 Recommendations 66 68 Bibliography 71 List of Tables 3.1 Calculated runout distance for the initialflowprofiles in Figure 3.7 vi List of Figures 1.1 1.2 Definition of geometric terrain parameters for empirical runout calculations Typical profile of slab avalanche initiation and subsequent granular flow 2.1 2.2 Profile view of the flow discretization in DAN Geometry and forces acting on a boundary slice in DAN 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15 3.16 3.17 3.18 Frontal speed vs. path width for constant initial flow depth Runout Ratio vs. runout zone spreading ratio Runout Ratio vs. variable or constant internal pressure coefficient Variation in the internal pressure coefficient with position Flow and deposit length vs. variable and constant internal pressure coefficient . . . . Mean flow and deposition depth vs. variable and constant internal pressure coefficient Possible initial flow profiles for a given volume Variable initial depth for a given starting position and initial frontal speed Flow thinning under gravity and associated frontal and tail speeds Average and maximum flow depths vs. initial depth Speed profiles for variable starting positions and given maximum speed Runout Ratio vs. maximum and average flow depth Deposit depth vs. flow depth Runout Ratio vs. scaled flow depth Runout Ratio vs. scaled flow length Scaled deposit length vs. Scaled flow length Runout Ratio vs. scaled flow length . Speed profiles and runout distances for a givenflowmass divided into different numbers of mass flow blocks 3.19 Deposit profile for the same total mass divided into different numbers of mass flow blocks 3.20 Runout for variable initial speed profiles 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 Comparison of measured and DAN frontal speeds for the Arabba avalanche of 21 December, 1997 Variation in runout distance with changing friction coefficient fi for the Arabba avalanche path Comparison of measured and DAN frontal speeds for the Madergrond avalanche of 17 January, 1985 Comparison of DAN versus the model of McClung (1990) for the Madergrond path . Comparison of measured and DAN frontal speeds for the Khibins (a) avalanche . . . Comparison of measured and computed frontal speeds for the Aulta avalanche of 8 February, 1984 Comparison of DAN with the model of McClung (1990) for the Aulta path Speed profiles for different starting positions for the Madergrond avalanche path . . 2 7 18 19 24 25 26 27 28 28 29 31 33 34 35 36 37 38 39 39 40 41 42 44 48 49 50 51 52 53 54 55 List of Figures vii 4.9 DAN runout vs. flow length and mean depth for the Aulta avalanche path 4.10 Comparison of measured and DAN speeds, starting in the middle of the path, for the Aulta avalanche of 8 February, 1984 4.11 Comparison of DAN starting in the middle of the path with the model of McClung (1990) for the Aulta path 4.12 Speed profiles for the Ariefa-Samedan avalanche calculated using DAN and a Voellmyfluid model 4.13 DAN frontal speed profiles for the Ariefa-Samedan avalanche of January, 1951, starting from rest and in the middle of the path with an input initial speed 4.14 Runout vs. flow length and mean flow depth for the Granite Creek 13 avalanche path in the North Cascades, Washington, USA 4.15 Upper- and lower-bound frontal speed profiles for the Granite Creek 13 data in Figure 4.14 4.16 Upper- and lower-bound frontal speed profiles for the Granite Creek 7 avalanche path 57 58 59 60 61 62 63 64 Vlll Acknowledgements I am grateful to the Natural Sciences and Engineering Research Council of Canada for the support of this research. Oldrich Hungr in the Earth and Ocean Sciences Department at the University of British Columbia provided the numerical model for this work, and his assistance is greatly appreciated. David McClung of the Department of Geography and Greg Lawrence of the Department of Civil Engineering at the University of British Columbia played invaluable advisory roles. 1 Chapter 1 A Brief Outline of Avalanche Dynamics Models 1.1 Introduction Human population expansion into mountainous areas has necessitated the study of flowing avalanche dynamics. Responsible land-use planning relies on hazard mapping techniques that require, as input, predictions of the stopping position and speed profile of an extreme avalanche. The placement and defense of civil infrastructure, residential buildings, and recreational facilities may depend on some estimate of avalanche speed and runout distance. In addition to the runout distance of an extreme event, the destructive potential, as a function of position, is important. The design of defense structures in avalanche terrain relies on the likely speed of an extreme event, coupled with the flow density, to provide impact pressure estimates. These necessary speed calculations are one of the primary outputs of a dynamics model. Dynamic modelling involves the selection of resistance parameters to calculateflowspeed and runout distance. As speed is a function of position, so is the associated impact pressure. Drawing appropriate boundaries between acceptable and prohibited locations for the placement of facilities then involves the determination of acceptable levels, if any, of potential impact pressure. A dynamics model provides one method for calculating runout distance; another involves empirical modelling. Empirical models use a data set of local extreme runout data to statistically predict the runout for a given avalanche path. This is generally done using either a least-squares method or an extreme-value distribution, fit to the data set. For a recent review and comparison of empirical models, see McClung (2001a). Figure 1.1 contains definitions of geometric parameters for empirical runout calculations, some of which have been adopted in this study. 2 Figure 1.1: Definition of geometric terrain parameters for empirical runout calculations. The alpha (a) angle is defined as the angle sighted from the extreme runout position of an avalanche path to the top of the starting zone. The beta point is the position where the slope angle first decreases to 1(T. The beta angle is the angle sighted from the beta point to the top of the starting zone. The Runout Ratio, as defined in McClung and Mears (1991), is the ratio of the horizontal distance from the beta point to the extreme runout position divided by the horizontal distance from the starting position to the beta point, or jj£. The total vertical fall of the path is denoted as H. An advantage of empirical runout prediction over dynamic modelling is the ability to quantify uncertainty. The variance of the data set used is crucial, however. Given that extreme events are rare, and that data from one mountain range cannot necessarily be extrapolated to another, it is often difficult to compile data sets large enough for statistical analysis. Empirical methods are built entirely on field data, thus they are verifiable. However, such methods are inappropriate for predicting runout on slopes that have irregular geometry. As discussed in McClung (2001a), it may not be possible to uniquely identify terrain parameters on some slopes. Additionally, empirical methods do not make any prediction of avalanche speed. If speed estimates are needed for impact pressure calculations, empirical methods of runout prediction need to be complemented with dynamic modelling. In applications, an empirical method can be used to specify an avalanche stopping position, and a dynamics model can then be used to calculate the avalanche speed, as recommended by McClung (1990). 3 1.2 Physical Description of Dense Flowing Avalanches Dynamic modelling follows from the pursuit of an accurate physical description of avalanche motion. Figure 1.2 shows the general character of a dry slab avalanche from the time that the snow slab is released to the stage of developed granular flow. Dense snow avalanche flow is an unstable, three-dimensional phenomenon with likely nonlinear rheology (Ancey and Meunier, 2004). Though a complete physical representation is in advance of theory and detailed observations, certain bulk characteristics can be deduced. Dynamic modelling necessarily begins with some knowledge of the physics of flowing snow. This knowledge comes from a combination of field experience, observations, experiments, and analogies with other earth processes. Indeed, dynamics models have historically been based on analogies with other physical phenomena (Ancey and Meunier, 2004). If dense flowing snow is considered a granular flow, then granular and continuum modelling approaches are both possible (Goodwin and Cowin, 1971). The granular approach describes a collection of individual finite particles (e.g. Dent, 1992; Hutter et al., 1987; Jenkins and Savage, 1983). The continuum approach represents material properties as continuous functions. In this manner, the flow can be divided arbitrarily without loss of generality. Both approaches rely on knowledge of the physical character of flowing snow. 1.2.1 Observations and Measurements of Dense Flowing Snow The extreme avalanches of interest here are generally released as cohesive slabs of snow. Extreme avalanches are defined as events of large magnitude relative to the size of the given avalanche path. Such avalanches occur with low frequency, or long return period. Depending on the context, the return period of a design event may range from as little as 30 years up to 100 years or more. A rapidly propagating shear fracture in a weak layer beneath a slab followed by a tensile fracture within the slab release a large volume of snow simultaneously. As this slab begins to slide downslope, it rapidly accelerates and breaks into many blocks of different sizes. If the initial snow slab is weak, most of the sliding blocks are broken up into small fragments as the avalanche flows downslope. These flows become 'fluidized,' with small snow and ice particles making up most of the mass. Conversely, a strong initial slab may not become entirely fluidized, and blocks on the order of 10-10 cm across may remain present on the surface of the deposit. In either 2 case, the resulting granular flow is composed of ice crystals, snow blocks, possibly liquid water, and 4 interstitial air. Internal shearing of granular flows is only possible if the slope angle is greater than the internal friction angle of the granular material. Tiefenbacher and Kern (2004) used a large chute, sloped at 32°, to measure the properties offlowingsnow. They observed an approximately 5cm basal shear layer with roughly constant velocity (plug flow) above this layer. This suggests that the internal friction angle of the snow was near to or greater than 32°. It is not clear, however, how the results from these small, experimental slides might scale to typical avalanche volumes. For extreme avalanches, the applicability of chute flow data is questionable. McClung and Schaerer (1985) found fluctuation in pressure records for dry snow avalanches on a 31° slope, suggesting that the snow was not yet in a locked state. This would indicate an internal friction angle of less than 31°. In any case, descriptions of plug flow have been made by others for snow avalanches and granular flows, especially in the runout zone when the slope angle decreases (e.g. Hutter et al., 1987; McClung and Mears, 1995; Bartelt et al., 1999; Louge and Keast, 2001). Schaerer and Salway (1980) observed dense flow in the bottom of an avalanche channel, an overlying zone of 'light,' or intermediate, flow, and an enveloping turbulent powder cloud. In pressure records, the dense and light parts of the flow were both observed to appear in waves rather than as continuous flow. The mid zone of light flow (1-5 m thick) was observed to overrun the dense part in the runout zone when the front of the dense flow came to rest. For this reason, the intermediate flow layer was claimed to be different than a saltation layer. Schaer and Issler (2001) confirm the presence of a similar layer above the dense flow, with particle sizes in the range 1-50 cm; however, they retain the phrase saltation layer to describe this component of the flow. McClung and Schaerer (1983) measured avalanche speeds in the middle of avalanche paths at Rogers Pass, British Columbia, Canada. They reported values in the range 12-62 m/s for the middle of the 38 paths observed. McClung and Schaerer (1993, p. 110) describe an upper limit envelope for maximum avalanche speed (hereafter referred to as ULE), based on measured speeds scaled with the total vertical fall of the respective avalanche path. For data at Rogers Pass, the U L E can be expressed as W = 1&>/H, where H is the total vertical fall of the path in metres, v max (1.1) is in m/s, and the constant 1.8 necessarily has units y/m/s. Alternatively, McClung (1990) described an upper limit envelope scaled with 5 the curvilinear path length, rather than the vertical fall, for avalanche speed data from Canada, Switzerland, and Norway. It may be more appropriate to express a maximum speed envelope as (1.2) where c is now a nondimensional constant and g is gravitational acceleration. In this manner, such an expression could be compared to analogous relations for other types of gravity currents. Equation 1.1 would then take the form (1.3) where c = 1.8/y/g = 0.6. Drawing on a fluid mechanics analogy, some authors have reported the Froude number (Fr) for flowing avalanches, defined as where v is flow velocity, g is gravitational acceleration, and L is a characteristic length. The Froude number is a measure of the prevalence of gravity waves in a flowing fluid (Wilcox, 2003). Values of Fr > 1 indicate supercritical flow; disturbances experienced by the flow front cannot propagate upstream. This result was postulated theoretically by Salm (1966) and Mellor (1968). Values of Fr for flowing avalanches have been cited in the range 6.2-8.5 (Schaerer and Salway, 1980), 4-12 for the experimental flows of Tiefenbacher and Kern (2004), and as high as 15 (Naaim-Bouvet et al., 2004). Visual observation of the dense, destructive core of a dry snow avalanche in motion is difficult due to the usual presence of an enveloping turbulent powder cloud. For this reason, impact pressure measurements have provided valuable insight into the nature dense flowing snow. McClung and Schaerer (1985) report peak impact pressures on the order of 10-100 kPa for speeds exceeding 30 m/s. They noted that dry avalanches have higher peak impact pressures, but lower average pressures, than wet avalanches. From measured speeds and duration of impact records, flow lengths for the reported events lie in the range 10-600 m. Higher peak impact pressures (on the order of 100 kPa) were correlated with the longest flow lengths (on the order of 100 m). Mean flow density was reported as 215 kg/m , and an upper bound on the effective flow density, defined as the peak impact pressure 3 P divided by the square of the velocity, could be as high as 700 kg/m . The volume fraction of solids 3 m 6 for the avalanches recorded was as high as 55%. These are important results for the formulation of dynamic modelling parameters, to be discussed below. 1.2.2 Frontal Ploughing, Basal Erosion, and Mass Entrainment As a flowing avalanche moves downslope, it encounters undisturbed snowcover ahead of the flow front. This snow can be ploughed above or below the advancing flow front, to be entrained into the dense core or turbulent powder cloud. The body of the avalanche, behind the front of the flow, can erode stationary basal snow and entrain this mass into the flow. The processes of ploughing, erosion and subsequent entrainment of existing snowcover are not well understood, however. Early dynamic modelling efforts neglected entrainment entirely, setting the mass of the flow constant (e.g. Voellmy, 1955; Perla et al., 1980; Ancey and Meunier, 2004). This assumption or simplification is still widespread. Schaerer (1974) claimed that entrainment made a negligible contribution to avalanches observed at Rogers Pass, British Columbia, Canada. However, more recent field observations have indicated that entrainment can add considerable mass to the dense flow from initiation to deposition. Mass balance measurements of four avalanches in the 1997-1998 winter in Switzerland (Sovilla et al., 2001) show a large variation between released mass and maximum flowing mass. Gauer and Issler (2004) illustrate several possible mechanisms of erosion of snowcover. At present, insufficient data exist to validate any of the numerous mechanisms of erosion or entrainment proposed by various authors. However, the likely importance of these processes on the force balance suggests that neglecting them is inappropriate. The uncertainty propagated by the assumption of no mass entrainment cannot currently be quantified. For extreme avalanche prediction, McClung (1990) argues that ploughing, entrainment, and deposition are implicitly included in his model. If a model is calibrated with data from extreme events, which are necessarily characterized by minimum resistance, then inertial effects such as entrainment, which would add considerable resistance to the flow, are implicitly accounted for without explicitly specifying them in the boundary conditions. Given the inappropriateness of neglecting entrainment for short return period avalanches, it can be argued that, at present, extreme events with minimal entrainment are the only ones that we can model. 7 Start of motion. Granular flow begins. Granular flow obscured by suspended material. Figure 1.2: Typical profile of slab avalanche initiation and subsequent granular flow (McClung and Schaerer, 1993, Figure 5.24). In the top frame, a cohesive snow slab is released all at once after a shear fracture in a weak layer in the snowpack followed by a tensile fracture in the snow slab. After a short time, the snow slab breaks up into many snow blocks of different sizes (middle frame). Snow blocks are further broken into small ice particles by shearing and internal deformation of the dense granular flow. 8 1.3 Relevant Granular Experiments and Models Granular experiments have provided an analogous basis for analyzing dense flowing snow. The parallels between granular chuteflowsand dense snow avalanche flow are numerous, which is advantageous for a number of reasons. First, granular experiments are relatively easy to carry out and observe compared to snow chute flows or experimental observation of natural avalanches. Additionally, the literature on granularflowsis extensive, covering a multitude of grain/fluid mixtures, grain sizes and types, and boundary conditions. Care should be exercised, however, in adopting granular theories for dense flowing snow. The base of a dense avalanche is composed of small, irregularly shaped ice crystals with variable physical properties. Nonetheless, homogeneous granular experiments provide a good basis upon which to build equations for flowing snow. The groundbreaking work of Bagnold (1954) measured shearing of suspended spherical particles in an annular shear cell. The grain concentration was varied from 13-62%. Bagnold identified two flow regimes: an inertial regime for high shear rates, and a viscous regime for low shear rates. The shear speeds relevant for dense flowing snow fall into the inertial regime, for which both the shear and normal stresses are proportional to the square of the shear strain. For snow avalanches, the density of the snow particles is on the order of 10 to 10 times more 2 3 dense than the interstitial fluid (air), in contrast to Bagnold's neutrally buoyant, fluid-suspended particles. Another difference between Bagnold's results and those expected for snow avalanches is in the nature of the grain collisions. For flowing snow, the collisions between particles are destructive (McClung, 2001b). These results imply that collisions between snow particles in avalanche flow must be the dominant mechanism for momentum transport. An important result discussed by Bagnold is the suppression of turbulence in the interstitial fluid for high grain concentrations. He defined a nondimensional parameter, now called the Bagnold number, which indicates the dominance of either viscous stresses or grain collision stresses in the flow. The inertial regime exists for Bagnold numbers of greater than about 400. Snow avalanche flows are characterized by Bagnold numbers on the order of 10 to 10 (McClung, 2001b). Given the 3 4 volume fraction of solids inferred by McClung and Schaerer (1985) for dense snow avalanches (up to 55%), turbulence in the interstitial air should similarly be suppressed. Savage and Jeffrey (1981) analysed the granular flow of smooth, hard, elastic spheres at high 9 shear rates. They restricted their analysis to cases in which the interstitial fluid plays a negligible role in the force balance. Their analysis is similar to a statistical-mechanical description of the liquid phase, with momentum transport dominated by particle collisions. They described three general sources of bulk stress: 1. dry friction/rubbing 2. momentum transport by particle motion 3. momentum transport by particle interactions In most cases, they assert, one of these processes will be dominant. For the high volume fraction of solids in their analysis, momentum transfer is dominated by interparticle collisions, and frictional rubbing is negligible. Additionally, in agreement with the experimental work by Bagnold (1954), they found that the shear and normal forces were coupled with the square of the shear rate. Jenkins and Savage (1983) found the same coupling of shear and normal forces for granular flow of identical, smooth, nearly elastic spheres. They assert that the ratio of shear to normal forces, however, depends on the boundary conditions. The coupling between the shear and normal stresses is dependent on the basal or free surface boundary condition. If some flux of energy is specified across the boundaries, the ratio of shear to normal forces is modified, but the coupling is preserved. For the case in which there is no boundary energy flux, the shear and normal forces are proportional to the square of the strain rate. The specification of boundary conditions thus becomes important in the formulation of equations of motion for granular flow and analogous avalanche dynamics models. Jenkins and Savage (1983) report that there is always some basal slipping, even over a roughened boundary, for granular flows. Louge and Keast (2001) assert that the nature of a granular flow is strongly dependent on the base over which it flows. Flat 'frictional' sliding surfaces produce different characteristic flows than inclines with a roughness height on the order of the size of the particles. On planar inclines, they observed a basal shear layer with overlying passive plug flow. The shear layer thickness was on the order of the diameter of the particles. The grains interact by a combination of 'enduring' frictional contacts and instantaneous collisional interactions. For steady, fully developed flow, the shear and 10 normal forces are coupled with the tangent of the slope angle via S/N = tan i>} (1.5) With a theoretical argument that the normal force at the base of the flow cannot be dissipated entirely by collisional forces, Louge and Keast (2001) write the total basal stress as the simple sum of an 'impulsive' and an 'enduring' term. This formulation preserves the observed coupling between the shear and normal forces, and allows the dominant mechanism of momentum transport to shift between particle collisions for high Froude number flow and sliding friction (at high enough solids fraction) for low Froude numbers. 1.4 Snow-Specific Models The most popular avalanche dynamics model to date was introduced by Voellmy (1955). He modelled flowing snow as an 'endless' fluid using open channel hydraulic theory. In the centre of an infinite, constant slope channel, a force balance was written that lumps the resistive terms at the base of the flow. This basal resistance was assumed to take the form of a simple sum of a Coulomb-like sliding friction and a turbulent resistance proportional to the square of the speed. In this formulation, there is no coupling between the shear and normal forces in the flow. Given values for the two resistance parameters and a flow depth, the terminal velocity of the flow is predicted by the equation V = [SH(sinB - ncosO)] ' 1 2 t (1.6) where £ is a turbulence coefficient, H is flow depth, 6 is the slope angle, and / i is the coefficient of sliding friction. The Voellmy formulation makes no prediction of where on the slope this maximum velocity occurs, however. If the position at which deceleration commences can be specified, then the runout distance is an analytical solution to Equation (1.6): 8 = V? 1* ^—, (1.7) [2g(u-cos9 - sin6) + jf£] 1 1.1 Louge and Keast (2001) use a for the slope angle. V is used here to avoid ambiguity with the a-angle of Figure 11 where Hp, the mean deposition depth, must be assumed. The need to specify the location at which deceleration begins potentially introduces subjective human judgement into the runout calculation. This undesirable consequence was addressed by Perla et al. (1980). Building upon the work of Voellmy, they derived a depth-averaged, centre-of-mass formulation for calculating maximum speed and runout. The equation for maximum avalanche speed on a constant slope with constant friction coefficients and tan 8 > fj. (so that the gravitational driving force is greater than the Coulomb frictional resistance in order to produce an acceleration) is given by V = {^-{sind - ncosd)] ' . (1.8) 1 2 t The resistance term Mg/D is equal to £H in the Voellmy formulation, and proportional to v . It is 2 assumed to have three components: a centripetal force, an entrainment term, and a drag or ploughing term. By dividing up an avalanche path into discrete segments of different slope angle, equation 1.8 can be applied to calculate the speed at the end of a segment. A momentum correction at changes in slope retains only the slope-parallel component of velocity for the new slope segment. In this manner the solution can be iterated downslope. If the velocity drops to zero before the end of a slope segment, the runout distance for that segment is interpolated. In this manner, it was no longer necessary to specify, or guess, the position at which the avalanche would begin to decelerate. A more complex continuum derivation was made by Norem et al. (1987) using a CriminaleEricksen-Filbey fluid analogy. As opposed to most previous models, this approach attempts to account for internal forces in the body of the avalanche. The normal stress in the flow was divided into terms describing effective pressure, pore pressure (which, for dry snow, should be negligible), and dispersive grain pressure. The shear stress terms included cohesion, Coulomb friction, and dynamic stresses. Though this approach is physically more descriptive, it requires the specification of numerous physical parameters that are unknown for either static or flowing snow. McClung (1990) emphasized that the observed rapid deceleration of avalanche flow in the runout zone results in the need for high precision data on runout distance and speed in avalanche dynamics models. McClung proposed scaling the speed profile of an avalanche by initially specifying the runout distance. This effectively divides the problem into two separate calculations: runout prediction and speed calculation. 12 Given the difficulty in accurate flow-depth predictions, and the uncertainty in the dynamics of internal deformation of flowing snow, the proposed approach of McClung (1990) is independent of depth. In the absence of entrainment, rough depth approximations can be made given initial slab fracture limits and path geometry. However, given the expected importance of entrainment in the mass balance and dynamics of flowing avalanches, excessive sensitivity to flow depth would be an undesirable feature of a dynamics model. McClung and Mears (1995) outline an approach that calculates the position of the leading edge, rather than the centre of mass, of an avalanche. Such an approach leads to more conservative estimates of runout as well as runup against adverse slopes. The formulation explicitly includes a passive internal pressure calculation that is a function of the local slope angle. A centre-of-mass model, or any other one-dimensional formulation, can only implicitly account for internal processes. A granular model for dense flowing snow was proposed by Hutter et al. (1987). However, they found the granular approach to be less helpful than expected. Numerous parameters in the analysis were difficult to specify, such as the diameter of particles, coefficient of restitution, fluctuation energy, and a granular temperature. The thickness of the fluidized basal shear layer must also be specified, and the results are sensitive to this uncertain input. A simplified granular approach made by Dent (1992) considered the theoretical shear flow of rough, inelastic spheres between parallel rough plates. Newton's laws of motion were applied to determine the two components of linear acceleration and the single component of angular acceleration for each particle at each time step. After numerical integration, the position and velocity of each particle were updated. An interparticle friction coefficient n described the rubbing friction between particles in enduring contact. Even when n = 0, there was significant dynamic friction due to momentum transfer by normal collisions between particles. The coupling between the shear and normal forces was still observed. At the highest rates of shear, the dynamic friction approached the value of the interparticle friction coefficient ft. Additionally, the dynamic friction was found to be independent of the thickness of the shear layer. According to Dent (1992), increasing normal stress leads to decreasing dynamic friction, defined as the ratio of shear to normal stresses. Thus, friction should be lower at the base of the flow, which is deeper and thus has higher overburden. This may explain the presence of a thin, localized basal shear layer for some avalanches. Additionally, increased shear speed corresponded to increased 13 dynamic friction. 1.4.1 Friction Parameters Many different frictional formulations and material rheologies have been proposed to describe dense flowing snow. The models discussed above all require some basic assumptions about the form of the different resistance mechanisms in avalanche flow. It is generally agreed that the majority of resistance is encountered at the base of the flow. Salm (1966) proposed that a general formulation for the basal resistance could be written as a sum of terms: R = o + a\v + a,2V , 2 0 (1-9) plus higher-order terms. The first term, independent of speed, is often taken as a Coulomb-type sliding friction. The second term is analogous to a linear viscous resistance, and the third term is an inertial term, often used to describe turbulent fluid resistance or centripetal effects. The variables a are assumed to be constant; however, in reality they may have a functional variation with speed or other physical parameters. Salm (1966) and Brugnot and Pochat (1981) use all three parameters of this type in formulating equations to describe dense flowing snow. Perla et al. (1980) point out, however, a mathematical redundancy in any multi-parameter model. An infinite number of parameter combinations can reproduce a given runout position. If additional data are available to supplement runout distance, such as measured speeds, the selection of two or more parameters is better confined, but the redundancy remains. Dent and Lang (1980) modeled snow avalanches as modified Bingham fluids, with viscous flow below a given yield stress. Lang et al. (1985) used a Bingham-type model for a basal shear layer with an overlying passive load, and introduced a 'fast-stop' exponential decay of velocity in the runout zone in order to reproduce the steep gradient of velocity observed as the flow decelerates. These analyses are based on ad-hoc proposals, with interesting qualitative behaviour but little effort on verification. Using a granular approach similar to that of Dent (1992), Salm (1993) derived the familiar Voellmy parameters ix and £ for the hypothetical shearing of identical, spherical snow blocks. Salm 14 assumed that avalanche flow consists of snow blocks of size on the order of 10cm, not individual ice particles. McClung (2001b) shows that the size of snow and ice particles in avalanche flow is more likely to be on the order of a few millimeters. Salm set the coefficient of restitution to be zero for these colliding snow clods, in contrast to Dent. However, as discussed above, observations indicate concentrated shearing near the base of granular flows, and field observations of avalanche debris indicate decreasing particle size with depth below the surface. Thus, it would seem likely that the majority of momentum transfer in dense avalanche flow is by the granular shearing of small, individual ice particles near the base of the flow. Furthermore, the coefficient of restitution for small ice particles is not likely zero. McClung (1990) has argued that the majority of the basal resistance is in the form of a Coulomb sliding friction. He related the basal sliding friction coefficient \x to both the cv-angle of the avalanche path and flow speed. If surface drag can be neglected, then the speed dependence of the friction is dropped, and \i is constant and equal to tan a. In this case, the friction coefficient is interpreted as the mean value of the total flow resistance from flow initiation to deposition. Ancey and Meunier (2004) underline the interpretation of fi as a bulk parameter that accounts for all internal resistive terms in the flow as well as the boundary conditions. They back-calculated values for this parameter based on observed events, and found / i to be non-linear in speed. Tiefenbacher and Kern (2004) calculated a bulk dynamic basal friction coefficient of 0.72 ± 0 . 1 1 from chute flow measurements. They claim this value to be high, compared to the range of previously used basal friction coefficients, because it takes into account internal resistance. 1.4.2 Internal Pressure Internal energy considerations are difficult to account for in the equations of motion for dense flowing snow. This is due to the difficulty in measuring and observing dense snow avalanches. Simplified centre-of-mass or other one-dimensional models cannot explicitly account for internal flow processes (McClung and Mears, 1995). Additionally, the constitutive equation for snow is unknown, making many modeling parameters, especially an internal pressure, particularly speculative. In fluid mechanics, the simplest form of internal pressure is hydrostatic. Hydrostatic pressure accounts for stresses when the strain rate is zero. For snow avalanche flow, the most widely used modification of the simple hydrostatic pressure formulation draws on a physical analogy from soil 15 mechanics. Rankine earth pressure theory describes 'active' and 'passive' states of static soil pressure. The active state is analogously linked to extending granular flow, while the passive state represents compressing flow. These states are manifested in a multiplicative coefficient of the given hydrostatic pressure equation. This concept was first applied to snow avalanche equations by Salm (1966), and has been adopted widely since then. The justification of the analogy between a static earth process and a dynamic granular-fluid flow has not been rigorously defended, however. Salm (1993) used a constant passive internal pressure in equations for calculating runout. However, he found agreement with existing data only for very large values of the passive pressure coefficient (A). Assuming a Voellmy resistance law, values of A = 10 and 20 provided better fits to observational data. In Salm's formulation of A, these values correspond to internal friction angles (<}>) of 55° and 65°, respectively. The more traditional value of A = 2.5 (<j) = 25°) was discarded as 'too small.' McClung and Mears (1995) discuss static values of <j> for snow in the range 50-60°, but for flowing snow they assert that these values are too high. Other investigators have adopted a pressure term that varies between active and passive states, depending on different flow conditions (e.g. Hungr, 1995; Sartoris and Bartelt, 2000). T6masson and Johannesson (1995) switch between passive pressure for compressive flow and zero internal pressure for extending flow. In their interpretation, passive internal pressure opposes longitudinal flow compression, but for extending flow the internal pressure plays no role. The granular theory of Louge and Keast (2001) entirely ignores the earth pressure analogy for longitudinal straining, and the resultant equations of motion nonetheless reproduce observed granular flow features. 1.5 Discussion The difficulty in measuring the properties of flowing snow has led to a great deal of speculation in the formulation of dynamics models. Necessarily, simplifications must be made when modeling a nonlinear, three-dimensional, unsteady phenomenon. Lumping unknown internal resistance terms into a bulk basal resistance formulation has been commonplace. Attempts at explicitly specifying internal avalanche dynamics, while qualitatively more physically descriptive, are faced with difficulty in thecalibration and selection of multiple uncertain parameters. The divergence in basal frictional formulations between different investigators is evidence of both 16 the uncertainty in flow dynamics and differing opinions on how to best simplify a complex process and yet produce meaningful and reliable results. The early experimental work of Bagnold (1954) suggests that turbulence of the interstitial fluid in a granular suspension will be suppressed given high enough volume fraction of solids; for dense flowing avalanches, this finding is expected to hold. Nonetheless, turbulent resistance parameters continue to be used in many dynamics models. Other speed-dependent resistance terms can describe centripetal effects or surface air drag, but these can often be neglected due to their small contribution to the force balance compared to a Coulomb-type basal sliding friction. Salm (1993) re-interprets the familiar Voellmy resistance parameter £ to account for collisional interactions between flowing snow granules. Recent advances in granular flow theory, however, (e.g. Louge and Keast, 2001) account for collisional interactions between individual particles without violating the observed coupling between shear and normal forces at the base of the flow. Ploughing and entrainment of snowcover ahead of a flowing avalanche remain largely uncertain processes. These effects have recently been shown to play an important role in avalanche mass balance for short return period events (Sovilla et al., 2001), and their effect onflowdynamics is expected to be significant. The common practice of neglecting these processes from avalanche initiation to deposition is questionable. However, attempts at explicitly accounting for ploughing, erosion, or entrainment are largely speculative and in advance of the necessary, but difficult, field measurements. Clearly the field of avalanche dynamics modelling has made progress in the last 50 years. Improved observational data and methods have provided valuable foundations upon which to build numerical models. A good deal of work remains in order to reduce, and especially manage, the uncertainty encountered in modelling such a complex phenomenon. Much of this uncertainty cannot be quantified, and often it is not discussed. The work that follows is an attempt to better manage the uncertainty encountered in avalanche dynamics modelling. The primary means of achieving this will be through limiting the domain of numerical modelling to avalanche runout zones. It is assumed that entrainment effects are small when the flow is decelerating compared to when the flow is accelerating in the starting zone and track. Correspondingly, the assumption that the flow mass will be constant in the runout zone for extreme avalanches is more reasonable than holding the mass constant over the entire avalanche path. 17 Chapter 2 Numerical Model Description and Features The numerical model ("DAN") used for the following analysis was developed for the unsteady flow of earth masses. The model was proposed and outlined by Hungr (1995). Though primarily intended for landslides, DAN was built to allow the option of selecting a variety of material rheologies (e.g. Voellmy, Bingham, and frictional) for the flow mass in question. Given the assumption that the majority of flow resistance in dense snow avalanche motion is in the form of a Coulomb basal sliding friction, the application of D A N for modelling such flows is straightforward. The flow mass is divided into a discrete number of flow blocks, as in Figure 2.1. The mass blocks are indexed i = 1 to n — 1, the boundary slices j = 1 to n, starting at the upslope boundary of the flow and incrementing downslope. Each mass block has the same volume, and is separated from neighbouring blocks by infinitesimal boundary 'blocks.' These boundary blocks are essentially infinitesimal cross-sections used for the application of the momentum equation. The mass blocks, which are contained longitudinally by these boundary blocks and laterally by the specifiedflowwidth, are held at constant volume, assuming no entrainment. The flow path is approximated as a rectangular channel, with width defined as the surface flow width and input by the user. The flow depth, following the theory of open channel flow, is interpreted as the hydraulic depth. This allows the modelling of irregular channel cross sections, which in practice range from convex alluvial fans to narrow, confined channels. DAN employs a Lagrangian numerical scheme, such that the curvilinear coordinate system moves downslope with the flow. This approach is believed to be more appropriate for unsteady flow equations than a fixed-grid, Eulerian numerical scheme. A first-order finite difference discretization scheme of the equations of motion allows for easy numerical computation of position and velocity 18 Figure 2.1: Profile view of theflowdiscretization in DAN. Boundary slices between massflowblocks are oriented normal to the slope. updates for each boundary block. Figure 2.2 shows the geometry of a mass block and the adjacent boundary slices. Newton's second law acting on a boundary block is expressed as f (M v )=Y.Fit i i (2- ) 1 where Mj is the mass of the block i and Uj is its velocity. The sum of forces on the right hand side can be expressed as Z Fi = Mig sim!>i +Pi-Ti, (2.2) where g is gravitational acceleration, Pi is the internal pressure, and Tj is the basal resistance. The mass of a boundary slice is Mi = pHiBids (2.3) where p is the material density, Hi and 23, are boundary slice height and width, respectively, and ds is an infinitesimal length along the slope. The internal pressure term is written as Pi = -kM ^-{l i9 + —)cos^i, (2.4) 19 Figure 2.2: Geometry and forces acting on a boundary slice in DAN. where k is the longitudinal pressure coefficient. The term a = vf/R is the apparent centrifugal c acceleration that arises in the non-inertial Lagrangian reference frame for a curved path with vertical radius of curvature R. In Rankine earth pressure theory, the pressure coefficient k takes values in the range 0.2 for a static, 'extending' (active) mass to 5.0 for a compressing (passive) mass. These static states are analogously linked to extending and compressing flow states in DAN. In a hydrostatic fluid, k is always equal to 1. For granular flow in a state where the slope angle ib is less than the internal friction angle <f> of the material, internal deformations do not occur. For snow avalanche flow in runout zones, the condition ip < cf> likely holds true in most cases. Therefore, plug flow is expected and passive internal pressure is appropriate (McClung and Mears, 1995). For the following analysis, the internal pressure is held constant at the passive limit k = 2.5, corresponding to <j> = 25°, consistent with the analysis of p McClung and Mears (1995). DAN can modify k based on the state of longitudinal strain in the flow. For a description of the strain calculation, and the associated update in k, refer to Hungr (1995). The basal resistance term T can take a variety of forms in DAN. For dense snow avalanche flow, 20 Coulomb-like basal friction is assumed as the dominant resistance parameter. This takes the form Ti = \iMi(g com/a + (2.5) where p, is the Coulomb friction coefficient. Substitution of Equations 2.2-2.5 into a first-order discretization of Equation 2.1 leads to a solution for the velocity of a boundary slice at a new time step: , Vi where A t is the time step size, and = V i FjAt - Mj + pHiBids ' - (2 6) is the velocity of block i at the previous time step. Once the new velocities have been computed for each block at the new time step from Equation 2.6, the displacements Si are computed as Si = S'i + ^(v i + vl). (2.7) The continuity equation arises from the constant volume constraint of the mass blocks (for the assumption of no entrainment and incompressible flow), and allows the calculation of the average depth of the flow of a mass block: hj = 2Vi (Si+i-SiXBi^ + Bi)' (2.8) where V} is the volume of the mass block. The height of each boundary block is calculated as the average of adjacent mass blocks: H = t (2.9) The front and tail of the flow are tapered off such that the height of the outer boundary slices are each set at half the height of the adjacent mass block: #l = |,#n=^. (2.10) At this point the solution is complete for the current time step, assuming the internal pressure 21 coefficient is held constant. The solution proceeds downslope until all blocks have come to rest. 2.1 The New Modelling Procedure DAN was designed, as are most mass flow models, to simulate the flow from initiation in a starting zone to rest at some position downslope. The initial conditions for the calculations are thus v = 0 at some initial position s at t — 0. The simulation continues until the flow decelerates to rest and a again v = 0 at some position s = sfinal at time t = tji i. na As discussed in the previous chapter, a good deal of uncertainty exists surrounding the mechanisms of erosion and entrainment of existing snowcover by a passing dense snow avalanche. As a result, beginning a numerical simulation in the start zone of an avalanche path and neglecting entrainment may not produce a realistic picture of the flow dynamics. Sovilla et al. (2001) report a large variation between released mass and maximum mass, illustrating the importance of changing mass balance for certain events. The proposed method of analysis presented here attempts to avoid the unnecessary propagation of uncertainty attributable to entrainment. As stated above, maximum running avalanches are necessarily characterized by minimum resistance, and thus minimum entrainment. However, the importance of entrainment, even in these extreme cases, still cannot be quantified. Given increasing snowcover with elevation, it can be expected that most of the mass to be entrained by a dry slab avalanche will occur in the upper reaches of the path. By starting the numerical simulations in the middle of the avalanche path, rather than at the top, it is postulated that much of the uncertainty surrounding entrainment can be avoided. By selecting the position at which the flow might be expected to attain maximum speed, the calculations can proceed downslope during the decelerating phase of the flow. The initial maximum speed of the flow will be based on measured avalanche speeds. Entrainment is considered negligible while the flow decelerates compared to the acceleration phase of the flow in the upper reaches of the path. It should be stressed that this analysis is limited to extreme events only, for which resistive forces are inherently minimal. If the bulk basal resistance term is calibrated with extreme data, then the model implicitly includes inertial effects such as entrainment and ploughing (McClung, 1990). Therefore, the neglecting of entrainment and ploughing of erodible snowcover, in the runout zone 22 only, seems justified. Improved field observations and knowledge of entrainment in the future will determine the viability of this assumption. Entrainment has two deceleration effects that would need to be specified. First, the mass added to the flow needs to be accelerated from rest. To conserve momentum, the velocity of the flow must decrease with the increased mass. Second, there is a basal frictional force associated with the erosion of basal snow. Both of these effects will be difficult to measure and specify for inclusion in future dynamics models. 23 Chapter 3 Theoretical Considerations on a Smooth Parabolic Path A sensitivity analysis was performed to illustrate the role of the different physical parameters used in DAN and their dependence on speed and runout. A hypothetical parabolic path, defined by the equation y = .0003a; — 1.2a; + 1200, was considered. The vertical fall of 1200 m is close to that of 2 many large avalanche paths that threaten facilities with potentially destructive extreme events. The thick line in Figure 3.1 shows this hypothetical path profile. The representation of an avalanche path as a parabola, or second-order polynomial, allows the determination of unique terrain parameters for use in statistical prediction of extreme runout. This is due to the smoothly-changing nature of a parabolic approximation of three-dimensional terrain. Third-order or higher polynomial representation may involve changes in curvature and redundancy or ambiguity in the location of terrain parameters such as the /3-point (see Figure 1.1). Though refined dynamic modelling allows more realistic characterization of actual avalanche paths, a parabolic path is a good starting point from which to explore the behaviour of DAN. For comparison with previous investigations, some simulations were started from rest. This was also done to isolate certain aspects of the model before introducing additional complexity. Thereafter, simulations were initiated in the middle of the path as proposed. The coefficient of Coulomb sliding friction \i was determined following the depth-averaged, continuum mechanical derivation of McClung (1990, Equation (9)). Assuming that the turbulent drag term is negligible, \i = tan a. In this formulation, /x is independent of speed and only dependent of the path geometry, fi is interpreted as the average value of friction over the avalanche path. 24 3.1 Changing Path Width Given the dependence of all the contributing forces in Equation (2.2) on flow width, unit-width representation of an avalanche path should result in the same numerical predictions as any constantwidth channel. This confirmation is shown in Figure 3.1, which contains a comparison of avalanche speed for different constant channel widths for the parabolic path. All other initial dimensions (initial flow depth, length, friction, and starting position) were held constant. Channels ranging from 1 to 50 m in width were investigated. The flow was released from rest in order to guarantee that the only parameter being varied between simulations was the width. 0 500 1000 1500 Horizontal Reach (m) 2000 Figure 3.1: Frontal speed for variable path width. The thick line is the parabolic path profile. Identical speed profiles are shown for widths of 1, 2, 5, 10, 25, and 50 m. The initial flow depth was held constant, and each simulation was started from rest to ensure identical initial conditions. Using similar reasoning as above, changing the flow density while holding all other parameters constant will produce identical speed and runout. This is due to the mass dependence of each term in Equation 2.2. Though density is important in impact force calculations, it does not influence the dynamics of a Coulomb basal friction model. This is beneficial given the uncertainty in measured or estimated flowing avalanche densities. To investigate the effect of changing width over the course of a given path, the runout zone width of the same parabolic path was altered. The width of the track was held constant until the slope angle decreased to 16° at x = 1500 m. From this point, the width was linearly increased or decreased 25 until the slope angle reached 0° at x — 2000 m, after which the width was held constant. Figure 3.2 shows scaled runout versus the ratio of runout zone width to track width. When compared to runout for a constant width channel, runout distance increases as width decreases in the runout zone. Conversely, runout decreases for increasing channel width, given the same speed input. Granite Creek 13 is one such avalanche path that has decreasing width in the runout zone. This path intersects Highway 20 in Washington, USA. The ratio of runout zone width to track width for this path is between 1/3 and 1/4, depending on where in the track the width is measured. D A N predicts that any paths with a spreading ratio less than 1 will produce longer runout than if the ratio was higher. This behaviour can be explained by the form of the continuity equation used in the model. As the flow cross-sectional area decreases, speed must increase in order to maintain the same mass flow rate. _J I I I , I I I L_ a. •5 n: 0.18 2 4 6 Runout zone width / Track width Figure 3.2: Dependence of Runout Ratio on width change in the runout zone. A spreading ratio greater than 1 corresponds to width increase in the runout zone; less than 1 indicates constriction in the runout zone. 3.2 Internal Pressure Formulation The interest of the present investigation is on the lower portion of avalanche paths, for which the slope angles will generally be less than the expected internal friction angle of flowing snow. Given 26 this assumption, the body of the flow will be in a locked, or plug flow, state. Under these conditions, the flow mass will likely be in a state of passive internal pressure. DAN allows this assumption to be tested by either allowing the internal pressure coefficient of Equation (2.4) to vary between the active and passive states or by simply fixing k = k . Figure 3.3 p shows the effect of the form of k on runout for the parabolic path. The simulations were started near the middle of the path, such that the position of maximum speed occurred at nearly the same place for each run. The maximum speed (~62 m/s) corresponded to the prediction of the upper limit envelope (ULE) for the 1200 m vertical fall of the parabolic path. 0.18 1 1 . 1 —I 1 • 0.17 — • — .2 1 ! or 0.18 — — • 2 4 Initial Flow Depth (m) 6 8. Figure 3.3: Dependence of Runout Ratio on the form of the internal pressure coefficient for three different flow depths for the parabolic path of Figure 3.1. ("•_), k = k ; (•), variable k. The 0.03 range in Runout Ratio shown corresponds to a difference in horizontal reach of less than 50 m for the parabolic path. p The difference in calculated runout between the k = k and variable k simulations is small in p Figure 3.3. This indicates that when the internal pressure term was allowed to vary, based on the gradient offlowdepth along the flow profile (Equation 2.4), the results did not differ much from those using the passive value k = 2.5. This is confirmed in Figure 3.4, which shows the value of k for each v boundary slice along the flow profile for four different locations as the flow travels downslope for a 5 m initial depth simulation on the parabolic path. The value of A; deviates little from its initial value for the upper half of the flow. As the flow reaches lower-angle slopes and decelerates, k becomes 27 more variable along the flow profile, and the average value of k increases. 800 1200 1S00 Horizontal Reach (m) Figure 3.4: Variation of the internal pressure coefficient along the flow profile for four different positions along the parabolic path (thick line). (+), internal pressure coefficient k within theflowfor different positions of theflowprofile on the path. The effect of variable internal pressure on other flow features is more pronounced. Figure 3.5 shows the influence of variable k on flow and deposit lengths. For k — k , both flow length at p maximum speed and deposit length are greater than for variable A;. The decreased flow length for variable k is correlated with increased flow and deposit depths, as seen in Figure 3.6. The influence of a variable internal pressure formulation, with active and passive limits, is overall relatively minor. DAN simulations at constant passive pressure predict slightly longer, and thus more conservative, runout. However, the difference is small. Additionally, the appropriateness of the analogy between active and passive states of static earth pressure and the internal pressure of a granular flow is unclear. Given the simpler form of the equations of motion and numerical implementation for constant passive internal pressure, k = k = 2.5 was implemented in all further p simulations. 28 Initial Flow Depth (m) Figure 3.5: Dependence of flow length on initial depth for (+), k = k ; '(•), variable k. Deposit length vs. initial depth is also shown for (o), k = k ; (o), variable k. p p Initial Flow Depth (m) Figure 3.6: Dependence of meanflowdepth at maximum speed on initial depth for (+), k = k ; (•), variable k. Also shown is mean deposition depth vs. initial depth for (o), k = k ; (o), variable k. p p 29 3.3 Shape of Initial Flow Profile Starting a numerical simulation in the middle of an avalanche path presents the challenge of specifying the initial flow profile. A first approximation would be to specify constant initial depth, approximating uniform plug flow at or near maximum speed. Impact pressure measurements of dense avalanche flow, however, have indicated peaks in pressure, often in the front half of the flow (McClung and Schaerer, 1985). This suggests that there may be a concentration of mass toward the front of the flow, or simply that the front of the flow is travelling faster than the tail. The quality of such flow measurements is not high enough, however, to give an accurate representation of the gradient of depth, or velocity, along the flow profile. As an hypothetical exercise, however, Figure 3.7 shows four different initial flow profiles. Each of these profiles was used in a simulation of runout in DAN. Each simulation was started near the middle of the path, and the maximum speed and location of this peak were held constant. Table 3.1 shows the runout predicted by DAN for each case. J i I 1 1 16 —| 12 — H 0 100 200 Position (m) 300 Figure 3.7: Possible initialflowprofiles for a given volume, (f), constant initial depth; (k), leading peak depth; (M), trailing peak depth; (•), symmetric parabolic depth. Depths are increased by a factor of 20 for visualization. The variation in calculated runout is more than 100 m, indicating a high sensitivity to the shape of the initial flow profile. The excessive (and likely unrealistically high) flow depths of Figure 3.7 30 Initial Flow Profile Leading depth peak (A) Constant initial depth (•) Symmetric parabolic depth (•) Trailing depth peak (•) Horizontal Reach (m) 1974 2003 2098 2100 Table 3.1: Calculated runout distance for the initial flow profiles in Figure 3.7 exaggerate this variance in calculated runout. However, as will be discussed below, the flow thins considerably by the time maximum speed is reached after the start of each simulation. Given impact pressure records, the profile with maximum depth at the rear of the flow is not likely realistic. Given the lack of high quality data from which to infer flow profiles, it is hard to make an argument for or against any of the other flow profiles, however. Given this lack of data, and the expectation of plug flow on the lower slopes of interest in this investigation, the initial depth for further simulations was held constant. This simplification seems to be a reasonable first approach, and less arbitrary than attempting to specify an initial depth gradient without better data. 3.4 Dependence of Flow Depth on Runout For a given starting position, the initial flow depth was altered for the parabolic path. The path width, initial flow length, and initial speed were held constant. Constant-depth longitudinal flow profiles (• in Figure 3.7) were used for the simulations. Figure 3.8 indicates the strong dependence of flow depth on both maximum speed and runout distance in the model. It also points out an interesting two-dimensional characteristic of DAN. Note the short period of acceleration after the initial speed was given to the flow at the first time step. This is an unrealistic characteristic of the simulation related to the thinning of the flow under its own weight, aided by the internal pressure formulation in DAN (see Figure 2.2 and Equation 2.4). In other words, the assumed initial flow profile and speed do not characterize a 'stable' DAN hydrograph. A hydrograph characterizes flow speed and depth along the entire flow profile. A stable hydrograph is a characteristic of a fully-developed flow that is travelling at approximately constant speed. However, accurate hydrograph data for snow avalanches do not exist. Even if such data were available, a measured hydrograph for a snow avalanche would not necessarily correspond to a stable 31 hydrograph in DAN. 0 300 600 900 1200 1500 1800 Horizontal Reach (m) Figure 3.8: Variable initial depth for a given starting position and initial frontal speed. (•), 10 m; 5 m; (+), 3 m; (A.), 2 m; (M), 1 m; (if), 0.5 m; (—), path profile. The acceleration phase of the speed profiles in Figure 3.8 was observed for all parameter combinations. Upon beginning the numerical calculations with a non-zero initial speed, the flow profile has a tendency to thin and elongate. This behaviour is coupled with the observed acceleration. For a smooth slope profile, as in the parabolic path used above, the depth profile elongates and approaches an approximately parabolic shape. Upon reaching maximum speed, the flow has settled into what might be termed a 'stable' flow profile. After this point the change of flow speed, depth and length are primarily correlated to changes in the path. The physically realistic portion of the numerical simulation thus begins at the point where the maximum speed is reached. The short period of initial acceleration, present in all simulations, is a result of a deviation in the user-input flow profile from what would be a stable hydrograph in DAN. For actual avalanche paths, minor undulations in slope and curvature retard the flow thinning caused by the internal pressure, and initial acceleration is minimal. This will be evident in Chapter 4. Figure 3.9 isolates the "thinning" effect that gives rise to the acceleration phase seen in Figure 3.8. In this case, a 100 m long, 5 m thick profile divided into 15 mass blocks was placed on a 0° flat slope. The basal friction coefficient was set at 0.3. With no other external forces or initial speed input, the flow thins due only to the internal pressure term. Over a period of about 8 seconds in the 32 simulation, the front and tail of the flow reach unrealistically high speeds of about 10 m/s before decelerating to rest. The length of the mass increases by about 80%. This alarming behaviour raises questions about the widely adopted formulation of internal pressure used in this model and many others. As evident in Equation 2.4, the internal pressure term depends on a gradient in flow depth, dH/ds, at a given position. Initially in Figure 3.9, the only positions along the mass profile that have non-zero values of dH/ds are at the ends. This is a result of the definition of the end blocks of any flow mass in DAN; the height for the first and last blocks are set at half of the height of the respective adjacent internal block. This tapering of the flow depth at the front and tail of the mass invokes a non-zero internal pressure via the depth gradient, which acts to displace the exterior boundary blocks. This, in turn, introduces a non-zero depth gradient between the next two interior blocks. The flow thus continues to deform under the action of this internal pressure, and stops when the basal friction coefficient exceeds the internal pressure. For fi — 0, the flow mass would continually thin without bound, approximating the classic dam-break problem in fluid mechanics. This thinning feature is present in DAN for any depth of flow, but is more pronounced for greater depths. This explains the higher maximum speed and runout distance for the highest initial depth in Figure 3.8. This also calls into question the appropriateness of the formulation of the internal pressure term, if only for an initially stationary mass. Though quasi-static alpine snow does thin under the influence of gravity through a process called creep (McClung and Schaerer, 1993), the timescale over which this happens is much longer than a timescale on the order of 1-10 s, as in Figure 3.9. It is not clear whether it is appropriate to draw an analogy between a static granular pressure and a flowing granular pressure, though this approach is quite standard in the literature on avalanche dynamics. Previous dynamics models have all started calculations from rest in the start zone, and have perhaps not noticed the implication of the internal pressure formulation brought to light by starting DAN in the middle of the path. Simple tests of the physical nature of a numerical model, such as shown in Figure 3.9, do not seem to have been used in order to argue for or against any similar pressure formulations. The only obvious alternative would be to abandon any attempt at describing and implementing an internal pressure term in the equations of motion, owing to the lack of a clear and convincing 33 10 - , so ISO Position (m) 200 Plow ttiinning under gravity • — • • Initial profile # — # — • Final prof la Frontal speed Tail speed Figure 3.9: Flow thinning under gravity and resultant frontal and tail speeds. physical justification for including such a term in a granular flow model. This was done by Louge and Keast (2001) for granular flows, without any loss of physical accuracy. If the volume fraction of solids is less than 50%, it is hard to justify a flowing active/passive pressure analogy with static soils that are in excess of 90% solids. Setting this argument aside, most current avalanche dynamics models include some form of internal pressure. In the absence of a clear alternative, and largely for comparative purposes, the pressure formulation as outlined in Equation 2.4 was retained, with the internal pressure coefficient k held constant at the passive value k = k = 2.5. p Figure 3.10 illustrates an apparent upper-bound on flow depths once the maximum speed has been reached. Both the average and maximum depth taper off for increasing initial depth. The average depth has an apparent upper-bound near 4 m, and the maximum depth is slightly more than 6 m. This does not necessarily correspond to physically realistic bounds on snow avalanche flow. However, it suggests that such limits on flow depth may exist and have a fundamental physical explanation in the balance of forces. The Upper Limit Envelope (ULE) determines the desired maximum speed of the flow in the middle of the path. However, due to the flow thinning and resultant non-physical acceleration, as described above in relation to the internal pressure formulation, the initial speed of the flow must be adjusted so that the desired maximum speed is reached. This requires a trial-and-error iterative procedure for initial speed. Figure 3.11 shows the different runout distances predicted from different 34 • < ++ • t + • + ! Initial depth (m) Figure 3.10: Average (+) and maximum (+) flow depths at maximum speed vs. initial depth. initial starting positions in the middle of the parabolic path. The maximum speed for each simulation is 62 m/s, as predicted by the ULE for a 1200 m vertical fall path. As might be expected, for a given maximum speed, increasing the horizontal starting position leads to increased runout distance. However, the starting positions vary by 400 m along the path, but the predicted runout distance varies by little more than 50 m. This indicates only a weak sensitivity of the runout distance to the starting position of the numerical simulation. The speed profile of a one-dimensional approximation of a Coulomb frictional material will achieve a peak at the position where sin ip — p, cos ip = 0, in other words when the driving force exactly equals the flow resistance. The point at which sin ip — p cos i/> = 0, or where the slope angle ip = tan~~ p, l thus divides what might be termed the 'acceleration' and 'deceleration' regions of the path. Upslope of the position where ip — tan p, _1 the gravitational driving force is greater than the basal frictional force, and vice versa in the deceleration region downslope. It seems reasonable to extend this one-dimensional reasoning to DAN in the selection of the starting position. If the maximum speed is expected to occur at or near the point where ip = tan p, _1 then the starting position in DAN can be chosen accordingly. However, in the two-dimensional description of the flow in DAN, each flow block can have a different speed. For instance, the front, tail, and centre of mass can all reach different maximum 35 Figure 3.11: Speed profiles for variable starting positions and given maximum speed. Frontal starting positions: (i), x=600 m; (M), x=700 m; (•), x=800 m; (k), x=900 m; (+), x=1000 m. (—), path profile. speeds at different locations along the path. Which flow block achieves the highest overall speed will vary depending on the path in question- Thus there is some ambiguity in which part of the flow should be taken as 'representative' for monitoring speed. Impact pressure measurements of avalanche flow (McClung and Schaerer, 1985) indicate that the peak pressure, and likely corresponding peak velocity, usually occurs in the first half of the flow. It seems reasonable, then, to manipulate the starting position of the simulation so that either the frontal speed or centre of mass speed coincides with the location where sin tp — p cos tp = 0. Schaerer and Salway (1980) found that an initial peak pressure in impact sensor monitoring could, on average, be correlated with frontal speed. Though uncertainty remains, for further simulations, the starting position was modified so that maximum frontal speed occurred as nearly as possible to the location where the slope angle ip = ton /i. _1 Once the above procedure was calibrated, the influence of flow depth on runout distance was investigated. The maximum and average depths at maximum speed were recorded and compared with the runout ratio. The qualitative nature of both the maximum and average depth groupings, in Figure 3.12, is similar. Runout distance seems to peak at a characteristic flow depth before tailing off. However, this is a small effect. For the parabolic path considered, the range in horizontal reach corresponding to the Runout Ratios shown in Figure 3.12 is less than 40 m. 36 0.18 — 0.17 I I i — | • - ; - + + CC "3 0.16 i j - + - r - • i • — Run o I i t 0.16 — ! 0.14 — 0 I | 2 ! • I | 4 1 Depth at maximum speed (m) | 6 1 Figure 3.12: Runout Ratio vs. maximum (+) and average (+) flow depth. Each symbol represents a unique, independent simulation on the parabolic path. Both maximum and averageflowdepth were recorded when the flow reached maximum speed. The Runout Ratio increases for increasing flow depth (both maximum and averageflowdepth), up to a certain threshold. Above this threshold (about 3.5 m for averageflowdepth and about 5 m for maximumflowdepth), the runout for subsequent simulations decreased for increasingflowdepth. Any change in depth as the flow decelerates could be an important factor influencing runout distance. Figure 3.13 shows the change in depth from the middle of the path where peak speed occurs to the depth of the deposit in the runout zone. As seen, the deposit depth is only weakly dependent on flow depth. For initial flow depths of less than 4 m, the resulting deposit depth increases. This is the case for both the maximum and average depth, indicating that the flow length decreases, or the flow piles up, as it slows down with decreasing slope. For initial flow depths of greater than 4 m, the deposit depth decreases slightly, indicating moderate flow thinning. It is important to remember that these results are for a unit-width channel. Therefore, the effects discussed above are independent of width variation along the path. Changing width along the path introduces an additional influence on flow dynamics. Variable width along a path, as seen in Figure 3.2, has a greater influence on flow speed and horizontal reach than the flow-depth effects discussed in this section. Passive or plug flow behaviour for certain flow depths is seen in Figure 3.14. The flow depth is scaled with the deposit depth and compared with runout. A depth ratio of 1 most closely ap- 37 E. £ • • • • + + . , 1 2 4 Flow Depth (m) 6 Figure 3.13: Maximum deposit depth vs. maximumflowdepth (+) and average deposit depth versus averageflowdepth (+). Each symbol represents a unique, independent simulation on the parabolic path. Both maximum and averageflowdepth were recorded when the flow reached maximum speed. The deposit depth shows little sensitivity to changingflowdepth. proximates passive, or plug, flow. For average and maximum depth ratios near 1, there is a peak in runout. Flows with increasing depth in the runout zone (depth ratio less than 1) did not rim out as far. Physically, this often corresponded to a flow profile with a thin front. The flow front first comes to a rest, followed by the tail of the flow piling up from behind. Depth ratios greater than 1, with thinning of the flow in the runout zone, also did not run out as far. 38 • t • + + ++ 0.5 1 1.5 Flow Depth at Max Speed / Deposit Depth Figure 3.14: Runout Ratio vs. scaled maximumflowdepth (+) and scaled averageflowdepth (+). Each symbol represents a unique, independent simulation on the parabolic path. A scaled depth of 1 represents no change in depth (for either maximum or average depth) as the flow decelerates to rest in the runout zone. This condition roughly corresponds to maximum Runout Ratio. 3.5 Dependence of Flow Length on Runout The flow length, scaled with the deposit length, does not exhibit the same dependence as the scaled flow depth when compared with the runout ratio. What might be termed a passive or plug flow in this sense, with a length ratio (flow length / deposit length) of 1, corresponded with shorter runout (Figure 3.15). The longest runout was for flow lengths that were about 25% longer than the deposit length. The explanation for this lies in the depth change of many of the flow blocks. Though the longest runout was for flows in which the maximum depth did not change, the average flow depth generally approached the maximum depth as the flow decelerated in the runout zone. This led to a slight decrease in flow length. When scaled with the curvilinear path length, the deposit length was roughly constant for changing flow lengths, as seen in Figure 3.16. However, the runout ratio increases for increasing scaled flow length, as evident in Figure 3.17. This was primarily a function of the increase in mass pushing the flow from behind as the front of the flow decelerated upon reached lower slope angles. 39 —i—i—i • • _i—i i • • • —i—i—i—i——i—i—i—i— ' 1.6 2 2.5 3 Flow length at maximum speed / deposit length 3.5 Figure 3.15: Runout Ratio vs. scaledflowlength. Each symbol represents a unique, independent simulation on the parabolic path. A scaled length of 1 represents no change in length as the flow decelerates to rest. Maximum runout was observed for a scaled length between 1 and 1.5, indicating a slight contraction in length on the longest running flows. 0.05 0.1 0.15 0.2 Flow length I Path length Figure 3.16: Scaled deposit length vs. scaledflowlength at maximum speed. Each symbol represents a unique, independent simulation on the parabolic path. The deposit length shows little sensitivity to theflowlength. 40 S 0.175 0.1 0.15 Flow length / Path length Figure 3.17: Runout Ratio vs. scaled flow length. Each symbol represents a unique, independent simulation on the parabolic path. Increasing flow length corresponded with increased runout. 41 3.6 Variable Number of Mass Flow Blocks DAN allows a given flow mass to be discretized into between 1 and 59 flow blocks for numerical simulation. Representing the flow with a single block corresponds to a rigid, sliding slab. Two blocks would represent that initial slab broken in two pieces, and so on. The greater the number of blocks that a given mass is divided into, the more fluid-like the flow behaves. Figure 3.18 shows the difference in calculated speed and runout distance for varying numbers of blocks, given the same initial mass and flow profile. The maximum speed was determined by the ULE, and the starting position of each simulation was modified to match the position of peak speed at the point where i/> = tan' 1 fi, as discussed above. 0 400 800 1200 1S00 Horizontal Reach (m) Figure 3.18: Speed profiles and runout distances for a givenflowmass divided into different numbers of massflowblocks. (•), 5 blocks; 15 blocks; (k), 25 blocks; (D), 40 blocks; (o), 50 blocks. As the number of blocks that a given mass is divided into increases, the flow thins more rapidly after the start of the simulation. In Figure 3.18, the simulation for the 5-block flow thins gradually, producing a smooth, rounded speed profile. For the 50-blockflowdiscretization, the flow thins more rapidly and produces a sharper initial acceleration. This is evident in the higher initial speed (about 45 m/s) and further downslope starting position necessary for this simulation to match the ULE maximum speed where tp = tan" p. 1 In Figure 3.18, the difference between the initial speed and the given maximum speed seems to be reaching a lower limit as the number of blocks discretizing the flow increases. The difference is small 42 between the speed profiles for the 40 and 50 block calculations. The apparent minimum mismatch between the initial speed and maximum speed in DAN is an isolation of the influence of the internal pressure formulation on the calculated flow behaviour. If the internal earth pressure analogy were neglected in DAN (i.e. Pi = 0 for all i in Equation 2.4), there would be no internal mechanism to promote flow thinning (see Figure 2.2) and the accompanying frontal acceleration on a smoothly changing path profile, such as the parabolic path considered in this chapter. Figure 3.19 shows the profile of the deposited mass in the runout zone for four of the simulations in Figure 3.18. Figure 3.19(a) and 3.19(b) show thinning at the front of the deposit. Figure 3.19(c) and 3.19(d) more closely resemble a plug flow mass, where the shape of the deposit is much the same as the flow profile. However, the deposit length and maximum depth of (c) and (d) are quite different. 1500 1600 1700 1600 1900 2000 Horizontal Reach (m) Figure 3.19: Deposit profile for the same total mass divided into different numbers of mass flow blocks, (a), 40 blocks; (b), 25 blocks; (c), 15 blocks; (d), 5 blocks. 43 The data in Figures 3.18 and 3.19 do not contain compelling enough information to suggest that a specific order of discretization is appropriate. This characterization of the flow seems largely arbitrary given the lack of observational data or a physical argument to use 15 blocks rather than 5 or 25, or any value in between. If complex avalanche terrain with small and/or rapidly changing radius of curvature is to be modelled, then a minimum number of flow blocks may be necessary to maintain numerical stability. This is not, however, a physical argument, even though the choice of the number of flow blocks influences the physical character of the flow. For numerical stability over quickly changing terrain, Hungr recommends using 15 blocks for simulations. As this value seems appropriate, given the expected passive or plug flow behaviour in avalanche runout zones, it was adopted for further simulations. For flow lengths of the order 10 m, 2 this leads to a mass block length on the order of 10 m. 3.7 Variable Initial Speed Profile Speed profiles of flowing avalanches are not well known. From impact force measurements (Schaerer and Salway, 1980; McClung and Schaerer, 1985), qualitative speed characteristics can be inferred. The first half of the flow mass is likely traveling near the maximum speed when the flow is in the middle of the track. The speed of the rear of the flow is probably lower, as indicated from pressure records that reveal lower impact pressure in the tail of the flow. These features axe difficult to quantify, however. It is assumed, therefore, that any gradient in speed along the flow length will have a peak near the front when the avalanche is in the track. This is valid assuming a smoothly changing path profile in the absence of abrupt width changes. In reality, this may not always be the case. As the front of the flow reaches regions of locally higher friction or lower slope angles, it may be possible for the middle or rear of the flow to be traveling faster than the front. Additionally, changes in path confinement can lead to higher speed in the rear of the flow. DAN has the capability of demonstrating these phenomena. However, for the purpose of assigning initial speeds to each flow block, a simple approach is desired as a first approximation. For each simulation, the highest speed is given to the flow front. The speed given to successive blocks is then decreased to a specified minimum at the tail of the flow. Two methods of decreasing 44 the initial speed from the flow front to the rear were investigated. Both a linear decrease and an exponential decrease in speed were attempted for the initial time step. The speeds were adjusted such that the rear speed was a given percentage of the frontal speed. The initial speed of the flow front was held constant, as were the initial starting position and depth of flow. Figure 3.20 shows the effect of different speed profiles on the runout ratio for the parabolic path. Runout was longest when each flow block received the same initial speed. As the speed of the tail decreased relative to the front, the runout distance decreased. If the rear of the flow is initially traveling at 50% of the speed of the flow front, the runout distance is about 75 m shorter than if the rear was initially at the same speed as the front. I O i £ I I I I I J U _ 0.16 - 0.14 - 0.6 0.7 0.8 Initial tail speed / Initial frontal speed 0.9 Figure 3.20: Runout for variable initial speed profiles. (•), linearly decreasing initial speed from front to tail offlow;(+), exponentially decreasing initial speed. 3.8 Discussion It is clear that many parameters influence the results in DAN. The motivation for starting the numerical simulations in the middle of the avalanche path is reasonable. The elimination of a good deal of uncertainty surrounding the dynamics of basal erosion and entrainment is desirable. Most previous dynamic modelling efforts have neglected these effects entirely, without discussing the propagation of uncertainty in calculations of runout distance and speed. Others have attempted 45 to model the processes of erosion and entrainment, though the formulations proposed still involve numerous assumptions about processes that are still not well understood. If erosion and entrainment are indeed as significant to the mass balance as recent observations suggest, then it will perhaps no longer be appropriate to neglect these effects in modelling scenarios. The approach employed here rests largely on the formulation of the maximum speed for recorded avalanches. A secondary assumption is that entrainment is a significant dynamic process in dense flowing snow. However, it is assumed to be most prevalent in the upper reaches of an avalanche path when the flow is still accelerating. Even if entrainment is minimum for extreme avalanches, this does not mean it can be neglected entirely. Given these assumptions, it seems promising to use a numerical model for only the lower half of an avalanche path. However, the approach employed here introduces unique complexities that introduce additional uncertainty. The net effect on uncertainty is difficult to quantify. Given the lack of high quality, high resolution data on the three-dimensional (or even twodimensional) nature of avalanche flow, it is difficult to calibrate the procedure for using DAN given this new technique. Spatial characteristics such as length, depth, and speed gradient have a varied but quantifiable effect on the results. Basic parameters such as the flow depth and length, or pseudoparameters such as the initial speed gradient or degree of flow discretization, each individually contribute to an uncertainty in runout on the order of 10-100 m. However, these uncertainties are not independent, and the net effect on the uncertainty in calculations is difficult to quantify. In the absence of better observational data, many of these flow characteristics must be assumed or qualitatively deduced. It is difficult, if not impossible, to compare the magnitude of uncertainty in neglecting entrainment when starting a dynamics model from rest with the proposed method here. Improved observations of avalanche behaviour will likely lead to reduced uncertainty with both methods; it will remain to be seen which approach will will provide better results. The challenge of model comparison and validation will be explored in the next chapter. 46 Chapter 4 Model and Procedure Validation Validation of any numerical model requires observational data of the physical process of interest. Often, analogous physical processes are used as a basis for model verification. This is often the case for dynamics models of flowing snow due to the relative scarcity of spatial and temporal snow avalanche data. Granular flow theory and experimental chute flow data are often used to calibrate and test dynamics models for flowing avalanches. A model designed for simulating extreme events is presented with an even greater challenge. Long return period avalanches are, by nature, rare in human experience. This fact is reflected in the scarcity of recorded observations for such events. In cases of recorded extreme avalanches, standards and quality of collected data vary. The most common reported feature of an extreme avalanche is the frontal stopping position. This may be accompanied by a size classification. More scarce are data on the speed of the flow, or geometric specifications of the avalanche deposit or slide path. To fully verify avalanche dynamics models of more than one dimension, higher data quality will be necessary than are presently available. Spatial measurements of the avalanche path must be coupled with spatial and temporal bounds of the flow. Such data is rare for dense snow avalanches. Visual observations of the dense core of a dry slab avalanche are usually obstructed by an enveloping turbulent powder cloud. Impact-pressure records or microwave measurements of avalanche flow provide a picture of the internal flowing core that is unavailable to the human eye. Such methods can give indications of spatial flow parameters such as depth and length, albeit with some uncertainty. Most avalanches recorded with such means, however, are not extreme events. Given the sensitivity of computed runout distance to spatial characteristics such as flow length, depth and width, as seen in Chapter 3, a full verification of the predictive potential of DAN is in advance of observational data. Nonetheless, the model can be compared against historical observations of some events. Some modelling parameters will remain poorly defined for the time being, 47 but reasonable bounds can be defined with which to compute upper- and lower-bounds on runout distance and speed. 4.1 Validation of D A N With Simulations Starting From Rest Comparison of DAN against existing models is one way to demonstrate the applicability of the model for snow avalanches. Examples in the literature are present that provide measurements of avalanche speed and runout distance, coupled with various model predictions or back-calculations. Before attempting to verify the method of beginning the calculations in the middle of the path, DAN was compared against other models and measurements of speed, beginning from rest in the starting zone. In this manner, it can be determined whether or not the model is appropriate for simulating flowing avalanches under conditions in which entrainment can be neglected. For the Arabba avalanche path in Italy, the 21 December, 1997, avalanche was chosen for comparison. The measured speeds were reported by Sovilla et al. (2001), and digitized for this analysis from Ancey and Meunier (2004), together with the interpolated path profile. The upper limit of the fracture zone was taken as the highest point of the recorded path profile, and the lower limit was correlated with the beginning of the interpolated velocity profile (Ancey and Meunier, 2004, Figure 5). Initially, the friction coefficient was chosen based on the path geometry, i.e. / i = tan a. The friction was then adjusted to reproduce the observed runout. Figure 4.1 shows the results. For /i = tan a — 0.66, DAN over-predicts the speed along the path, and at the end of the path the front of the flow is still traveling in excess of 10 m/s. For ji = 0.674, the runout distance is matched. Similar to the analysis in Ancey and Meunier (2004), three digits of precision were often needed to reproduce a given runout position in DAN. For this reason, a dynamics model is not the best means for predicting runout. Rather, specifying runout using an empirical method and then applying a dynamics model for scaling avalanche speed is more appropriate for applications (McClung, 1990). The maximum speed, and position of this peak, are closely matched for both values of p in Figure 4.1. In the runout zone, however, the calculated speed is greater than that measured. This deviation is due to the tail of the flow catching up to the front as the slope angle decreases in the 48 Figure 4.1: Comparison of measured (•) and computed frontal speeds for the Arabba avalanche of 21 December, 1997. The thick line is the path profile. (—), speed profile for \i = tan a = 0.66; (-), observed runout is matched for fj, = 0.674. In the runout zone, the predicted speed in excess of measured speed may be a result of either the constant model friction being less than the true friction or changes in flow width that were not accounted for here. DAN simulation. This gives an extra push to the front of the flow, and relatively high speed is carried far into the runout zone. This feature could be corrected with increasing friction as the flow decelerates, as suggested by McClung (1990) The upper limit envelope (ULE) for maximum avalanche speed predicts a maximum speed of 42 m/s for an extreme avalanche given the same vertical fall as the Arabba path. Starting from rest in DAN with fi — tan a, the maximum calculated speed is about 26 m/s. This is a considerable difference. Without further information about the Arabba path, such as changing width or surface roughness, there is no obvious explanation for this discrepancy. A sensitivity analysis of runout distance to the friction coefficient fi was performed for the Arabba path. This allows comparison of DAN predictions with the model of Ancey and Meunier (2004, Figure 12). Figure 4.2 shows these results. For a given value of p, DAN predicts greater runout distance than the center of mass model of Ancey and Meunier (2004). This is no surprise, as the runout distance in DAN is reported as the frontal stopping position, necessarily ahead of the center of mass. Qualitatively, the shapes of the two profiles in Figure 4.2 are similar, however. This reflects similar Coulomb frictional approaches for both cases. 49 0.6 0.65 0.7 0.75 0.8 0.85 0.9 Figure 4.2: Variation in runout distance with changing friction coefficient p for the Arabba avalanche path. (—), Ancey and Meunier (2004); (t)> DAN. A similar analysis, starting from rest, was carried out for the Madergrond avalanche path in Switzerland. Speed data for three events were reported in Gubler et al. (1986). The 17 January, 1985 avalanche was chosen for comparison, and the speed data and path profile again digitized from figures in Ancey and Meunier (2004). Figure 4.3 shows the results. Similar to Figure 4.1, DAN predicts runout greater than that measured for p = tan a. This again points to the danger in using a dynamics model to predict runout. The computed maximum speed for the lower friction run in Figure 4.3 is close to the measured maximum speed, but the shape of the two speed profiles is quite different. The calculated maximum speed (46 m/s) is near that predicted by the ULE (49 m/s). For p, = 0.576, the observed runout distance is matched by DAN. The speed profiles in the start zone and runout zone agree well, but DAN fails to predict the sharp peak in speed in the middle of the path. There are two possibilities for this discrepancy. First, the model uses constant friction. In reality, the friction is likely speed-dependent (McClung, 1990). An initial low value of friction would lead to higher calculated speeds in the middle of the path. Increasing friction would then lead to the sharp observed deceleration of avalanches in the runout zone. Second, the calculations in Figure 4.3 assume a unit-width path. In reality, most avalanche paths have variable width and confinement as 50 2600 0 200 400 600 800 1000 1200 1400 1600 Horizontal Reach (m) Figure 4.3: Comparison of measured (•) and DAN frontal speeds for the Madergrond avalanche of 17 January, 1985. The thick line is the path profile. (—), speed profile for fx = tan a = 0.466; (- -), observed runout is matched for fi — 0.576. a function of position. The sensitivity of DAN to changing width was demonstrated in Chapter 3. The difference in speed estimates for a unit-width versus a true-width calculation in DAN could, in some cases, be significant. For the Madergrond simulations in Figure 4.3, note the difference in the position of peak speed between the two DAN simulations. For the higher friction simulation, the position of maximum speed is about x = 450m, whereas for the \x = t a n a = 0.466 simulation, the peak frontal speed occurs beyond x = 1100m. This feature highlights the difficulty later encountered when attempting to manipulate the position of peak speed by starting the calculations at different points along a given path. For a smooth, gradually changing path profile, there are no sub-peaks of locally high frontal speed. However, for a path profile that has frequent changes in curvature, the position of peak speed is largely dependent on the choice of starting position and the friction coefficient. The Madergrond avalanche path was also modelled by McClung (1990). A comparison of McClung's model results with DAN is shown in Figure 4.4. As predicted by McClung, the use of constant mean friction in DAN for this case required a separate, lower value of \i (n = 0.2) in the start zone in order to allow the slide to accelerate. Similar to McClung's constant-friction run, DAN predicts speed in the middle of the path which underestimates the measured speed. The sharp deceleration 51 of the speed profiles in the runout zone are well matched for each case, however. 2600 —I 0 1 1 400 1 I 800 I Horizontal Reach (m) I 1200 1 L 1600 Figure 4.4: Comparison of DAN versus the model of McClung (1990) for the Madergrond path. (—), speed measured by Gubler et al. (1986); (+), variable-friction model of McClung (1990) without a surface drag term; (O), constant mean friction model of McClung (1990) without a surface drag term; (k), DAN frontal speed. Two examples of avalanche paths for which DAN under-predicts runout for p = tan a are shown in Figures 4.5 (Khibins path) and 4.6 (Aulta path). The Khibins speed data were reported in Kotlyakov et al. (1977), and the Aulta speed data originated in Gubler et al. (1986). Path profiles and speed measurements were again digitized from Ancey and Meunier (2004). In Figure 4.5, the runout distance is reproduced for p, = 0.449, which is less than 3% smaller than tan a. The speeds in the start zone and runout zone are well matched, but DAN does not reproduce the two measured peaks of speed in the lower half of the path. These peaks could be related to slide path features such as changing width, but the absence of such data makes any remarks purely speculative. The maximum speed could be reproduced with a speed-dependent friction formulation. However, the two distinct peaks of measured speed would not likely be captured by a unit-width model using a gradually changing friction formulation. The maximum speed calculated by DAN is almost 20% lower than that predicted by the ULE in this case. This could in part be due to the scaling of maximum speed with the total vertical fall of the path (McClung and Schaerer, 1993) rather than with the total curvilinear path length, as 52 Figure 4.5: Comparison of measured (•) and DAN frontal speeds for the Khibins (a) avalanche. The thick line is the path profile. (—), speed profile for fi = tan a = 0.462; f- .-), observed runout is matched for \i = 0.449. proposed by McClung (1990). Unfortunately, reliable path length data was not available in DAN, given the model's spline-fit method of approximating and smoothing path profiles. In Figure 4.6, the runout distance is matched for a value of /x about 7% less than tan a. The calculated speeds are well in excess of measured speeds, and about 30% greater than that predicted by the ULE. Qualitatively, however, the shapes of the speed profiles are similar. The sharp acceleration and deceleration phases of the measured speed are well represented by the model. The higher calculated speeds in Figure 4.6 could be the result of the absence of a surface drag resistance term in DAN. For the high speeds produced by the large vertical fall of this path, a surface drag term would limit the maximum speed, perhaps into better agreement with observed speed data. Figure 4.7 shows a comparison of DAN with the variable-friction model of McClung (1990) (including surface drag) for the Aulta path. The addition of a surface drag term in McClung's model does not bring the calculated speed in the middle of the path into any better agreement than DAN without surface drag. Both models agree well with the measured speed in the deceleration phase, however. 53 Horizontal Reach (m) Figure 4.6: Comparison of measured (•) and computed frontal speeds for the Aulta avalanche of 8 February, 1984- The thick line is the path profile. (-), speed profile for p = tan a = 0.414; (- -), observed runout is matched for p. — 0.384. 4.1.1 Full-Simulation Discussion For the historical avalanches considered, DAN was able to reproduce the recorded runout position for friction values close to p — tan a. However, this does not imply that the model should be used to predict runout. Rather, given a specified runout position, the value of p that will reproduce this runout may be expected to be close to tan a. Calculated and measured speeds agreed well for the start zone and runout zone in most cases. Differences between predicted and measured speeds appeared for a number of reasons. For the Aulta path, the calculated speed in the middle of the path was higher than the measured speed. This could be amended, in part, by including a drag resistance term at the top of the flow. This was done by McClung (1990), which brought his calculated speed into better agreement with measured speed. In other cases, such as the Madergrond and Khibins paths, DAN predicts speeds which underestimate the measured speed. McClung (1990) addressed this problem by introducing a variable friction formulation, with low initial friction that increases as the flow travels downslope. Such a formulation in DAN would likely lead to better agreement with the maximum measured speed in the middle of these paths. Caution should be exercised in applying such an analysis, however, while modelling the avalanche path with unit-width. The sensitivity of DAN to changing path width, as 54 Figure 4.7: Comparison of DAN with the model of McClung (1990) for the Aulta path. (—), speed measured by Salm and Gubler (1985); (- - -), variable friction model of McClung (1990), including surface drag; (k), DAN with constant friction fi = 0.446 and no surface drag term. seen in Chapter 3, should warrant an accurate characterization of the path width before introducing additional complexity in the friction formulation. It is not likely that a simple, variable friction model would reproduce the sharp, local peaks of speed that were measured in the Madergrond and Khibins paths. Such peaks could arise, however, from changes in path confinement, which might be easily captured by a variable-width model with constant friction. 4.2 Validation of the Proposed Method In order to avoid making unnecessary assumptions about entrainment, the model was started near the middle of the avalanche path with some initial speed v > 0. The initial speed was adjusted such that once a stable flow profile was reached, the maximum speed was equal to the ULE prediction for the given vertical fall. Alternatively, the maximum speed could be scaled with the curvilinear path length given a confident measure of this length. For the selection of the friction coefficient, two possibilities emerged. The first was to select the friction coefficient according to fi — tan a. The second was to run a simulation from rest in the start zone, and adjust the /z such that the specified avalanche stopping position is reproduced. This value 55 of \x could then be used in a simulation starting near the middle of the path. In either case, given a constant value of /x, the starting position and initial speed for the simulation was adjusted such that the desired runout position was reached and the maximum flow speed corresponded to either the prediction of the ULE or a local, measured maximum speed. It was found that there was some redundancy in the method of selecting the starting position of the simulation. In many cases, several different starting positions in the middle of the path reproduced approximately the same runout. An example of this is seen in Figure 4.8. The friction coefficient p was held constant, and the initial speed was adjusted such that the maximum stable flow speed was as near as possible to that predicted by the ULE. The two simulations that were started highest on the slope produced speed profiles that were very similar in the runout zone. The simulation starting the farthest downslope, however, has visibly higher speed for much of the runout zone. The horizontal reach for this case is slightly farther than the other two, however. Horizontal Reach (m) Figure 4.8: Speed profiles for different starting positions for the Madergrond avalanche path. Arrows indicate the position on the slope that each calculation began from. Each calculation had the same maximum speed and constant friction. To constrain the selection of the starting position, the position of the peak calculated speed was matched as closely as possible to the slope position that corresponds to ip = tan" 1 /x, as in Chapter 3 for the hypothetical parabolic path. Figure 4.8 shows the difficulty in manipulating the peak speed position for paths such as Madergrond, especially if the average track slope angle is near 56 tan* 11. For the Madergrond path, it is not possible to uniquely identify such a position. Given the 1 slightly undulating nature of the track, many points along the path exist that satisfy the condition -0 = t a n - 1 fi. Though the qualitative shape of each of the speed profiles in Figure 4.8 is similar, different sub-peaks of the profile can contain the maximum speed, given different starting positions. 4.2.1 The Aulta avalanche path and measured events The Aulta avalanche path in Switzerland has a fairly smooth profile, and there are measured speed data with which to compare model simulations. Given the lack of data on flow length, however, a variety of length and depth combinations were attempted. Neither the mean flow depth nor length (recorded at the position of maximum speed) could be held constant while changing the other variable; the two parameters are not independent. In general, the flow depth decreased and the flow length increased from their respective initial values after the simulation began. Upon reaching a stable flow profile, at or near maximum speed, further changes in depth and length were due to changing path geometry. However, some 'memory' of the previous shape of the flow profile was sometimes observed, often in the form of a wave-like disturbance in flow depth along the length of the flow profile. This wave could travel either upslope or downslope, depending on changing slope angle, curvature, width, or friction. Such waves have their primary influence on the dynamics in the form of the internal pressure term. The magnitude of this term is proportional to the gradient of depth between adjacent blocks. A downslope traveling wave had a tendency to increase frontal speed and runout distance, and upslope waves brought the tail of the flow to rest sooner and generally increased the deposition depth. It should be noted that these waves are distinct from shallow surface waves that cannot propagate upstream for Froude numbers greater than 1. D A N does not have the ability to simulate such surface disturbances; the observed waves are due to a changing gradient of flow depth along the total flow length. Figure 4.9 shows the dependence of runout on mean flow depth and flow length. In general, increased depth and length led to increased runout. For mean flow depth approaching 5 m and flow length nearing one third of the path length (~800 m), DAN under-predicts the observed runout by about 40 m. This discrepancy is less than 2% of the observed stopping position, in relation to the total path length traversed. This exercise demonstrates the difficulty that will be encountered if higher-dimensional numerical models are used to predict avalanche runout distances. 57 The maximum speed predicted by the ULE for the Aulta path is 57 m/s, which is less than the maximum measured speed of 65 m/s. This could again be due to the speed scaling with the vertical fall of the path rather than the path length traversed (McClung, 1990). However, given a measured maximum speed that is greater than the ULE prediction, the measured speed was used as the target for the numerical calculation. This was done for the Aulta avalanche of 8 February, 1984, as shown in Figure 4.10. For the ULE maximum speed, it was not possible to reproduce the observed runout distance from any starting position. Using the measured maximum speed of 65 m/s, however, it was possible to nearly reproduce the measured stopping position. The speed profile in the deceleration phase was better represented by the higher maximum speed calculation. Another event for the Aulta path was reported by Salm and Gubler (1985). McClung (1990) used this event to show the sensitivity of his variable friction model to values of surface air drag. Figure 4.11 shows a comparison of the DAN calculation for this path with McClung's calculation including 58 2800 1 — — — —I— — — —I—I—I 1 1 1 1 1 1 1—1 I I I I I I I I I I I Horizontal Reach (m) Figure 4.10: Comparison of measured and DAN speeds, starting in the middle of the path, for the Aulta avalanche of 8 February, 1984- (•), speed measured by Gubler et al. (1986); (- -), DAN calculation for maximum speed 57 m/s, determined by the upper limit envelope for maximum avalanche speed (McClung and Schaerer, 1993); (—), DAN calculation for maximum speed 65 m/s, corresponding with the maximum measured speed. The observed runout was ~2450m. an air drag term. Because the DAN maximum speed is calibrated with actual measured speeds, the calculation implicitly includes the effects of air drag. Thus, the DAN method described here avoids the necessity of selecting and calibrating more than one resistance parameter for numerical simulation. 4.2.2 The Ariefa-Samedan avalanche of January, 1951 An extreme avalanche that has frequently been used to test avalanche dynamics models is the AriefaSamedan event of January, 1951 in Switzerland. Bartelt et al. (1999) used this example to demonstrate the characteristics of their Voellmy-fluid model, which calculates a single flow cross-section from the top of the slope until the flow depth and velocity both drop to zero. The fracture zone limits (upper and lower elevation and width) were documented for the event. Following the Swiss Guidelines, for the purpose of comparison, the initial fracture depth was set at lm. A comparison of results, first starting from rest, is shown in Figure 4.12. The friction coefficient for DAN was set at a constant value fi = 0.38, which, for comparison, reproduced the runout of the Voellmy-fluid model. 59 J 80 T o I I L 1 J r~ 500 I I I I I I I I i T i | i i i i | i—i—i—i—|—i—r*i—r 1500 2000 1000 I 2500 Horizontal Reach (m) Figure 4.11: Comparison of DAN starting in the middle of the path with the model of McClung (1990) for the Aulta path. (—), speed profile measured by Salm and Gubler (1985); (t), DAN speed profile; (O), McClung (1990) speed profile including a surface air drag term. In the sharp deceleration phase DAN closely reproduces the observed speed. Prior to the sharp deceleration, DAN predicts speed slightly greater than measured speed. The most striking feature of the Voellmy-fluid model in this case is the low speed and slow deceleration in the lower reaches of the path. A high value of turbulent friction limits the maximum speed to less than 40 m/s in the middle of the path. Thereafter, the low value of the basal sliding friction coefficient produces a gradual deceleration to rest. This behaviour is inconsistent with observations of the sharp deceleration of snow avalanches in the runout zone. For the Ariefa-Samedan path, DAN calculates a speed profile that is qualitatively more accurate. In the runout zone, the sharp deceleration is consistent with measured speed data, as seen previously in this chapter. For the purpose of defining impact pressure zones for hazard mapping, the Voellmyfluid model applied to a path such as Ariefa-Samedan will give a much less conservative, and likely unrealistic, prediction of speed in the runout zone. An interesting feature of the DAN calculation is in the discrepancy between front and rear maximum speed for this path. The ULE predicts a maximum speed of 46 m/s for this path. The frontal maximum speed predicted by DAN almost reaches this value, but the tail of the flow reaches a maximum speed in excess of 60 m/s. This may not be surprising, however, given a closer look at 60 Horizontal Reach (m) Figure 4.12: Speed profiles for the Ariefa-Samedan avalanche. (- -), maximumflowspeed for the Voellmy-fluid model (Bartelt et al, 1999); (+), DAN frontal speed; (M), DAN tail speed. The thick line is the path profile. the reported limits of the fracture zone for the January, 1951 event (Bartelt et al., 1999). The lower limit of the fracture zone was said to extend down to 2000 m a.s.l. Therefore, the front of the flow travels less than 300 m before reaching the clear transition point at which the slope angle decreases (at about x — 700m). The tail of the flow, however, travels over 3 times as far before the slope angle declines, allowing more time for sustained acceleration. As the tail of the flow catches up to the front, the flow compresses and the flow front is pushed forward. This causes the maximum frontal speed to occur well into the runout zone. Given the recorded fracture limits for this event, the front and tail speed profiles predicted by DAN seem qualitatively appropriate. Starting the DAN calculations in the middle of the path, compared to a simulation beginning from rest, is shown in Figure 4.13. In each case, the runout distance was matched as closely as possible to the January, 1951 avalanche that ran about 1620 m. If the maximum speed predicted by the ULE is used as the maximum speed for the DAN calculation, then the simulation can begin far into the runout zone. Over about the final 300 horizontal metres of the path, the speed for this calculation (X in Figure 4.13) matches the sharp deceleration of the DAN simulation started from rest. If the starting position is chosen such that the position of peak speed occurs at the slope position 61 Horizontal Reach (m) Figure 4.13: DAN frontal speed profiles for the Ariefa-Samedan avalanche of January, 1951. (- -), starting from rest; starting downslope, matching the position of peak speed at the slope position that corresponds to ip — tan" p.; (X), starting further downslope and matching the peak speed with the Upper Limit Envelope prediction. 1 corresponding to tp = tan -1 p, then the speed profile looks quite different. The speed in the runout zone decelerates sharply, but for a given runout zone position, this speed is much less than that calculated by either the simulation started from rest or the simulation started further downslope. The observed stopping position is only reached as a consequence of the tail of the flow catching up to the decelerating front. This briefly contracts the flow length and propels the front of the flow further into the runout zone. This feature is not present in the simulation started further downslope. It is noteworthy that for the Ariefa path, for the calculation with maximum speed located where ip — tan -1 p, the starting position is only about 100 m further downslope than the reported lower limit of the fracture zone. As a result, the DAN simulation does not avoid much uncertainty associated with neglecting entrainment. Additionally, to reproduce the observed stopping position, the maximum frontal speed for this calculation needed to be much higher than the simulation started further downslope. 62 4.3 Examples of D A N Calculations for Surveyed Avalanche Paths Several avalanche paths in the North Cascades mountains in Washington, USA, were surveyed for this analysis. Information collected included slope angle, approximate width of the last extreme avalanche, cross-slope shape, and location of the highway. One avalanche path of interest from a dynamic modelling perspective is Granite Creek 13. This path is characterized by a decrease in path width with distance in the runout zone. Figure 4.14 plots runout distance for different pairs of flow depth and length. For even the smallest size avalanches modelled, the calculated runout distance is well in excess of the location of the highway. The runout distance varies by about 150 m for the parameter pairs tested. Figure 4.14: Runout vs. flow length and meanflowdepth for the Granite Creek 13 avalanche path in the North Cascades, Washington, USA. The highway is located at about 2200 m horizontal reach. Figure 4.15 shows speed profiles that correspond to the upper and lower bound estimates of 63 runout in Figure 4.14. The width profile for the path is also shown. The decreasing width from the avalanche track to the runout zone is responsible for the sustained high speed of the avalanche front, even as the slope angle decreases. The maximum speed, predicted by the ULE, was the same for each event. _ 0 400 100 —I—-J 1 800 1 1— 1200 I I 1600 __l I 2000 i I Horizontal Reach (m) Figure 4.15: Upper- and lower-bound frontal speed profiles for the Granite Creek 13 data in Figure 4-14> well as the path width. The arrow indicates the approximate location of the highway. a s Contrary to the profile of Granite Creek 13, Granite Creek 7 is much wider in the runout zone than in the track. For this path, DAN predicts runout that stops short of the highway for an extreme event, given a friction coefficient selected from the a-angle of the path and a maximum speed predicted according to the ULE of McClung and Schaerer (1993). Figure 4.16 shows upper and lower bound speed profiles for this avalanche path, obtained in a similar fashion to the data in Figure 4.15. 64 Horizontal Reach (m) Figure 4.16: Upper- and lower-bound frontal speed profiles for the Granite Creek 7 avalanche path, as well as the path width. The arrow indicates the approximate location of the highway. 4.4 Discussion The proposed method of starting an avalanche dynamics model in the middle of the path, rather than at the top, is able to reproduce measured runout positions. The qualitative behaviour of the DAN speed calculations agrees with observations in many cases. The uncertainty introduced by this method, which is unique from traditional dynamics models, is mostly related to spatial characteristics of avalanche flow. Many of these characteristics, such as flow length and depth, are approximately constrained by field observations. The sensitivity of runout calculations to variations in these spatial parameters was explored. It is difficult to say how the different uncertainties involved are compounded into the final speed calculation for a given runout position. Currently, many calculations must be carried out with different spatial flow parameters, given constant friction and a given maximum speed, in order to establish upper and lower bounds on predicted speed in the runout zone. As experimental observations and 65 recording standards improve for snow avalanches, the uncertainty in many of the relevant parameters will likely be reduced. It is further believed that spatial characteristics of avalanche flows will be easier to measure, correlate, and implement into a dynamics model than the processes of basal erosion of existing snowcover and subsequent entrainment of mass into the core of a dense flowing avalanche. These basal processes are complex, three-dimensional phenomena that are not likely to be accurately captured by simple linear additions to an existing dynamics model. Though it is desirable from a practical standpoint to build a model that physically corresponds to reality, and is complete in the sense that every element of reality is captured by an element of the model, this is currently not practical. This fact has motivated the pursuit of a modelling method that reduces the unquantifiable uncertainty in neglecting or arbitrarily specifying entrainment. 66 Chapter 5 Conclusions Extreme dry snow avalanches, which are characterized by maximum speed and minimum friction, are the only avalanches that one can currently expect to model. The dynamics of short-returnperiod avalanches are more strongly dependent on ploughing and entrainment of existing snowcover. Application of a dynamics model for such events is inappropriate, especially when the given model is calibrated with extreme event data, as in this study. The typical method of avalanche dynamics modelling begins calculations in the starting zone of an avalanche path. This requires assumptions about entrainment of snow in the mass balance and resistance of the flow. Often entrainment is neglected entirely, and the mass of the flow is held constant from avalanche initiation to deposition. I have asserted that this assumption is inappropriate and thus provides the impetus for an alternative approach. Thus, a new procedure for modelling the speed profiles of extreme avalanches has been developed. It involves starting the calculations of a dynamics model in the middle of an avalanche path with an initial speed as input. This initial speed comes from a conservative upper limit on maximum measured avalanche speeds. It is assumed that entrainment is negligible as the avalanche decelerates, relative to the acceleration phase higher on the slope. By restraining the model domain of interest to the deceleration phase in avalanche runout zones, this procedure avoids dealing with much of the uncertainty introduced by making arbitrary assumptions about the physics or significance of entrainment. A specified avalanche stopping position is necessary in order to scale avalanche speed in the runout zone (McClung, 1990). Given a stopping position and an upper bound flow speed, both determined from appropriate local and regional data, the new modelling procedure produces a profile of avalanche speed for either the flow front, tail, or center of mass. The mass of the flow is thus held constant in the model calculations as the flow decelerates. In 67 the formulation of the model DAN (Hungr, 1995) for this study, both the driving and resistance forces are proportional to mass. Thus, the model sensitivity to other physical parameters such as now depth, length, and width was explored, given a range of reasonable values of the total mass. The flow mass, it should also be noted, is dependent on the density of the flow, which is still somewhat poorly defined. It was therefore seen as more appropriate to explore the sensitivity to flow depth, width and length than to total flow mass. The model uses a constant Coulomb basal friction formulation as the only flow resistance. In reality, speed-dependent surface drag limits the maximum avalanche speed in the middle of the path. However, since the new procedure is calibrated based on measured maximum speeds, this effect is implicitly accounted for. The friction coefficient is chosen either from the alpha angle of the avalanche path (McClung, 1990, Equation (9) with no speed-dependent drag) or from the friction necessary to reproduce the given runout position with a traditional dynamics calculation beginning from rest in the start zone. The results are similar for both methods of selecting the friction coefficient, though selecting the friction based on the path geometry is preferable to having a dependence on the traditional dynamics calculation procedure. Other required inputs for the new modelling procedure are the initial spatial characteristics of the flow. The depth profile and flow length must be input at the start of the calculation. The width of the flow is allowed to vary, but must also be input by the user based on the geometry of the path in question. Of these three flow characteristics, I have shown that the calculations are the most sensitive to changing flow width, so care should be taken to accurately designate the flow width along the avalanche path. I also suggest that variable-width avalanche dynamics calculations should be implemented before introducing additional complexity with a variable friction formulation. It was found that there was redundancy in the selection of the starting position for the model calculations in some cases. Multiple starting positions in the middle of the Madergrond avalanche path, given the same maximum speed, for example, were able to reproduce the same runout distance (see Figure 4.8). The speed profiles in the runout zone for such cases were often quite similar, however, indicating little sensitivity to the selection of where in the middle of the path to start the calculations. The effect of the internal pressure formulation in DAN was explored. Many avalanche dynamics models have adopted a variable internal pressure formulation, analogous to static earth pressure, with 68 active and passive limits. For this study, the internal pressure was held constant at the passive, or compressive, limit. This is assumed to be more appropriate in a typical runout zone as an avalanche decelerates. Figure 3.3 demonstrates that the difference in calculated runout position between a constant passive and variable active/passive internal pressure formulation is small. Figures 3.9 demonstrates that this internal pressure formulation produces excessive and unrealistically rapid thinning of a mass placed on a 0° slope. In this sense, the internal pressure promotes flow thinning if no unique boundary conditions are placed on the front and tail of the mass. If the internal earth pressure analogy were neglected in DAN (i.e. Pi = 0 for all i in Equation 2.4), there would be no internal mechanism to promote flow thinning. The new modelling procedure was able to reproduce the observed sharp deceleration of avalanche frontal speed in runout zones. This is worthy of remark, because most decisions about land use planning and hazard assessment are made in the lower reaches of avalanche paths. Therefore, a dynamics model should reproduce the observed behaviour of avalanches as they decelerate. 5.1 Recommendations Some possible areas for future improvement of observational standards were identified in this study. To date, there has been little attention paid to the longitudinal length of snow avalanches. Aside from order of magnitude estimates (e.g. Perla et al., 1980), I know of no recorded observations of flowing avalanche length. Given the sensitivity of DAN calculations to flow length, this matter should be given further attention. Whether in absolute terms, or scaled with the length of the path, better constraints on the likely length of an extreme avalanche will improve the predictions of multi-dimensional dynamics models. Similarly, the depth profile of the dense core of a flowing avalanche remains largely uncertain. The maximum depth can often be inferred after an event from vegetation damage, compressed snow on the uphill side of trees, or lateral snowcover disturbances in confined paths. However, the location along the profile at which this maximum depth occurs, on average, is unknown and may be quite variable. Laboratory turbid gravity currents may provide some clues as to what the dense core may be shaped like. Descriptions of the 'coherence' of a gravity current and its connection to the shape of the dense profile, such as in Ilstad et al. (2004), may be appropriate avenues for further exploration. 69 This study assumed a constant basal sliding friction coefficient. This is more correctly interpreted as the average value of a friction coefficient that is likely to be a function of speed, and perhaps other factors. A changing frictional formulation, such as in McClung (1990), could be implemented in DAN in a future study. Such a formulation in DAN could lead to better agreement with measured speed for some avalanches in some paths. Caution should be exercised in applying such an analysis, however, while modelling an avalanche path with constant or unit width. The sensitivity of DAN to changing path width, as seen in Chapter 3, should warrant an accurate characterization of the path width before introducing additional complexity in a friction formulation. Sharp, local peaks in measured speed could arise from changes in path confinement, which might be easily captured by a variable-width model with constant friction. The frictional formulation of Louge and Keast (2001) retains the coupling between the shear and normal forces in a granular flow, while describing different components of basal friction. The two primary components account for frictional rubbing between particles and interparticle collisions, with the weight of each component's contribution to the total friction being dependent on the Froude number of the flow. Though such an approach is probably not yet justified for snow avalanches given the large uncertainties involved, it seems an intriguing possibility for future exploration. The method of selection of the friction coefficient for this study could be quantitatively analyzed. This would involve a comparison of a angles for events calculated in DAN with the a angles of the paths, as defined by the previous extreme avalanche. It was not determined in this study whether, on average, DAN had a tendency to over- or underestimate runout if the friction was set according to p = tan a. Theoretically, for a one-dimensional model, setting the friction coefficient according to n = tan a should produce a slide that exactly reaches the point defined by the a angle in the runout zone. Many simulations could be compared to statistically determine the behaviour of DAN relative to this theoretical result. Any deviation may provide a hint as to how to adjust the friction coefficient with speed. However, such a deviation could also be an artifact of the manner in which the initial speed is assigned to the flow, or where exactly along the path a simulation is started. The location along the path at which a simulation is started could be further explored. Starting DAN in the middle of the avalanche path excludes a good portion of the path over which the flow would likely entrain snowcover. The length, or area (given knowledge of the path width), of the path 70 that was excluded from the calculations could be quantified. If it could be claimed that, say, 50% of the path was excluded from the analysis, this would provide more strength to the method proposed in this study than if only, on average, 25% of the path was excluded. The mass of an extreme avalanche can be statistically estimated using regional snowpack and meteorological records. Such a method is outlined in Schaerer and Fitzharris (1984). Given that both the gravitational driving force and frictional resistance are proportional to mass in a Coulomb basal friction model, the mass by itself does not influence the calculations. The spatial distribution of that mass, however, is very important. A statistical prediction of the mass of an extreme avalanche could provide a valuable additional constraint on the range of flow lengths and depths for use in the model. High-quality data sets of spatial and temporal observations of extreme avalanches, against which the proposed method could be further calibrated and refined, are not currently available. Useful data in the future would include the flowing length of the dense avalanche core, the flowing width and depth, the cross-sectional shape of the avalanche path, and speed profiles of the front and tail of the flow. In the meantime, simplifications and assumptions must be made in order to carry out any meaningful calculations. Given the proposed approximations and assumptions for this study, the new modelling procedure using DAN proved to be a promising approach to calculating the speed profiles of extreme snow avalanches. 71 Bibliography Ancey, C. and Meunier, M . 2004. Estimating bulk rheological properties of flowing snow avalanches from field data. Journal of Geophysical Research, 109. F01004, doi:10.1029/2003JF000036. Bagnold, R. A. 1954. Experiments on a gravity-free dispersion of large solid spheres in a Newtonian fluid under shear. Proceedings of the Royal Society of London, (225):49-63. Bartelt, P., Salm, B., and Gruber, U. 1999. Calculating dense-snow avalanche runout using a Voellmy-fluid model with active/passive longitudinal straining. Journal of Glaciology, 45(150):242254. Brugnot, G. and Pochat, R. 1981. Numerical simulation study of avalanches. Journal of Glaciology, 27(95) :77-88. Dent, J. D. 1992. The dynamic friction characteristics of a rapidly sheared granular material applied to the motion of snow avalanches. Annals of Glaciology, 18:215-220. Dent, J. D. and Lang, T. E. 1980. Modeling of snow flow. Journal of Glaciology, 26(94). Gauer, P. and Issler, D. 2004. Possible erosion mechanisms in snow avalanches. Annals of Glaciology, 38:384-392. Goodwin, M . A. and Cowin, S. C. 1971. Two problems in the gravity flow of granular materials. Journal of Fluid Mechanics, 45:321-339. part 2. Gubler, H., Miller, H., Klausegger, G., and Suter, U. 1986. Messungen an Fliesslawinen, Rep. 41. Eidg. Inst, fur Schnee- und Lawinenforsch. Davos, Switzerland. Hungr, O. 1995. A model for the runout analysis of rapid flow slides, debris flows, and avalanches. Canadian Geotechnical Journal, 32:610-623. 72 Hutter, K., Szidarovsky, F., and Yakowitz, S. 1987. Granular shear flows as models for flow avalanches. In Salm, B. and Gubler, H., editors, Avalanche Formation, Movement and Effects. (Proceedings of the Davos Symposium, September 1986) IAHS Publ. no. 162. Ilstad, T., Elverh0i, A., Issler, D., and Marr, J. G. 2004. Subaqueous debris flow behaviour and its dependence on the sand/clay ratio: a laboratory study using particle tracking. Marine Geology, 213:415-438. Jenkins, J. T. and Savage, S. B. 1983. A theory for the rapid flow of identical, smooth, nearly elastic, spherical particles. Journal of Fluid Mechanics, 130:187-202. Kotlyakov, V. M . , Rzhevkiy, B. N., and Samoylov, V. A. 1977. The dynamics of avalanching in the Khibins. Journal of Glaciology, 19:431-439. Lang, T. E., Nakamura, T., Dent, J. D., and Martinelli, Jr., M . 1985. Avalanche flow dynamics with material locking. Annals of Glaciology, 6. Louge, M . Y . and Keast, S. C. 2001. On dense granular flows down flat frictional inclines. Physics of Fluids, 13(5):1213-1233. McClung, D. and Schaerer, P. 1993. The Avalanche Handbook. The Mountaineers, Seattle, WA. McClung, D. M . 1990. A model for scaling avalanche speeds. Journal of Glaciology, 36(123):188-198. McClung, D. M . 2001a. Extreme avalanche runout: A comparison of empirical models. Canadian Geotechnical Journal, 38(6):1254-1265. McClung, D. M . 2001b. Superelevation of flowing avalanches around curved channel bends. Journal of Geophysical Research, 106(B8):16,489-16,498. McClung, D. M . and Mears, A. I. 1991. Extreme value prediction of snow avalanche runout. Cold Regions Science and Technology, 19:163-175. McClung, D. M . and Mears, A. I. 1995. Dry-flowing avalanche run-up and run-out. Journal of Glaciology, 41(138):359-372. McClung, D. M . and Schaerer, P. A. 1983. Determination of avalanche dynamics friction coefficients from measured speeds. Annals of Glaciology, 4:170-173. 73 McClung, D. M . and Schaerer, P. A. 1985. Characteristics of flowing snow and avalanche impact pressures. Annals of Glaciology, 6:9-14. Mellor, M . 1968. Avalanches. Technical Report CRSE III-A3d, Cold Regions Research & Engineering Laboratory. Naaim-Bouvet, F., Naaim, M . , and Faug, T. 2004. Dense and powder avalanches: momentum reduction generated by a dam. Annals of Glaciology, 38:373-378. Norem, H., Irgens, F., and Schieldrop, B. 1987. A continuum model for calculating snow avalanche velocities. In Salm, B. and Gubler, H., editors, Avalanche Formation, Movement and Effects. (Proceedings of the Davos Symposium, September 1986) IAHS Publ. no. 162. Perla, R., Cheng, T. T., and McClung, D. M . 1980. A two-parameter model of snow-avalanche motion. Journal of Glaciology, 26(94): 197-207. Salm, B. 1966. Contribution to avalanche dynamics. In International Association of Scientific Hydrology Publ. No. 69, pages 199-214. Symposium at Davos, Switzerland, 5-10 April 1965. Salm, B. 1993. Flow, flow transition and runout distances of flowing avalanches. Annals of Glaciology, 18:221-226. Salm, B. and Gubler, H. 1985. Measurement and analysis of the motion of dense flow avalanches. Annals of Glaciology, 6:26-34. Sartoris, G. and Bartelt, P. 2000. Upwinded finite difference schemes for dense snow avalanche modeling. International Journal for Numerical Methods in Fluids, 32:799-821. Savage, S. B. and Jeffrey, D. J. 1981. The stress tensor in a granular flow at high shear rates. Journal of Fluid Mechanics, 110:255-272. Schaer, M . and Issler, D. 2001. Particle densities, velocities and size distributions in large avalanches from impact-sensor measurements. Annals of Glaciology, 32:321-327. Schaerer, P. A. 1974. Friction coefficients and speed of flowing avalanches. In IAHS-AISH Publ. No. 114, pages 425-432. Union Geodesique et Geophysique Internationale. Association Internationale des Sciences Hydrologiques. Commission des Neiges et Glaces. Symposium, Mecanique de la niege, Actes du colloque de Grindelwald, Avril 1974. 74 Schaerer, P. A. and Fitzharris, B. B. 1984. Estimation of the mass of large snow avalanches. Canadian Journal of Civil Engineering, 11:74-81. Schaerer, P. A. and Salway, A. A. 1980. Seismic and impact-pressure monitoring of flowing avalanches. Journal of Glaciology, 26(94): 179-187. Sovilla, B., Sommavilla, F., and Tomaselli, A. 2001. Measurements of mass balance in dense snow avalanche events. Annals of Glaciology, 32:230-236. Tiefenbacher, F. and Kern, M . A. 2004. Experimental devices to determine snow avalanche basal friction and velocity profiles. Cold Regions Science and Technology, 38:17-30. Tomasson, G. G. and Johannesson, T. 1995. A simple dynamic model for computation of snow avalanche runout, unpublished. Voellmy, A. 1955. Uber die Zerstorungskraft von Lawinen. Schweiz. Bauztg., 73:159-165, 212-217, 246-249, 280-285. Wilcox, D. C. 2003. Basic Fluid Mechanics. DCW Industries, 2nd edition.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Dynamic modelling of extreme speed profiles of dry...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Dynamic modelling of extreme speed profiles of dry flowing avalanches Borstad, Christopher Paul 2005
pdf
Page Metadata
Item Metadata
Title | Dynamic modelling of extreme speed profiles of dry flowing avalanches |
Creator |
Borstad, Christopher Paul |
Date Issued | 2005 |
Description | A brief review of different approaches to avalanche dynamics modelling is considered, drawing on observations and knowledge of avalanche behaviour and analogies with other physical phenomena. The numerical model used for this analysis, DAN, is described with respect to the formulation of the equations of momentum and continuity. Avalanche path geometry is approximated as a rectangular channel with the possibility of variable width. The flow mass is discretized longitudinally into a variable number of constant-volume blocks. The uncertainty surrounding the dynamics of frontal ploughing of undisturbed snow, basal erosion, and mass entrainment is discussed. In order to avoid this uncertainty, which is implicit when modelling a slab avalanche beginning from rest in the starting zone, the calculations in DAN begin near the middle of the slope where the flow is expected to achieve maximum speed. Knowledge of the maximum probable speed for an extreme avalanche is used to supply an initial, non-zero speed to the flow near the middle of the path, and entrainment is assumed negligible as the flow decelerates. A sensitivity analysis explores theoretical considerations surrounding the modelling of a flow mass in multiple dimensions, given this new procedure as well as traditional from-rest approaches. The model is shown to be more sensitive to changes in flow width than flow length or depth. DAN calculations are compared against examples of recorded extreme avalanche data, with generally good agreement. The new modelling procedure was able to reproduce the observed sharp deceleration of avalanche frontal speed in runout zones. This is worthy of remark, because most decisions about land use planning and hazard assessment are made in the lower reaches of avalanche paths. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-05 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0063320 |
URI | http://hdl.handle.net/2429/17501 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2006-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2006-0011.pdf [ 4.62MB ]
- Metadata
- JSON: 831-1.0063320.json
- JSON-LD: 831-1.0063320-ld.json
- RDF/XML (Pretty): 831-1.0063320-rdf.xml
- RDF/JSON: 831-1.0063320-rdf.json
- Turtle: 831-1.0063320-turtle.txt
- N-Triples: 831-1.0063320-rdf-ntriples.txt
- Original Record: 831-1.0063320-source.json
- Full Text
- 831-1.0063320-fulltext.txt
- Citation
- 831-1.0063320.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0063320/manifest