Introduction

In North America and Europe, recreational, residential and commercial use of alpine terrain has required an evolution of avalanche control techniques in order to limit or mitigate the damage caused by avalanches. These techniques consist of a variety of procedures. One of the more active measures involves the initiation of avalanches with explosives, thereby either avoiding large avalanches or reducing the potential of the slope to avalanche unexpectedly. This technique is widely used in ski areas and over threatened highways. Over 100000 explosive charges are detonated annually for avalanche control (Reference PerlaPerla, 1978). The type of explosive, size of the charge and the style of delivery are primarily determined by trial and error. Information leading to a more efficient determination of how to employ explosives in releasing avalanches is thus desirable. When an explosive charge is detonated in or closely above the snowpack, inelastic deformation in the immediate vicinity of the explosion takes place, forming a crater and a region of cracks surrounding the crater. Outside of this region, inelastic deformation of the snow becomes insignificant, but attenuation of the stress waves still occurs through geometric spreading, porous structure effects and effects of ice structure/pore material interaction (Reference BrownBrown, 1980a, b). It is the propagation of stress waves in the region outside the crater which this investigation addresses.

The continuum theory of mixtures is based on rational thermodynamics. The fundamental premise is that the constituents making up the mixture can be modeled as superimposed continua such that each point in the mixture is occupied by a material point of each constituent. Balance statements are then postulated for each constituent. They contain additional supply terms to account for the interactions between the constituents. In the last 20 years much effort has been devoted to preserving the generality of the theory, in particular with respect to the constitutive equations. This generality has resulted in complex and cumbersome constitutive relations which inhibit the application process. Reference Bedford and DrumhellerBedford and Drumheller (1983) cite this as one of the factors why applications of the theory have been so slow in developing. There is not a single comparison of theoretical results with experimental data in extensive reviews by Reference Atkin and CraineAtkin and Craine (1976) and Reference Bowen and EringenBowen (1976).

Applications of the continuum theory of mixtures to model snow has likewise received little attention. Reference Decker and BrownDecker and Brown (1985) and Reference Liston, Brown and DentListon and others (1993) used it to model turbulent flow of a mixture of snow and air. Reference Adams and BrownAdams and Brown (1989, 1990) applied a continuum mixture theory to model physical phenomena associated with snow metamorphism. There is little evidence of modern mixture theory being applied to the propagation of stress waves in snow. Reference Bowen and ReinickeBowen and Reinicke (1978) have investigated general wave fronts, shock waves and harmonic plane waves in binary mixtures of elastic materials, but these investigations have remained on the theoretical level. Reference AtkinAtkin (1968) examined the propagation of small amplitude waves for a fictitious mixture of an isotropic elastic solid and an inviscid fluid with an early version of a continuum mixture theory. Reference JohnsonJohnson (1978) investigated elastic stress waves in snow using an early porous-media model developed by Reference BiotBiot (1956). Biot’s theory was specifically formulated to model the linear elastic behavior of a fluid-saturated porous media rather than being a specialized case of a more general mixture theory. While his work was a significant contribution to porous-media modeling there are certain technical problems with it (Reference Bowen and EringenBowen, 1976). Specifically, the constitutive relations for the momentum supplies do not satisfy certain thermodynamic constraints.

Theory

Balance of mass requires that mass be accounted for at all times. This holds for each constituent and for the complete mixture. In this case we consider the mixture to consist of the elastic ice phase and the air and vapor phase. For the a^th constituent the local form of the mass-balance principle may be expressed as

and for the mixture

In the above, ĉ_a is the mass supply which determines the rate at which that particular constituent is acquiring mass from the other constituents. In this study, ĉ_a is assumed to be zero, as mass exchange is not considered. In the above, ρ is the material density of the mixture, ρ_a is the density of the a^th constituent, and ẋ and x’_a are the velocity of the mixture and the a^th constituent, respectively.

