Dynamic models for impact-initiated stress waves through snow columns

Abstract

The objective of this research is to model snow's response to dynamic, impact loading. Two constitutive relationships are considered: elastic and Maxwell-viscoelastic. These material models are applied to laboratory experiments consisting of 1000 individual impacts across 22 snow column configurations. The columns are 60 cm tall with a 30 cm by 30 cm cross-section. The snow ranges in density from 135 to 428 kg m⁻³ and is loaded with both short-duration (~1 ms) and long-duration (~10 ms) impacts. The Maxwell-viscoelastic model more accurately describes snow's response because it contains a mechanism for energy dissipation, which the elastic model does not. Furthermore, the ascertained model parameters show a clear dependence on impact duration; shorter duration impacts resulted in higher wave speeds and greater damping coefficients. The stress wave's magnitude is amplified when it hits a stiffer material because of the positive interference between incident and reflected waves. This phenomenon is observed in the laboratory and modeled with the governing equations.

1. Introduction

Understanding how snow responds to dynamic, impact loading has relevance to polar aircraft landing (Napadensky, 1964), cryo-seismological monitoring (Johnson and others, 1993; Podolskiy and Walter, 2016), snowplowing (Wakahama and Sato, 1977), vehicular travel on snow (Brown, 1979a), avalanche release (Schweizer and others, 1995; Thumlert and Jamieson, 2014) and other cold regions problems. Although the viscoelastic nature of snow has long been recognized, many modeling efforts treat snow as an elastic material. In this paper, we observe how columns of snow respond to impacts experimentally. Then, we compare elastic and viscoelastic models of the snow's behavior. The chief motivator for this work is avalanche release, however the modeling techniques extend to those applications noted.

One of the first studies exploring snow's response to impact involved explosively driving a metal plate into a column of dense snow (~500 kg m⁻³) and observing the deformation with a streak camera (Napadensky, 1964). The observed stress waves were later modeled using a constitutive law based on pore collapse (Brown, 1979b, 1980a) and a different constitutive law based on neck growth (Brown, 1980b, 1980c, 1980d).

Stress waves generally attenuate through two mechanisms: material and geometric (Kolsky, 1963). Material damping is the energy loss due to internal friction, heat generation, plastic deformation, cracking, internal reflections from inherent heterogeneity, etc., and involves the material absorbing some of the stress wave's energy. On the other hand, geometric damping is attributed to the stress wave expanding and losing energy density as its surface area grows. Attenuation of stress waves in one dimension is entirely attributed to material damping. Geometric damping is relevant in two and three dimensions but is out of the scope of this paper.

Other studies have focused on how acoustoelastic signals transmit through snow. Some early studies (Ishida, 1965; Lang, 1976) were focused on sonic attenuation. Others involved the application of Biot's model (Biot, 1956) which describes dilatational waves transmitting through both the air and ice, as well as a distortional wave through the ice (Johnson, 1982; Capelli and others, 2016). This group of studies generally involved a wave source of insufficient magnitude to cause snow failure (e.g. breaking pencil lead) but provides estimates of snow's elastic modulus as was done by Gerling and others (2017).

In addition to elasticity, the viscous nature of snow has long-been observed (Bader and others, 1939). The Burger's linear, four-element viscoelastic model has been applied (Yosida and others, 1956; Shinojima, 1967) and nonlinear effects have been included (Bader, 1962). Regardless of the particular constitutive relationship, snow's viscoelastic properties have been shown to vary across orders of magnitude (Shapiro and others, 1997). Thus, any attempt to characterize and model snow's behavior should be at a relevant scale to the problem at hand.

Motivated by avalanche release, a variety of field-based studies have measured the snow's response to dynamic loads in naturally occurring mountain snowpacks. In these experiments, the most controlled and repeatable loading methods have involved dropping a known weight from a known height as performed with the Rammrutsch (Schweizer and others, 1995) and a drop hammer (Thumlert and Jamieson, 2015). Additionally, stability test loading schemes have been implemented to impact the snow. These include hand taps during compression tests/extended column tests (van Herwijnen and Birkeland, 2014; Thumlert and Jamieson, 2015; Griesser and others, 2023) and Rutschblock loading (Schweizer and others, 1995; Schweizer and Camponovo, 2001). More realistic avalanche triggers and less repeatable loading methods have included skiing with a knee dip, skiing and falling and snowmobiling (Thumlert and others, 2012; Thumlert and Jamieson, 2014). Observations in these studies have included the use of a camera with markers in the snow pit's side wall (Schweizer and others, 1995; van Herwijnen and Birkeland, 2014) and force sensors within the snowpack (Schweizer and others, 1995; Schweizer and Camponovo, 2001; Thumlert and others, 2012; Thumlert and Jamieson, 2014, 2015; Griesser and others, 2023). Some of these studies (Schweizer and Camponovo, 2001; Thumlert and Jamieson, 2014, 2015) modeled the stress within the snow by idealizing the snow as an elastic, semi-infinite half space under a line load. Although force was measured dynamically, singular peak force measurements were drawn from these timeseries data and fit to a static model. These studies motivate our work to model the dynamic nature of stress wave transmission tuned with experiments executed in a laboratory-controlled environment.

2. Theoretical background

In an elastic model, stress, σ, is directly related to strain, $\epsilon$, by only the elastic modulus, E. Negative stress is compressive, and positive stress is tensile.

(1)$$\sigma = E{\rm \epsilon }$$

In a Maxwell-viscoelastic model, the constitutive relationship depends on the stress rate, $\dot{\sigma } = {\rm d}\sigma /{\rm d}t$, strain rate, $\dot{\epsilon } = {\rm d}\epsilon /{\rm d}t$, and viscosity, η, in addition to the elastic modulus, E, and current applied stress, σ. It takes the form,

(2)$$\dot{\sigma } + \displaystyle{E \over \eta }\sigma = E\dot{{\rm \epsilon }}{\rm \;}$$

The equations of motions for these materials are derived by inserting their respective constitutive equations into the balance of linear momentum. An elastic material's equation of motion in 1D is referred to as the wave equation (Graff, 1975). This formulation neglects body forces.

(3)$$\displaystyle{{\partial ^2w( {z, \;t} ) } \over {\partial t^2}} = c^2\displaystyle{{\partial ^2w( {z, \;t} ) } \over {\partial z^2}}$$