The second balance principle is the linear-momentum balance principle. It requires that linear momentum for each constituent of the mixture and also of the complete mixture obey Newton’s second law. In local form, this principle may be expressed as

and for each constituent and for the mixture

is the momentum supply which represents the exchange of linear momentum which occurs as the constituents interact with each other. Τ is the stress tensor, and the terms following the density in the above two equations are the accelerations.

Summing the constituent momentum-balance equations and requiring that they sum to the mixture-balance relation leads to the following equations:

These are the balance equations that are used to model stress wave propagation through snow. The term u_a is the diffusion velocity, which is the velocity of the a^th constituent relative to the mixture velocity. In addition to these, constitutive relations for each constituent are needed to close the problem. The ice phase is modeled as a linear elastic material, since interest here is in wave propagation beyond the crater where inelastic effects dominate. The air and vapor phases are modeled together as a single ideal elastic fluid. The subscripts g and s are labels for the gas and the solid phase, respectively, and the subscript R denotes the reference state. The derivations of the constitutive equations are both tedious and complex, adding little insight to the application process. Thus they will simply be presented. For a rigorous development of the constitutive equations presented here, the reader is referred to Reference Bowen and EringenBowen (1976). The linearized version of the constitutive equations are

w_g and w_s are, respectively, the displacement components for the gas and solid phases. E _s and E _gare the infinitesimal strain tensors for the solid and gas phases, respectively. I is the identity tensor; Π_gR and Π_sR are constants which represent the pressure in each constituent when there is no strain measured relative to the reference state. The coefficients σ_sg, σ_gs, λ_gs, μ_s, λ_g, λ_B are material constants, and σ_sg and σ_gs are coupling coefficients which arise because of a dependence of the momentum supplies on the strain gradients. They account for the local interaction of the solid and the fluid that occurs even in static situations. The term (λ_B − σ_gs ) is analogous to the lame parameter λ in elasticity. The coefficient is the shear modulus for the solid phase, and (σ_sg + λ_gs ) accounts for the dependence of the stress in the solid due to the strain of the fluid, (σ_gs + λ_gs ) accounts for the stress in the fluid due to the strain of the solid. The term (λ_g − σ_sg ) is related to the modulus of elasticity for the fluid. The Stokes drag coefficient, ξ, arises due to the momentum supply dependence on the relative velocities of the solid and gas phases. It is a result of the second law of thermodynamics that ξ > 0, λ_g > 0, μ_s > 0 and

(Reference Bowen and EringenBowen, 1976).

Substituting the constituent equations ((7)-(9)) yields the isothermal linear equations of motion for the solid and gas phases. The resulting system of equations can be decomposed into dilational and rotational components by using a Helmholtz decomposition of each constituent displacement vector, w_a (Reference AtkinAtkin, 1968).

where ϕ_a is a scalar function of time and space, and ψ_a is a vector function of time and space. Physically the dilational part represented by ϕ_a corresponds to a change in volume (i.e. the displacement due to the hydrostatic stresses). Compression of a unit cube would be an example of this. The rotational part represented by ψ_a describes the displacement due to the deviatoric stresses, thus the same unit cube would go through a rotation and distortion without undergoing a change in volume. In harmonic wave motion the displacement due to the dilational disturbance is parallel to the wave propagation direction, while the displacement due to the rotational disturbance is perpendicular to the direction of the wavefront motion. Substituting into the balance equations and ignoring body, the force results in the following system of equations for the solid and gas constituents:

Presently values of the coupling coefficients σ_gs and σ_sg are not known, and therefore it is necessary to consider the special case of σ_gs = σ_sg = 0. This constitutive assumption is known as an ideal mixture and corresponds to assuming that the equilibrium free energy of the a^th constituent is independent of the deformation of the other constituents. At the very least this assumes that the strains in one constituent have an insignificant effect on the stresses in the other constituents (Reference Bowen and EringenBowen, 1976, p. 99).

The stress tensor for an ideal elastic fluid has the form Τ = Π(ρ)I, where Π is the thermodynamic pressure and is a function of density only. An equation of state defines the thermodynamic pressure Π. From mixture theory the isothermal linear equation of state for an elastic fluid (Reference Bowen and EringenBowen, 1976) is

The coefficient λ_g is equivalent to the bulk modulus of elasticity for a linear isothermal elastic fluid.

Tables 1, 2, and 3 show the calculated values for the material constants. The fluid properties exhibit small changes with changes in density for the range considered in this study, and the thermal effects are ignored. Therefore the fluid properties are approximated as constants.

It remains to determine values for the friction coefficient ξ. Reference BiotBiot (1956) writes this coefficient as

where b is the ratio of the total frictional force between the fluid and the solid, per unit volume of bulk material, to the average fluid velocity relative to that of the solid in steady-state flow (Poiseuille flow) (Reference Deresiewicz and RiceDeresiewicz and Rice, 1962). F(k) is a frequency-dependent correction factor which is a measure of the deviation from Poiseuille flow friction, k is a nondimensional frequency parameter, and a is a characteristic size of the pore cross-section, ω is the circular propagation frequency given by ω = 2πf = 2π/Τ, where Τ is the period.

The function F(k) is written as

in which

Here berk and beik denote the real part and the imaginary part of the Bessel-Kelvin function of the first kind and order zero, and ber’k and bei’k are the derivatives of berk and beik with respect to the variable k. For pores which behave like circular tubes of diameter d, b may be written as b = 32μϕ/d² (Reference BiotBiot, 1956). For small values of k (k ≪ 1), F(k) reduces to a form having real and imaginary parts F_R and F_I (Reference Deresiewicz and RiceDeresiewicz and Rice, 1962):

For large values of k(k ≫ 1), F(k) takes on the approximation (Reference Deresiewicz and RiceDeresiewicz and Rice, 1962)

For the purposes of application, a harmonic propagating disturbance is considered. For a dilation wave this is represented by the scalar ϕ_a:

n, η and ω are the unit propagation vector, the circular wave number, and the circular frequency, respectively. In general the wave number η is complex with real and imaginary parts η_R and η_I . It is also dependent on properties of both constituents. The frequency ω is assumed to be real valued. Substituting the complex form of η in Equation (23) yields

This shows that the amplitude of the wave attenuates and that the attenuation is proportional to the imaginary part of the wave number. Thus η_I is called the attenuation coefficient. The propagation velocity is given by c = ω/η_R.

A dispersive media is one in which the wave propagation velocity is a function of frequency. The relationship which yields this frequency-dependent velocity function is termed the dispersion relationship. Thus determination of η is required in order to evaluate the attenuation and the propagation velocity. The attenuation coefficient and circular frequency are found by substituting Equation (24) into the momentum balance equations ((11) and (12)) to obtain a characteristic equation which can be solved to yield a solution for η in terms of the circular frequency ω. The solution procedure is restricted to two frequency ranges (ω < 300 Hz and ω > 2000 Hz) in order to simplify the solution of this characteristic equation. Actually the low-frequency range is of more interest, since stress waves generated by explosives quickly lose their high-frequency content. This implies that by the time the stress levels have attenuated to produce an elastic wave, the frequency content is quite low, usually less than 300 Hz.

However, we will here consider both high-frequency and low-frequency regimes.

In each frequency range two dilatational solutions emerge. These are referred to as dilatational waves of the first and second kind. In the low-frequency range, one of the solutions is of a diffusive nature and does not represent a true wave, while the second solution is wave-like in nature: we will discuss both of these low-frequency dilatational waves. Both solutions are also found in the high-frequency range, and both of these waves are found to behave like true waves.