where w is the displacement in the z direction, t is time and c is the elastic-dilatational wave speed, simply referred to as wave speed hereafter. In this study, z is defined as positive in the upwards direction. The wave speed depends on the elastic modulus, E, and density, ρ.

(4)$$c = \sqrt {\displaystyle{E \over \rho }} $$

The equation of motion for a Maxwell-viscoelastic material in 1D is referred to as the damped wave equation (Davis, 2000).

(5)$$\displaystyle{{\partial ^2w( {z, \;t} ) } \over {\partial t^2}} + \displaystyle{E \over \eta }\displaystyle{{\partial w( {z, \;t} ) } \over {\partial t}} = c^2\displaystyle{{\partial ^2w( {z, \;t} ) } \over {\partial z^2}}$$

The additional term in this equation determines the attenuation of displacement. So long as the coefficient, E/η, is positive, the magnitude of displacement will decrease as the wave travels, effectively dampening it. Notice that as the viscosity increases (η → ∞) the material approaches a purely elastic material, and the damped wave equation reduces to the elastic wave equation.

The Maxwell-viscoelastic wave speed is the same as that of purely elastic wave speed (Davis, 2000). This wave speed, c, is not to be confused with the particle velocity, ∂w/∂t. In a homogeneous continuum, the wave speed is constant, whereas the particle velocity varies in space and time as the wave passes through. These governing equations are exemplified in the Supplementary Material Figure S1.

3. Methods

3.1. Modeling

In the spirit of developing the simplest, effective model, a linear-elastic model is first evaluated followed by the two-element, Maxwell-viscoelastic model. The elastic model contains no mechanism of energy dissipation, whereas the viscoelastic model does.

3.1.1. General problem formulation

In the model domain, the snow has a cross-sectional area, A, density, ρ, elastic modulus, E, viscosity, η, and extends from z = 0 at the base to z = H at the top (Fig. 1). The interior of the domain is governed by the wave equation in the elastic case and the damped wave equation in the Maxwell-viscoelastic case. Both partial differential equations are second order in space and time, so two initial conditions and two boundary conditions are needed for a solution.

Figure 1. The modeled domain of the laboratory tests. A time-dependent load −F(t) is applied to the top of the snow column. The snow has a cross-sectional area, A, density, ρ, elastic modulus, E, viscosity, η, and extends from z = 0 at the base to z = H at the top. The base is modeled as a spring with spring constant, k _eff.

There is no displacement, w, and no velocity, ∂w/∂t, when t = 0, leading to the two initial conditions:

(6)$$w( {z, \;0} ) = 0\;$$

(7)$$\displaystyle{{\partial w( {z, \;0} ) } \over {\partial t}} = \dot{w}( {z, \;0} ) = 0$$

The top boundary is subject to a force function, F(t), which is idealized as a Gaussian, which represents a smooth, continuous pulse of disturbance.

(8)$$F( t ) = F_{\,peak}{\rm e}^{{-}1/2{( {( {t-t_{\,peak}} ) /( {d/6} ) } ) }^2}$$

F _peak is the peak force and t _peak is the time at which the peak force occurs. Since 99.7% of the curve's magnitude occurs during six std dev. of a Gaussian, the duration of the force curve, d, is in the denominator of the exponent and divided by six. For example, a 12 ms impact is modeled as a Gaussian with a std dev. of 2 ms. A visualization of this idealization is shown in Supplementary Material Figure S2. This force function is translated into a boundary condition by dividing the cross-sectional area to form a stress function and applying the appropriate constitutive relationship.

In the elastic case, the top boundary condition is:

(9)$$\displaystyle{{F( t ) } \over A} = E\displaystyle{{\partial w( {H, \;t} ) } \over {\partial z}}$$

In the viscoelastic case, the top boundary condition is

(10)$$\displaystyle{\partial \over {\partial t}}\displaystyle{{\partial w( {H, \;t} ) } \over {\partial z}} = \displaystyle{1 \over A}\left({\displaystyle{{F( t ) } \over \eta } + \displaystyle{{\dot{F}( t ) } \over E}} \right)$$

and,

(11)$$\dot{F}( t ) = \displaystyle{{\partial F( t ) } \over {\partial t}} = {-}\displaystyle{{F_{\,peak}{\rm e}^{{-}1/2{( {( {t-t_{\,peak}} ) /( {d/6} ) } ) }^2}( {t-t_{\,peak}} ) } \over {{( d/6) }^2}}$$

At the lower boundary, the snow rests on the base subassembly, which is modeled as a linear, elastic spring with spring constant, k _eff.

(12)$$w( {0, \;t} ) = \displaystyle{{\sigma ( {0, \;t} ) A} \over {k_{eff}}}$$

The aluminum plates are assumed to be rigid, so the spring constant, k _eff, is defined entirely by load cell deformation. According to the data sheet (TE Connectivity, 2018), each load cell deforms 0.05 mm at 445 N. Thus, the spring constant, k _lc, for each is 8.9 × 10⁶ N m^–1. Since there are four load cells in parallel, the effective spring constant, k _eff, is 4k _lc.

(13)$$k_{eff} = 4k_{lc} = 3.56 \times 10^7\;{\rm N\ }{\rm m}^{{-}1}$$

This boundary value problem is solved via two numerical methods, finite difference and finite element, to verify their solutions.

3.1.2. Finite difference method

The finite difference method involves discretizing the time and spatial domains, approximating partial derivatives with finite difference operators, and solving for displacement, w (Langtangen and Linge, 2017). At each time step, displacements are solved for in the spatial domain. Using the calculated displacements, the strain, strain rate, particle velocity and particle acceleration are ascertained by using various finite difference operators. Matlab is employed to implement the finite difference method. The implementation is detailed in the Dryad repository in the ‘FD method’ directory (Verplanck, 2024).

3.1.3. Finite element method

Abaqus, a commercial software, is used to implement the finite element method. The modeling space is specified as 2D and planar. Even though it is a 1D problem, a 2D space is used because Abaqus does not contain a 1D space option. The snow column is created as a deformable, wire part. The part is meshed into 100 equal sized truss elements (T2D2) spanning the column height, the same spatial discretization as the finite difference method. A dynamic, explicit solution is generated using a time step of 1 × 10⁻⁹ s. This time step is chosen by iteratively decreasing it until a stable, consistent solution is generated. Details of the finite element implementation are found in the Dryad repository in the ‘FE method’ directory (Verplanck, 2024).

3.2. Laboratory experiments

Laboratory experiments were conducted in the Subzero Research Laboratory at Montana State University. The structural testing chamber has a substantial concrete floor that provides a near-ideal base for impact experiments due to its negligible deformation under dynamic loads in this study.

The experimental design was inspired by the compression test (Jamieson and Johnston, 1996; Greene and others, 2016), a test used by avalanche practitioners to assess the instability of a snowpack at a specific field location. The test involves isolating a 30 cm by 30 cm column of snow and loading the top of it with progressively harder hand taps. Specifically, ten taps hinging from the wrist, ten taps hinging from the elbow and ten taps hinging from the shoulder. The laboratory loading sequences, in this study, are analogous to the hand-tap loading but done in a more repeatable manner with a dropped mass from a prescribed height.

Snow columns were made by sifting snow (4.75 mm opening size) into a break-away mold. The snow was left to sinter for 24 h at −5°C. After sintering, the mold was removed. Then, the snow was cut horizontally to be 60 cm tall and cut vertically to create two columns. In some tests, colder temperature testing was executed by allowing the column to sinter an additional 24 h at −10°C followed by another 24 h at −15°C.

One column is dedicated to snow property measurements and the other for impact testing (Fig. 2a). The snow was characterized according to the American Avalanche Association's standard (Greene and others, 2016). The density measurements were made at six heights using a 100 cm³ rectangular density cutter. The snow columns ranged in density from 135 to 428 kg m⁻³. Grain forms consisted of decomposing and fragmented particles and rounded grains.

Figure 2. Depictions of the laboratory experiments. (a) The two laboratory snow columns. The snow column on the left was dedicated to snow property measurements, and the one on the right was for impact testing. Wireless acceleration sensors were embedded side-by-side at heights of 10, 30 and 50 cm. The plate-system sensor assembly (b) was used to measure force and acceleration at the top and base of the column. The impact loading was executed by dropping weights on the impact subassembly.

Furthermore, penetration resistance was measured with both a thin-blade penetrometer (Borstad and McClung, 2011) and SnowMicroPen^® (SMP) (Schneebeli and Johnson, 1998). The thin blade penetrometer is a hand-held instrument resembling a paint scraper that records peak force upon insertion. It is inserted horizontally into a vertically exposed wall with the wide blade dimension transverse. The SMP is a cone penetrometer driven by a motor that records force as it travels through the snow. On each test day, six thin blade and three SMP measurements were made in the characterization column. The thin blade measurements were made in even increments spanning the height of the column (every 7–9 cm). The SMP was driven from the top surface at locations which did not interfere with thin-blade or density measurements. The penetration resistance from the thin-blade penetrometer is denoted as R _TB and that of the SMP is denoted as R _SMP. After impact testing, a single SMP measurement was done on the test column. The SMP provided a precise measurement of column height because it was rigidly mounted to an aluminum test stand. The detailed characterization of each snow column can be found in Supplementary Information Table S1 (Verplanck, 2024) and heterogeneity analysis in Figure S3.

During impact testing, both force and acceleration were measured at the top and base of the test column. A plate-system sensor assembly (Fig. 2b) was responsible for acquiring these data and was comprised of two subassemblies: the impact subassembly and the base subassembly. Each subassembly had four load cells (TE Connectivity FC-2231-0000-0100-L) and a three-axis accelerometer (Analog Devices ADXL356). These sensors were mounted between two square pieces of aluminum (30 cm × 30 cm). In each subassembly, the plate above the load cells was 6.4 mm (1/4″) thick and the plate below the load cells is 12.7 mm (1/2″) thick. The accelerometers are intended to measure acceleration at the boundaries of the snow column. Thus, they were mounted to the plates which contact the snow. The accelerometer in the impact subassembly was mounted on the top surface of the 12.7 mm thick, lower plate. The accelerometer in the base subassembly was mounted to the underside of the 6.4 mm thick, upper plate. The subassemblies were wired to a single data acquisition system (National Instruments cDAQ-9188) logging at 30 kHz. The impact subassembly had a guide rod for the dropped masses.

Acceleration was measured within the column using wireless acceleration sensors (Lesser, 2023) sintered into the snow during construction. These sensors were configured to record at 10 kHz with a ±10 g range. Each height of 10, 30 and 50 cm has two sensors side by side, for a total of six locations within the snow column.

The loading was done by dropping acetal (Delrin^®) masses from various heights onto the impact plate. Two loading methods were implemented: short duration and long duration. The intent of the short-duration loads was to generate a sharp (i.e. short pulse length) stress pulse through the snow with peak magnitudes which span the measurement range of the load cells. The intent of the long-duration loads was to mimic data on hand-tap measurements commonly used by practitioners. Impact duration was increased by placing a foam cushion to act as a buffer between the dropped mass and aluminum plate. The mass and drop heights were increased for long-duration loads to reach peak forces similar to that of hand taps. The potential energy associated with a drop is calculated by mgh, where m is the mass, g is the acceleration due to gravity and h is the height from which the mass is dropped. Higher masses and higher drop heights are associated with more potential energy. The three levels of long-duration impacts have lower peak forces, but higher potential energies than their short-duration counterparts. The drop heights and load characteristics for these impacts are shown in Table 1 and visualized in Supplementary Material Figure S3.

Table 1. A comparison of the two loading methods

The mean ± 95% confidence interval is calculated from all the relevant drops that are in the analyzed dataset. The loading sequence is inspired by the compression test and in the long duration tests an additional 20 drops are done after the first 30 to investigate the influence of the prior loading. The potential energy is calculated from the mass, drop height and acceleration due to gravity.

The loading sequence was inspired by the hand tap loading of a compression test. That is, a sequence of ten low impacts, followed by ten medium impacts, followed by ten high impacts. In the long duration tests, these 30 impacts were followed by ten more medium drops followed by ten more low drops to investigate the influence of the initial 30 drops. A total of 1000 individual impacts were analyzed across 22 snow column configurations.

3.3. Determination of model parameters