This same procedure is repeated for the shear wave, which can be defined in terms of the vector ψ defined in Equations (10), (12) and (14). Assuming a form similar to Equation (24), substituting in Equations (12) and (14), and following a procedure as before, results in only one solution which is both damped and dispersive. Therefore the solution of the wave-propagation problem yields three types of high-frequency waves, two dilatational waves and one shear wave, all of which are damped and dispersive. In the low frequency range two types of waves are found (one dilatational and one shear), and these are also damped and dispersive.

Results

Figure 1 shows a good comparison between wave-velocity values from Reference SmithSmith (1965) and theoretically calculated propagation velocities. Smith’s data for the elastic moduli of snow were determined by sonic-wave techniques which correspond to low-amplitude high-frequency waves. To make a comparison, Smith’s elastic moduli were used in the high-frequency limit of the theoretical wave speeds. These values were then compared with the measured sonic-wave velocities. Note here that the first dilation-mode wave does not show up in the experimental data. This is due to the fact that the first dilation mode is highly attenuated (which will be shown in the following figures) and therefore is difficult to detect with measurement techniques.

Fig.1 . Experimental and theoretical propagation velocities.

+ experimental longitudinal wave velocity

− theoretical longitudinal wave velocity

* experimental shear wave velocity

− theoretical shear wave velocity

— ρ = 410 ∼ ρ=440 — ρ = 508 — ρ = 551

One should not read more into Figure 1 than it deserves. The excellent fit between theory and data is due to the fact that the material coefficients calculated by Reference SmithSmith (1965) to match his wave-speed data were used in the solutions developed here. As a consequence, one would expect the theory presented in this paper to give wave speeds very close to those measured by Smith.

Figure 2a and b show velocities as a function of frequency for the range of densities considered in the low- and high-frequency regimes. This figures indicates that the propagation velocity changes with frequency and density.

Fig.2 . Propagation velocity for the first dilation mode, a, in the low-frequency regime; b, in the high-frequency regime.

The asymptotic behavior accurately corresponds to that predicted by the limiting cases of low and high frequency considered earlier. Figure 3 shows that the velocity for the first dilation mode approaches its high-frequency asymptote from above. This behavior is characteristic of partial differential equations which describe a coupling of scalar or vector wave fields (Reference AtkinAtkin, 1968)

Fig.3 . Attenuation coefficient for the first dilation mode. a, in the low-frequency regime; b, in the high-frequency regime.

The wave velocities were found not to be strongly dependent on frequency for all wave types in the high-frequency range, whereas the frequency dependence was more pronounced in the low-frequency range (less than 300 Hz). This was the general trend for all wave types, so wave velocity will not be discussed further in order to concentrate on attenuation characteristics. It suffices to say that velocities for the second dilation mode are greater than those for the first dilation mode. This is in agreement with the “slow” and “fast” dilation waves found in the porous-media models of Reference BiotBiot (1956) and that of Reference AtkinAtkins (1968).

Figure 3a and b show the attenuation coefficients for the first dilatation mode in the low- and high-frequency ranges, and Figure 4a and b show attenuation coefficients for the second dilation mode. Results indicate that the attenuation coefficient increases with increasing frequency and decreases with increasing density. Attenua-tion coefficients for the shear wave are shown in Figure 5a and b. These plots indicate that the attenuation coefficient also increases with increasing frequency and decreases with increasing density. Comparison with the three waves shows that the first dilatational wave is the most highly attenuated of the three waves.

Fig.4 . Attenuation coefficient for the second dilation mode. a, in the low-frequency regime; b, in the high-frequency regime.

Fig.5 . Attenuation coefficient for the shear wave, a, in the low-frequency regime; b, in the high-frequency regime.

Figures 6 and 7 show the normalized amplitude attenuation of all three waves for three different frequencies. In Figure 6 the dilation wave modes are plotted. Note that the first dilation mode is highly attenuated relative to the second mode. For ω = 600 rads s⁻¹ the amplitude of the first dilation mode decreases very quickly to nearly zero and thus does not show up on the plot. Also the amplitude for the second dilation mode for the frequencies l0rads s⁻¹ and 100 rads s⁻¹ are so similar that they appear as one line on the plot.