Wave speed, c, is the column height divided by the elapsed time for the wave to travel from the top to the base. On test days when the total height change was on the order of the levelness of the columns (≤5 mm), the column height was assumed constant. Three out of the 22 tests resulted in a change in height >5 mm and a piecewise linear interpolation estimated the height for each impact on these test days (Supplementary Material Text S1 and Fig. S4). The elapsed time was calculated as the difference in time of signal onset from the top to the base. The onset was determined by an autoregressive approach that calculates the minimum Akaike information criterion (AIC) of the continuous wavelet transform of the signal (Kurz and others, 2005; Kalkan, 2016). Other methods of determining the signal onset are discussed in Supporting Information Figure S5.

The elastic modulus, E, was calculated from wave speed, c, and density, ρ, by rearranging Eqn (4).

(14)$$E = c^2\rho $$

After calculating elastic modulus, the viscosity, η, was found by iteratively running the finite-difference, viscoelastic model until the modeled peak force at z = 0 matched the measured peak force from the base plate, within 1 N. A linear interpolation root finding algorithm was used for this calculation (Chapra and Canale, 2006). Finally, the damping coefficient was found by dividing the elastic modulus, E, by the viscosity, η.

3.4. Regression approach

Multiple linear regressions were performed to find empirical relationships between predictors (i.e. measurable snow properties: density, penetration resistance and temperature) and model parameters (E, η). Separate regression coefficients were found for the short-duration loading and long-duration loading due to their distinct differences illustrated in Section 4.2. Penetration resistance as measured by the SMP, R _SMP, was considered separately from resistance as measured by the thin blade penetrometer, R _TB because of the similarity of these measurements. Grain type and hand hardness were not included as predictors because they lack numeric values associated with them. Grain size was not used as a parameter due to the low resolution of the measurement method (crystal card with 2 mm grid). Other snow properties, such as specific surface area, were not used because they were not made directly.

A stepwise approach was employed to determine the recommended regressions following the guidance of Weisberg (2014). In general, there are 2^p − 1 possible linear combinations of predictor variables (p), excluding interaction terms and the case where all predictor coefficients are zero. In this application, there were three predictor variables: density (ρ), penetration resistance (R _TB or R _SMP) and temperature (T); thus, 2³ − 1 = 7 regressions to compare for each model parameter. Various criteria for selecting a regression could be used. One such criterion is the corrected Aikake Information Criterion, AIC _c (Sugiura, 1978). This criterion was chosen for model selection because it corrects for small sample sizes and balances model complexity with goodness of fit. A lower AIC _c value implies a better model. After running the seven initial regressions, if the equation with the lowest AIC _c score contained multiple terms, then interactions were considered. Interactions were not considered earlier as to not violate the marginality principle (Weisberg, 2014). This process is exemplified in Supplementary Material Text S2 and Tables S2 and S3. Multicollinearity concerns are addressed in Text S2 and Figures S6 and S7.

3.5. Acceleration data processing

As a validation procedure, modeled acceleration data were compared with measured acceleration data. Since no acceleration data were used in the parameter determination, they represent an independent dataset. Three metrics were acquired from the acceleration recordings: the first minimum, first maximum and reverberation time. See Figure 3 for examples of measured acceleration using the plate system and a wireless accelerometer. The first minimum and maximum were used as metrics rather than absolute minimum and maximum because they were consistently present in the recorded data and were reliably compared with the modeled first minimum and first maximum. For example (Fig. 3), the first minimum and absolute minimum are not the same values at the top of the column, but at the other two locations shown they are.

Figure 3. Examples of measured acceleration from an individual impact. For each recording, the first minimum, first maximum and reverberation time are calculated to be used as validation metrics. Acceleration at the top and base are measured with the plate-system sensor assembly. Acceleration within the snow is measured with wireless acceleration sensors. Only data from one wireless sensor are shown in the figure for simplicity. The wireless sensors are not time synchronized with each other, nor with the plate-system sensor assembly.

In this study, the wave onset was found using same AIC-picker method as when calculating wave speed (Kurz and others, 2005; Kalkan, 2016). The first minimum acceleration is the first local minimum after the wave onset. Following the first minimum is the first local maximum. The local minima and maxima were found with a minimum separation of 5 ms to other local minima and maxima, respectively. The end of the wave was defined as the timestamp when the mean of the absolute value of the following 10 ms does not exceed 2.5 m s⁻². This value is twice the resolution of the wireless acceleration sensors (Lesser, 2023). A variable threshold based on three times the std dev. of the noise prior to wave onset was considered. However, a combination of low noise levels prior to the signal and long tails of subtle vibrations (e.g. Fig. 3, 30 cm above base and at the base) led to a fixed threshold based on the sensor resolution to be a more reliable method of determining the end of the wave.

This process was carried out for each acceleration recording. Each individual impact contained eight acceleration readings: one at the top, six within the snow and one at the base. One thousand individual impacts are executed in the laboratory, resulting in 8000 possible acceleration recordings. Of the 8000 possible recordings, 6600 are analyzed in the validation procedure. The gap in recordings is attributed to the use of only two wireless sensors, one at 10 cm and one at 50 cm, on the last two test days (the other sensors were in use for a different study) and issues such as poor wired connections, low batteries and wireless sensors not responding to remote control.

4. Results

4.1. An example of an individual impact

Force is measured at the top and base of the column (Fig. 4). In this example, notice the wave reaches the base while the top is still being loaded (location a in Fig. 4). Since the base plate subassembly rests on a stiff concrete floor, much of the incident wave is reflected resulting in an amplification of force (location b). After the impact ends, the positive force measurements at the base (location c) indicate a rebound due to the minimal deformation of base plate subassembly on the concrete floor. After this initial rebound, a few oscillations occur, and the dynamic force dissipates (location d). These data are used to determine wave speed and viscosity for each of the 1000 individual impacts as described in Section 3.3.

Figure 4. An example of the measured force data from an individual impact. The blue line is the recorded force at the top of the column, the red line is the recorded force at the base of the column. The difference in time of wave onset from the top to the base is used to calculate a wave speed and elastic modulus. The peak compressive forces are used to calculate viscosity. A reflection off the laboratory floor causes amplification of stress at the base.

4.2. Average model parameters: short duration vs long duration

The wave speed, elastic modulus, viscosity and damping coefficient vs SMP penetration resistance are displayed in Figure 5. SMP penetration resistance is chosen as the snow property to plot against, rather than density, since unique SMP measurements are made on all the test days, but some densities were repeated across test days. Relationships between penetration resistance and density are explored in Supplementary Material Figures S6 and S7. Higher wave speeds, c, are measured with both shorter duration impacts and increased penetration resistance (Fig. 5a). Thus, elastic moduli, E, show a similar trend since the elastic modulus depends on wave speed (Fig. 5b). The viscosities, η, are generally larger for longer duration impacts and higher penetration resistance (Fig. 5c). The difference in parameters is most pronounced in the damping coefficient, E/η, which is roughly an order of magnitude larger for the short duration impacts (Fig. 5d) which does not show a clear relationship with penetration resistance.

Figure 5. The ascertained model parameters: (a) wave speed, (b) elastic modulus, (c) viscosity and (d) damping coefficient plotted against SMP penetration resistance. In all plots, the marker represents the mean, and the error bars are one std dev. in each direction. The std dev. is the statistical spread across the repeated impacts. An investigation into measurement uncertainty is found in Supporting Information Text S4 and Figure S8. The parameters are grouped by duration of impact.

4.3. Regression coefficients

After performing the stepwise approach discussed in Section 3.4, the regression coefficients are found based on the lowest AIC_c and are shown in Table 2. The regressions are permutations of Eqn (15) where the type of penetration resistance, R (thin blade, R _TB, or SnowMicroPen, R _SMP) is specified in the ‘Penetration resistance type’ column. All considered regressions are included in the dryad repository in the ‘Regressions’ folder (Verplanck, 2024) because snow measurements are often sparse and do not contain all the measured predictors from these laboratory experiments. Thus, one could still estimate E and η without necessarily possessing the recommended measurements in Table 2. The only interaction included in Eqn (15) is density-penetration resistance because that is the only one to have resulted in a lowest AIC_c score. The unit for density, ρ, is kg m³, penetration resistance, R, is N, and temperature, T, is °C.

(15)$$parameter\approx a_0 + a_1\rho + a_2R + a_3T + a_4\rho R$$

Table 2. A stepwise approach based on the lowest AIC_c score is employed to determine the linear regression for each parameter

Separation regressions are run for each impact duration (long or short), resulting in four parameters: long and short duration elastic modulus, E _long and E _short, and long and short duration viscosity, η _long and η _short. The type of penetration resistance, as measured by the SMP, R _SMP, [N], or thin blade, R _TB, [N], is included in the ‘penetration resistance type’. All considered regressions can be found in the Dryad repository in the ‘Regressions’ folder (Verplanck, 2024).

Notice the parameters for short duration loads (E _short, η _short) only have one predictor term (R _SMP in both cases). The parameters for the long duration loads (E _long, η _long) have more predictor terms. This is due, in part, to having a larger sample size of long duration tests (n = 17) than short duration tests (n = 5), a factor which influences the AIC _c value. Although more samples would make the regressions more robust, the regressions are still of utility as demonstrated in the following validation section. Furthermore, the number of terms in these regressions is in accordance with multicollinearity concerns, a topic discussed in Supplementary Material Text S3 and Figures S6 and S7.

4.4. Validation: acceleration

Each of the 22 column configurations (Supplementary Material Table S1) has three drop heights, which leads to 66 unique model runs each for the elastic and viscoelastic models. The models are run using the average impact durations and magnitudes of the peak force from Table 1 and the constitutive parameters (E, η) are calculated using the regressions in Table 2. The same metrics are extracted from these simulations as are pulled from the measured data. That is, first minimum, first maximum and reverberation time at the five different heights of acceleration measurements: base, 10, 30, 50 cm and top (Fig. 2).

To exemplify the validation procedure, consider the same impact as shown in Figure 4. First, E _long and η _long are calculated using density, thin blade penetration resistance and temperature according to Table 2. These parameters along with the relevant impact duration and magnitude are used to run the model. The measured and modeled accelerations are shown in Figure 6. There are ten measured trials (drops 21–30 on 7 December 2022) with which to compare these modeled results. The mean and std dev. of the measured first min/max are shown at the same time as the modeled first min/max. One of the ten trials (drop 23 as used in the demonstration example) is also plotted for reference. Each peak measured value is aligned with the modeled peak to synchronize the time axes.

Figure 6. An example of modeled and measured acceleration. The three metrics used for validation are first minimum, first maximum and reverberation time. The error bars are calculated using the mean and std dev. across the ten repeated impacts of the same drop height and duration (long-duration impacts, drops 21 through 30 in Table 1). Both the finite difference and finite element methods are used to solve the governing equations. This plot was made using the finite difference method, although nearly identical results are acquired using the finite element method (Supplementary Information Fig. S9). The acceleration decreases with depth because the laboratory base exhibits little motion compared with snow. Stress wave attenuation also contributes to a decrease in acceleration in the measurements and the damped wave equation, but not in the wave equation. The wave equation predicts the column to be perpetually in motion after impact.

Ratios between the modeled and measured metrics are calculated to make a quantitative comparison. The average ratios for all the experiments are shown in Table 3. A value greater than 1 implies the modeled value is greater than the measured value, and a value between 0 and 1 implies the modeled value is less than the measured value.

Table 3. Validation table comparing acceleration metrics

Each value is the median ratio of measured to modeled across all relevant drops.

Generally, the modeled acceleration values are greater in magnitude than the measured acceleration values. On the other hand, the modeled durations are shorter than the measured durations. The wave equation has no mechanism of energy dissipation, so the snow column theoretically continues to vibrate forever with the elastic model. The damped wave equation (i.e. viscoelastic) modeled values are closer to measured values except by the first minimum metric, in the long duration case. Furthermore, the measured metrics and viscoelastic modeled metrics generally grew closer together further down the snow until the base measurements.

In addition to the quantitative metrics, validation of the viscoelastic model is supported by the shape of acceleration curves, particularly with respect to the lower boundary condition. The base plate subassembly is modeled as a stiff spring which causes the calculated stress wave to reflect. The shape of this modeled acceleration is more similar to measured acceleration than if the column were to be modeled as semi-infinite, with no lower boundary. Consider our example case (7 December 2022) viscoelastically modeled with two different lower boundary conditions: laboratory base and semi-infinite. These modeled acceleration curves halfway down the 60 cm tall column are plotted in Figure 7 along with the measured recordings.

Figure 7. The shape of the measured acceleration curve resembles that of the laboratory base boundary condition, especially compared with the semi-infinite column. Although the model overpredicts magnitude of measured acceleration, the similarity in shape indicates a reflection off the base. The measured data are from the left accelerometer, 30 cm height and drop 23.