Fig.6 . Normalized amplitude attenuation for the two dilation waves with density of 410kgm⁻³.

— 1st dilation mode at 10 rads/s

— 1st dilation mode at 100 rads/s

— 2nd dilation mode at 10 rads/s

•••• 2nd dilation mode at 100 rads/s

2nd dilation mode at 600 rads/s

Figure 7 shows the amplitude attenuation for the shear wave. All three waves exhibit increasing amplitude decay with increasing frequency, the first dilation mode being the most attenuated and the shear wave being the least attenuated.

Fig.7 . Normalized amplitude attenuation for the shear wave with density of 410kgm⁻³.

— 10 rads/s ···· 100 rads/s — 600 rads/s

Summary and Discussion

In a porous media consisting of a linear elastic solid and an elastic fluid, two dilational modes and a shear wave can propagate. The first dilational wave mode is associated with the fluid in the pore space while the second dilational wave mode and the shear wave are associated with the solid elastic material. The first dilation wave mode associated with the fluid phase is characterized by a diffusive nature at low frequencies, nondispersive character at high frequency, slow wave speed and high attenuation. The second dilation wave mode exhibits little attenuation with nondispersive character in both the low-frequency and high-frequency limits. Propagation velocity for the second mode is the fastest of all three waves. The shear wave travels faster than the slow dilation wave but slower than the second dilation wave. In the low- and high-frequency limits the shear wave propagates without dispersion. Like the second dilatation wave, the shear wave is not highly attenuated.

Wave velocities of the three wave types increase with increasing frequency and increasing density, while the attenuation coefficients increase with increasing frequency and decrease with increasing density for all three waves. High-frequency theoretical velocities show good agreement with sonic-wave velocities, and attenuation coefficients exhibit behavior that is known to exist in snow.

As the model presented here is a simple representation of a complex physical phenomena, it could be expanded in several directions. More experimental data is required to obtain values for the coupling coefficients which would allow the stress of one constituent to depend on the strain of the other constituent. Existing data for snow is primarily confined to the higher densities. While this may be adequate for maritime and polar climates, it is inappropriate for continental climates where the snow-pack has much lower densities. Formulating the governing equations in spherical coordinates would allow geometric spreading to influence the behavior of the model and would be a more realistic approximation of the wave behavior which occurs in avalanche control. This investigation modeled the void spaces as parallel cylindrical tubes; a more complex representation of the void spaces would allow void spaces to be sinuous or to take on different cross-sectional shapes or diameters. Scattering effects due to the porous geometry would then play a greater role in the attenuation characteristics.

The results reported here do allow some conclusions regarding explosives characteristics which would be desirable for effective control of avalanches. One obvious finding is that “slow” explosives (low-frequency) may be more effective than “fast” explosives (high-frequency), since wave amplitude attenuates more quickly as the frequency increases. Waves with a low frequency would tend to propagate further into a snow cover than a high frequency wave. However, it should be kept in mind that fast explosives usually contain much more energy and produce significantly higher pressures, thereby negating the advantages of slower blasts. The trick would be to use slow explosives that can also generate the needed blast energy. These results may explain the observed effectiveness of aerial blasts and the gas-operated control measures developed in Europe. Both of these techniques deliver intense, low-frequency loads to the snow-cover surface, which allow the pressure waves to penetrate significant distances into the snow cover. It is not clear to us that this would translate to stress waves in wet snow, as our observations have been restricted to dry snow, and the advantages of slow explosives in wet snow may not be significant. More research needs to be done on this problem.

Acknowledgement

The work reported here was funded by the Army Research Office under grant no. DAAL 03-92-G-0310. The authors wish to express their appreciation for ARO’s support.