The reflection off the base results in a distinct shape of acceleration. The first minimum value is lower in magnitude than the first maximum value because the reflected wave decelerates the incident downwards-moving wave. After the loading has finished the column rebounds leading to positive acceleration. Then some oscillations occur, and the column comes to rest. In the semi-infinite case, the downwards-moving wave passes through the column with no interference resulting in a shape that is similar to the derivative of a Gaussian, but the first minimum is larger than the first maximum because of the viscous element. If modeled elastically, the shape would be identical to the derivative of a Gaussian.

4.5. Measured and modeled stress in laboratory columns

In addition to acceleration, stress is modeled throughout the column and measured at the boundaries. Stress in the example case used throughout this paper is shown in Figure 8 (7 December 2022, drop 23). These measured data and model results are not compared as a validation procedure because the measured stress data are used to create the regression of model parameters. Thus, comparing modeled and measured stress values would be considered overfitting. Despite this fact, it is worth noting that the damped wave equation agrees better with measured values at the base of the column than the purely elastic wave equation.

Figure 8. A comparison of modeled and measured stress in the example case, drop 23 from 7 December 2022. Stress is measured at the top and base, but not within the column. Finite difference and finite element methods are employed to numerically solve the governing equations. This plot was made using the finite difference method, although nearly identical results are acquired using the finite element method (Supplementary Information Fig. S10). The reflection off the concrete floor causes an increase in stress with depth, a result that is obvious with the wave equation but still present with the damped wave equation and measured values. The damped wave equation agrees well with measurements, but since these measurements were used to determine model parameters, they are not used in the validation procedure.

In this example, the modeled stress magnitude increases (Fig. 8) with depth whereas the modeled acceleration decreases with depth (Fig. 6) which can be explained by the boundary conditions. The base of the column is modeled as a stiff spring resulting in minimal displacement, and thus minimal acceleration. The increase in stress with depth is a result of the reflected-incident wave interference. The top boundary is modeled as Gaussian applied force which is effectively zero after the loading event, leading to a free top boundary. A free boundary exhibits maximal displacement and thus maximal acceleration.

4.6. Verification: finite element vs finite difference

The validation procedure was executed to determine the degree to which the model agreed with experimental data; to determine if the model was implemented properly, a verification procedure was executed. In our study, model verification is done by comparing the finite difference and finite element solutions. The 66 permutations are quantitatively compared for both elastic and viscoelastic constitutive relationships. Both acceleration and stress are compared at the base, 10, 30, 50 cm and top by the first minimum and first maximum metrics. The mean difference in metrics is 0.4% and is attributed to numeric error. This small error and overlapping results in Supplementary Information Figures S9 and S10 support model verification.

4.7. Model application examples

4.7.1. The influence of impact duration on model results

Consider two impact curves with the same peak force, but different durations applied to a semi-infinite column of snow. If the wave propagation is modeled using a purely elastic model, then the peak remains constant in both cases. If using the viscoelastic model, then, not only do both waves decrease in magnitude as they travel down the column, but the short-duration impact's peak diminishes at a shallower depth than long-duration impact.

As a theoretical example, consider a semi-infinite snow column with a 30 cm × 30 cm cross-section, density of 300 kg m⁻³, temperature of −10°C, thin blade penetration resistance of 8 N and SMP penetration resistance of 2 N. This column is subject to two different impacts. The two impacts both have a peak force of −500 N, but one has a 1 ms duration, and the other has a 10 ms duration. The peak stress as the wave travels down the column in these cases is shown in Figure 9, generated with the finite difference method.

Figure 9. Comparison of the influence of load duration and modeling methods. When modeling the snow elastically, the peak stress does not change no matter the duration. When modeling the snow viscoelastically, the peak stress diminishes. The shorter duration impact diminishes at a shallower depth than the longer duration impact.

When using the elastic model, peak compressive stress does not change as the wave travels down the column, no matter the impact duration. Using the viscoelastic model, both waves attenuate, and the shorter-duration impact diminishes at a shallower depth. This example demonstrates the more realistic results when using the Maxwell-viscoelastic model compared with the purely elastic model.

4.7.2. Comparing theoretic base reflections

To demonstrate the utility of this model, suppose there are four columns of snow: three rest on different bases and one is semi-infinite. The three bases are made of granite (ρ = 2700 kg m⁻³, E = 3.2 × 10¹⁰ Pa), ice (ρ = 917 kg m⁻³, E = 1.05 × 10¹⁰ Pa) and glacial till (ρ = 1800 kg m⁻³, E = 1.0 × 10⁸ Pa). The granite values are from Karagianni and others (2017) and Smithson (1971); the ice values are from Petrovic (2003); the glacial till values are from Bowles (1996) and Schueler and Holland (2000). The snow column is 1 m tall and has a 30 cm × 30 cm cross-section. The base material has the same cross-section and extends infinitely downwards. The snow column has a density of 300 kg m⁻³, temperature of −10°C and average thin blade penetration resistance of 10 N. The column is subjected to a load similar to a tap hinging from the shoulder; using the terminology from this paper, that is, a long duration impact from a drop height of 20 cm. The long duration regression is run on these snow properties to estimate elastic modulus (2.0 × 10⁷ Pa) and viscosity (5.5 × 10⁴ Pa s).

The effect of different bases is illustrated in Figure 10. The peak compressive stress is plotted throughout the column in the four different cases. The finite difference method was used to generate this solution.

Figure 10. Theoretic peak compressive stress throughout homogeneous columns of snow on different bases. Granite has the highest elastic modulus and thus has the greatest reflection, followed by ice and glacial till.

In the case where the column base is more of the same snow, there is no reflection. So, the peak stress simply decreases in magnitude as the stress wave travels downward. At a depth of 1 m, the snow has attenuated 48% of the peak stress. The glacial till, ice and granite have larger elastic moduli than snow; thus, the compressive stress wave is reflected and there is positive interference between the incident and reflected waves. This effect is most pronounced near the base of the snow column. The granite acts as a near-perfect reflector and the peak stress near the base of the column in this case is 1.98× greater in magnitude than the case without a reflection. Ice's peak stress at the base is 1.94× greater than the semi-infinite column. Importantly, in the granite and ice cases, the peak stress at the base of the column is greater than the peak applied stress at the top. The glacial till is an intermediary case; here, the peak stress is 1.66× greater than the case without a reflection.

5. Discussion

5.1. Validation discussion

One reason the modeled and measured metrics are not always 1:1 is because of measurement error. Capacitive, MEMS (micro-electromechanical system) accelerometers were chosen because they have been used to measure shock waves in snow from explosives (Binger and Miller, 2016). The specific model in our study (Analog Devices ADXL 356) was selected because of its configurable measurement range ± 10/±40 g, operating temperature down to −40°C, analog output and relatively high bandwidth (2.4 kHz). Although, the capacitive, MEMS accelerometers may not capture the true peak acceleration (Measurement Specialties, Inc., 2017). Another measurement error leading to smaller acceleration measurements is that the z-axes of the accelerometers are not perfectly aligned with the snow column because of manual sensor placement. A third cause of lower measured acceleration values is the energy loss between the snow and accelerometer. Although the wireless sensor housing material is intended to match the mechanical impedance of snow (Lesser, 2023), the transition from snow to accelerometer is not perfect. Likewise, the sensors may continue to vibrate longer than the snow itself due to internal and external resonances which would increase measured reverberation time.

The load cells (TE Connectivity FC-2231-0000-0100-L) were chosen due to their operating temperature down to −40°C, analog output, relevant measurement range, piezoresistive strain gauge and ease of implementation into a sandwiched-plate system. Despite these strengths, the data logger recorded force at a fixed rate of 30 kHz and the true peak forces may be slightly undermeasured, leading to inaccuracies in the model parameters. This contribution is small in comparison to other measurement uncertainties, particularly wave speed (Supporting Information Text S4).

Besides measurement error, another source of validation discrepancy is the mathematical idealization of a physical system. In the model, the impact force is idealized as a perfect Gaussian. In the laboratory experiment, the dropped mass bounces off the top plate and strikes again. These subsequent strikes are not modeled with the Gaussian. The time between strikes is smaller in the short duration impact experiments which contribute to the larger discrepancy in reverberation times observed in these experiments.

The modeled lower boundary condition assumes the deformation of the base-plate subassembly is entirely attributed to elastic deformation of the load cells and the piece of aluminum above the load cells is perfectly rigid. In reality, the 6.4 mm (¼″) thick piece of aluminum is not perfectly rigid and its deformation is expected to be greatest at the center where the accelerometer is mounted due to the four-point support from the load cells. This could explain the larger measured than modeled acceleration values at the base of the column for the long duration loads. The inverse observation of the short duration loads may be explained by insufficient force for significant plate bending.

Lastly, the Maxwell-viscoelastic constitutive model is an idealization. It is an improvement over a purely elastic model, yet it is not perfect. It assumes linear behavior and treats the snow as a continuum. Microstructure is ignored, and bulk behavior is modeled. The pressure wave through the gaseous-pore space is blurred with the pressure wave through the ice lattice. Micro-scale cracking of grain bonds and viscous/plastic deformation of ice grains are modeled as viscous deformation of the continuum.

In summary, the damped and elastic wave equations typically model acceleration values which are not 1:1 with measured values. The model values typically overestimate measured magnitudes of acceleration and underestimate reverberation times. The elastic wave equation predicts the snow column to vibrate forever, whereas the damped wave equation does not. The overestimation of magnitude and underestimation in reverberation may be attributed to both limitations of the sensors and model idealizations. The error is reduced at greater depths and the magnitude of acceleration approaches 1:1 at a depth 50 cm (Table 3). According to McClung and Schaerer (2006), the average crown thickness is ‘about two-thirds of a meter’ from a sample of 200 dry slab avalanches. In our study, the agreement in measured and modeled acceleration magnitudes at a depth of 50 cm demonstrates the utility of the viscoelastic model at common depths of slab avalanches.

5.2. Model limitations

The focus of this study's experiments is on the snow that is on the denser, harder side of that observed in seasonal snowpacks. Greene and others (2016) define a hard slab as having density of ≥300 kg m⁻³ but notes that, informally, hard slabs generally have a hand hardness of one finger or greater. Although 15 of the 22 snow column tests had densities <300 kg m⁻³, they all had hand hardness of one finger or greater. So, by Greene and others’ (2016) definition, we cannot strictly say the dataset is mostly hard slabs but, rather, on the harder side of the spectrum observed in seasonal snow. Our study's dataset is focused on snow columns that exhibit minimal deformation (<5 mm) under these impacts. The model parameters should not be applied to snow with properties significantly different than those which are tested in this study.

Likewise, the parameters should not be applied to loading situations which are significantly different from those in this study. The loading methods are discretely binned into two categories: short-duration impacts and long-duration impacts. Dispersive waves are waves that travel at different speeds depending on their wavelength. Although the impacts in our study do not have a wavelength, since they are not periodic signals, the impact duration (i.e. pulse length) is analogous to a wavelength. Since wave speed was observed to depend on impact duration, as illustrated in Figure 5a, dispersive wave transmission may be present in this study. The impacts in our work are treated categorically rather than as a spectrum of possibilities. A more complete model would include the dispersive effects of arbitrary pulse length.

Damage is neglected in the model. Occasionally, trends appeared between model parameters and subsequent impacts, but these trends were inconsistent and conflicting. With little to no observable height change, damage is presumed to take place at the microstructural level. Although our study is focused on bulk behavior, SMP measurements were made before and after impacts but did not show a discernable change in penetration resistance.

5.3. Relevance for avalanche practitioners

One of the findings is that snow's response is dependent on the impact duration. Higher wave speeds were observed in the short-duration impacts and these short-duration impacts attenuated at shallower depths. For the same peak force, a shorter duration impact will have less momentum than a longer duration impact since the momentum is the area under the force-time curve. These higher momentum impacts will transmit deeper into the snowpack. It has recently been suggested that one could achieve a similar force dropping a ski pole to load the snow in a stability test as tapping with a gloved hand (Sedon, 2021). Although this is plausible in terms of peak applied force, the duration of impact would be longer with a gloved hand and hence, the stress would transmit to greater depths with the hand tap. Thus, the snow's response would not be equivalent by these two loading methods.

Another finding is that stiff surfaces reflect and amplify the stress wave. This effect was observed by executing laboratory experiments on a substantial, concrete floor. The amount of the incident stress wave that is reflected depends on the elastic modulus of the stiffer material. A material like granite or ice is theorized to reflect more of the stress wave than glacial till. The scope of work in this study is limited to one spatial dimension and does not include observations of layered slabs, weak layers, crusts or ice layers. Thus, more work needs to be done to quantify this phenomenon in different layers observed within a mountain snowpack. From a theoretical standpoint, if the initial crack size in a weak layer was determined from a stress-based criterion (Reiweger and others, 2015), then one would expect to see larger initial cracks if the weak layer happened to rest on a stiff surface such as an ice layer.

An initial goal of this study was to better understand the effect of repeatedly loading the snow, as is the case in a compression test. Any apparent trends in model parameters with subsequent loads were inconsistent and conflicting, thus, a damage parameter was not included in this work. Most of the tests (19 out 22) resulted in a change in height on the order of the levelness of the columns (≤5 mm), so, damage is presumed to take place at the microstructural level. A study focused on microstructural change is recommended as future work to better understand the effect of repeated loads.

At the time of this study, there were limited data on how hard people tap during a compression test, particularly with respect to impact duration. Thus, the long-duration impacts that were intended to mimic hand taps were determined by an informal analysis of the authors' hand taps. Currently, there is an article in review that quantifies the hand taps of 286 avalanche practitioners (Toft and others, 2023) and the long-duration hand taps used in our study appear to be of similar peak force and duration to that of an average tap. Due to this similarity, the parameters developed for long-duration impacts in this study could be used to model the compression test.

5.4. Comparison to previous work

Much of the previous avalanche-focused work in observing snow's response to dynamic loads has taken place in the field and observation methods have been at lower sampling rates. The cameras used by Schweizer and others (1995) and van Herwijnen and Birkeland (2014) were 25 and 240 fps, respectively. The experiments which measured the time series stress response recorded at rates of 22 Hz to 2 kHz (Schweizer and others, 1995; Schweizer and Camponovo, 2001), 8–160 Hz (Thumlert and others, 2012; Thumlert and Jamieson, 2014, 2015) and 100 Hz (Griesser and others, 2023). Executing experiments in a laboratory-controlled environment with sensors recording at 10 and 30 kHz allows for improved observation of snow's response to impact loading.

The modeling aspect of our work intends to develop the simplest, effective model for dynamic stress wave transmission through snow in one dimension. Hence, we first consider a one-parameter elastic model, followed by a two-parameter viscoelastic model. In both cases, snow is treated as a homogeneous continuum. Others have applied more complex models such as Biot's or Burger's. Although these models may capture more nuance of snow's deformation, a limitation of them is the larger number of parameters necessary for a solution.

One study (Capelli and others, 2016) measured wave speed through columns of snow by initiating a signal with pencil lead fracture and using acoustic emissions sensors to detect travel time. They observed wave speeds of 300–950 m s⁻¹. These values are greater than but overlapping with the short-duration impacts in our study (300–500 m s⁻¹) and faster than our study's long-duration impacts (100–300 m s⁻¹). A pencil lead fracture is a shorter duration pulse than either of the loading methods in our study and agrees with our observation of shorter pulse lengths traveling at higher speeds. This trend is also in agreement with applications of Biot's model to snow that show higher frequency waves traveling at higher speeds, particularly for the waves transmitting through the pore space which are thought to carry most of the energy (Johnson, 1982; Capelli and others, 2016). Furthermore, these studies showed attenuation increases with frequency, a conclusion similar to ours of the inverse relationship between damping coefficient and impact duration.

A recent study determined elastic modulus of snow over a similar density range found it to range from 10 to 340 MPa and noted that ‘the true (high frequency, small strain) elastic modulus of snow is fairly high and should be distinguished from studies which are likely effected by the viscoplasticity of ice’ (Gerling and others, 2017). The elastic moduli determined in our study ranged from 1.5 to 94 MPa, overlapping but lower than that by Gerling and others (2017). Furthermore, our ascertained elastic moduli showed load dependence and were lower for the long-duration impacts stemming from the lower observed wave speeds from these impacts.

The viscosity in our work is effectively a parameter to include damping. To compare our ascertained viscosities with viscosities associated with a Kelvin–Voigt or Burger's viscoelastic model would not constitute an equivalent comparison. The Maxwell model has been applied to creep studies of snow and termed ‘axial viscosity’ by Mellor (1974). It is determined in uniaxial, constant compressive stress experiments and calculated by dividing the constant applied stress by the observed strain rate. This axial viscosity has been estimated at 10⁸ to 10¹⁴ Pa s (Mellor, 1974) over a similar density range to our study. The disagreement between our paper's viscosities (10³ to 10⁵ Pa s) and these ‘axial viscosities’ is attributed to substantially different loading methods and inherent meaning of the model parameter. The ‘axial viscosity’ is a parameter to describe creep, and viscosity in our study is a parameter to describe damping.

Although the elastic moduli and viscosities found in this study are rooted in well-established mechanics, they should not be applied to loading situations and snow types significantly different than presented here. They should be considered as parameters for a practical model of stress wave transmission, rather than ‘true’ elastic modulus as defined by Gerling and others (2017) or a viscosity associated with snow creep.

6. Summary and conclusions

Impact forces of different durations are applied to homogeneous snow columns in the laboratory and the snow's response is measured and modeled. Elastic and Maxwell-viscoelastic model parameters are determined using the measured force data and snow density. Multiple linear regressions are used to correlate measured snow properties with ascertained model parameters. As a validation procedure, modeled and measured acceleration values are compared. As a verification procedure, finite difference and finite element solutions are compared. Two theoretical model applications are explored: a semi-infinite snow column subject to impacts of different duration and a snow column resting on different base materials. This work leads us to the following conclusions:

• • The shorter-duration impacts result in higher wave speeds. Because of the wave speed's dependence on impact duration, different elastic moduli are calculated for the two impact durations.

• • The shorter duration impacts attenuate at shallower depths. The Maxwell-viscoelastic model includes a damping term whereas the elastic model does not. The ascertained viscosities are also determined separately for the two impact durations.

• • The laboratory's substantial concrete floor led to the reflection and amplification of stress waves at the base of the column. This phenomenon was captured in the numerical models. The magnitude of stress amplification depends on the elastic modulus of the base material; a granite base is theorized to amplify the stress more than ice and glacial till, respectively.

• • The Maxwell-viscoelastic constitutive relationship provides a two-parameter model for snow's response impact and is an improvement over an elastic model because it includes a damping term. The parameters determined in this study can be used to model snow's response to impacts of similar magnitude and duration.