Elastic wave velocities in snow and firn

Kohnen (Reference Kohnen1972) developed the following empirical relation between density and fast P-wave velocity based on measurements by Kohnen and Bentley (Reference Kohnen and Bentley1973) near New Byrd Station in West Antarctica:

(1)$$\rho( z) = {\rho_{\rm ice}\over 1 + [ ( V_{P, {\rm ice}}-V_P( z) ) /2.25 \, {\rm km\, s^{-1}}] ^{1.22}},\; $$

where ρ(z) and V _P(z) are the density and P-wave velocity at depth z, and ρ _ice and V _P,ice are the density and P-wave velocity of ice.

Diez and others (Reference Diez2014) gave a similar empirical relation between density and S-wave velocity based on measurements at Kohnen Station, West Antarctica:

(2)$$\rho( z) = {\rho_{\rm ice}\over 1 + [ ( V_{S, {\rm ice}}-V_S( z) ) /0.95 \, {\rm km\, s^{-1}}] ^{1.17}},\; $$

where V _S,ice and V _S(z) are the S-wave velocity of ice and at depth z, respectively. Relationships between density and velocity such as those in Eqns (1) and (2) have the potential to provide improved monitoring of the surface mass balance of glaciers and ice sheets using airborne or satellite altimetry by allowing the measured surface elevation to be corrected for compaction of snow.

Since snow and firn are porous media, their elastic wave velocities may be calculated using Biot theory (Biot, Reference Biot1956a, Reference Biot1956b, Johnson, Reference Johnson1982b; Sidler, Reference Sidler2015). Denoting angular frequency by ω, fluid viscosity by η and fluid density by ρ_f, the crossover from low to high frequency behavior occurs at a frequency ω_c at which the viscous skin depth, $\sqrt {2\eta /\rho _{\rm f}\omega }$, is approximately equal to the pore radius, a, (Johnson, Reference Johnson, Burridge, Childress and Papanicolaou1982a):

(3)$$\omega_{\rm c} = {2\eta\over \rho_{\rm f} a^2}.$$

Capelli and others (Reference Capelli, Kapil, Reiweger, Or and Schweizer2016) list pore diameters in snow in the range 0.2–0.35 mm. For air with η = 1.7 × 10⁻⁵  Pa s, ρ _f = 1.3 kg m⁻³ (Sidler, Reference Sidler2015), Eqn (3) gives ω_c < 3 kHz for this range of pore diameters.

At high frequency, this theory predicts that there is a fast compressional wave corresponding to the solid and fluid (air) moving in phase, a slow compressional wave corresponding to the solid and fluid moving out of phase, and one shear wave (Johnson, Reference Johnson, Burridge, Childress and Papanicolaou1982a). At low frequency, the fast compressional and shear waves propagate, but the slow wave is described by the solution of a diffusion equation. These predictions are consistent with observations of P-wave velocities that increase with density and are faster than in air when sources in direct contact with the ice frame are used, but that are slower than in air when sources not in contact with the ice frame are used (Capelli and others, Reference Capelli, Kapil, Reiweger, Or and Schweizer2016).

In the high frequency limit of the theory, the velocities of the fast and slow compressional wave and shear wave are given in the Appendix following Johnson (Reference Johnson, Burridge, Childress and Papanicolaou1982a). To illustrate the nature of these predictions, the empirical expressions for the required properties of snow as a function of porosity used by Sidler (Reference Sidler2015) based on prior information available in the literature are employed. Denoting porosity by ϕ, density and bulk modulus of the solid (ice) by ρ_s and K _s, and bulk and shear moduli of the ice frame by K _b and μ_b, the empirical relations used by Sidler (Reference Sidler2015) are as follows:

(4)$$K_{\rm b} = K_{\rm s}( 1-\phi) ^{30.85/( 7.76-\phi) },\; $$

(5)$$\mu_{\rm b} = {3K_{\rm b}( 1-2\nu_{\rm b}) \over 2( 1 + \nu_{\rm b}) }.$$

The following empirical relation for Poisson's ratio ν_b of the frame was used by Sidler (Reference Sidler2015):

(6)$$\nu_{\rm b} = 0.38-0.36\phi.$$

Sidler (Reference Sidler2015) assumes a value K _s = 10 GPa. Equation (5) then gives μ _s = 2.61 GPa. An additional parameter required is a geometrical quantity τ > 1 that increases as the tortuosity of the pore space increases. Sidler (Reference Sidler2015) used the expression derived by Berryman (Reference Berryman1980a) for the case of isolated spherical solid particles in the fluid:

(7)$$\tau = {1\over 2} \left( 1 + {1\over \phi}\right) .$$

Figure 1(a) shows the predictions of Biot theory in the high frequency limit obtained using these expressions.

Fig. 1.

Fast (full curve) and slow (dashed curve) P-wave velocities as a function of density as predicted by Biot theory. The circles in (a) show the data of Smith (Reference Smith1965) (red) and Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974) (blue), while the green stars show the data of Oura (Reference Oura1952), Ishida (Reference Ishida1965), Marco and others (Reference Marco, Buser, Villemain, Touvier and Revol1998) and Gudra and Najwer (Reference Gudra and Najwer2011). Panel (b) shows a detailed view of the fast and slow P-wave velocities for density in the range 170–190 kg m⁻³. The dash-dotted and dotted curves in (b) show the stiff frame predictions (K _b,  μ _b ≫ K _f).

Figure 1(b) shows a detailed view of the fast and slow P-wave velocities for density in the range 170–190   kg m⁻³, where the curves appear to cross in the upper subplot. The dash-dotted and dotted curves in the lower figure show the predictions of the theory in the stiff frame limit (K _b,  μ _b ≫ K _f, where K _f is the bulk modulus of the fluid). The fast and slow P-velocities in this limit are given by Johnson (Reference Johnson, Burridge, Childress and Papanicolaou1982a) as:

(8)$$V_{P, {\rm fast}}^{\rm stiff} = \sqrt{{K_{\rm b} + 4\mu_{\rm b}/3\over ( 1-\phi) \rho_{\rm s} + ( 1-1/\tau) \phi\rho_{\rm f}}},\; $$

(9)$$V_{P, {\rm slow}}^{\rm stiff} = {V_{\rm f}\over \sqrt{\tau}}.$$

A comparison of the fast and slow compressional velocities in Figure 1(b) with these expressions shows that at low density, when the ice frame has low bulk and shear moduli, the slow compressional wave travels mainly on the ice frame, while the fast compressional travels mostly in the air phase. At higher density, the fast compressional wave travels mostly on the ice frame, while the slow compressional wave travels mostly in the air phase. Above a density of 190 kg m⁻³ the fast wave is accurately represented by the stiff frame limit given by Eqn (8) which is essentially independent of τ because of the low density of air in comparison with that of ice.

The most complete published datasets with measurements of density, fast P- and S-velocities are those of Smith (Reference Smith1965) and Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974). The red and blue points in Figure 1(a) show these data, which are seen to be in reasonable agreement with the velocity of the fast compressional wave predicted by Biot theory using the empirical relationships of Sidler (Reference Sidler2015). The green stars show the data of Oura (Reference Oura1952), Ishida (Reference Ishida1965), Marco and others (Reference Marco, Buser, Villemain, Touvier and Revol1998) and Gudra and Najwer (Reference Gudra and Najwer2011). The data of Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974), Oura (Reference Oura1952), Ishida (Reference Ishida1965), Marco and others (Reference Marco, Buser, Villemain, Touvier and Revol1998) and Gudra and Najwer (Reference Gudra and Najwer2011) were taken from Capelli (Reference Capelli2021). It is seen that as density increases the green stars fall below the velocity of the slow compressional predicted by Biot theory using the parameters of Sidler (Reference Sidler2015). This suggests that τ is greater than that predicted by Eqn (7), and that the tortuosity of the pore space in snow and firn is greater than that assumed in deriving this equation. This is consistent with the observations of Schwander and others (Reference Schwander, Stauffer and Sigg1988) who found high tortuosity factors for firn samples, as compared to a sphere pack, indicating the existence of dead-end channels and/or relatively narrow passages in the pore channels. The absence of measurements of the slow wave at high density probably results from the high tortuosity of the pore space at high density.

Figure 2 compares the predictions of the empirical relations given by Eqns (1) and (2) with the measurements of Smith (Reference Smith1965) and Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974). The dotted curves correspond to P- and S-velocities V _P,ice and V _S,ice of ice estimated using density of ice of 915 kg m⁻³ given by Kohnen (Reference Kohnen1972) and bulk modulus K _s = 8.90 GPa and shear modulus μ _s = 3.46 GPa of isotropic polycrystalline ice calculated from the single crystal elastic stiffnesses of Gammon and others (Reference Gammon, Kiefte, Clouter and Denner1983) listed in Table 1 using the Hill (Reference Hill1952) average of the Voigt (Reference Voigt1910) and Reuss (Reference Reuss1929) bounds. The elastic moduli of Gammon and others (Reference Gammon, Kiefte, Clouter and Denner1983) were chosen for this comparison since these are the most recent measurements of the single crystal elastic stiffness coefficients of ice in Table 1. Also shown by the dash-dotted curves in Figure 2 are the predictions of Biot theory using the empirical relations of Sidler (Reference Sidler2015) listed above.

Fig. 2.

Comparison between measurements of fast P-velocity (a) and S-velocity (b) in snow and firn measured by Smith (Reference Smith1965) (red points) and by Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974) (blue points) with fast P-velocity and S-velocity obtained using the empirical relations of Sidler (Reference Sidler2015) (dash-dotted curves). The dotted curves show the predictions of Eqns (1) and (2) with P- and S-velocities V _P,ice and V _S,ice of ice computed using a density of ice of 915 kg m⁻³ given by Kohnen (Reference Kohnen1972) and bulk modulus K = 8.90 GPa and shear modulus μ = 3.46 GPa of isotropic polycrystalline ice calculated from the single crystal elastic stiffnesses of Gammon and others (Reference Gammon, Kiefte, Clouter and Denner1983). The full curves show the predictions of Eqns (1) and (2) using a reduced shear modulus μ = 3.20 GPa for isotropic polycrystalline ice.

Table 1.

Single crystal elastic stiffnesses (GPa) of ice measured by Jona and Scherrer (Reference Jona and Scherrer1952), Green and Mackinnon (Reference Green and Mackinnon1956), Bass and others (Reference Bass, Rossberg and Ziegler1957), Brockamp and Querfurth (Reference Brockamp and Querfurth1964), Dantl (Reference Dantl1968), Bennett (Reference Bennett1968) and Gammon and others (Reference Gammon, Kiefte, Clouter and Denner1983) compiled by Gusmeroli and others (Reference Gusmeroli, Pettit, Kennedy and Ritz2012) together with computed Voigt (V), Reuss (R) and Hill (H) estimates K _V, K _R, K _H, μ _V, μ _R, μ _H, M _V, M _R, M _H, ν_V, ν_R and ν_H of the bulk K, shear μ and compressional M = K + 4μ/3 moduli and Poisson's ratio of an isotropic polycrystalline aggregate of ice crystals

It is seen in Figure 2(a) that Eqn (1) with bulk modulus of isotropic polycrystalline ice obtained from the single crystal elastic stiffnesses of Gammon and others (Reference Gammon, Kiefte, Clouter and Denner1983) agrees best with the measurements of Smith (Reference Smith1965), while the fast P-velocity obtained using the empirical relations of Sidler (Reference Sidler2015) agree best with the measurements of Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974). Figure 2(b) shows that Eqn (2) with shear modulus μ _s = 3.46 GPa of isotropic polycrystalline ice obtained from the single crystal elastic stiffnesses of Gammon and others (Reference Gammon, Kiefte, Clouter and Denner1983) overestimates the measured S-velocities, while Biot theory using the empirical relations of Sidler (Reference Sidler2015) underestimates the S-velocities at high density. The shear modulus at zero porosity calculated from Eqns (4) and (5) used by Sidler (Reference Sidler2015) is μ _s = 2.61 GPa, which is significantly less than any of the shear moduli for isotropic polycrystalline ice given in Table 1. The full curve in Figure 2(b) shows the prediction of Eqn (2) for a shear modulus μ _s = 3.20 GPa. At high density, this is seen to give a much better prediction than that calculated from the data of Gammon and others (Reference Gammon, Kiefte, Clouter and Denner1983) and is the value used in the rest of this paper.

Poisson's ratio ν may be calculated from V _P and V _S as follows:

(10)$$\nu = {V_P^2-2V_S^2\over 2( V_P^2-2V_S^2) }.$$

Figure 3 compares Poisson's ratio ν as a function of fast compressional velocity computed by combining Eqns (1) and (2) using bulk modulus K _s = 8.90 GPa and shear modulus μ _s = 3.46 GPa of isotropic polycrystalline ice calculated from the single crystal elastic stiffnesses of Gammon and others (Reference Gammon, Kiefte, Clouter and Denner1983) (dotted curve) and by using the empirical relation of Sidler (Reference Sidler2015) given by Eqn (6) (dash-dotted curve). The full curve shows the predictions of Eqns (1) and (2) using a reduced shear modulus μ = 3.20 GPa for isotropic polycrystalline ice. It is seen that Eqns (1) and (2) may lead to very small and even negative values of Poisson's ratio at low density, while ν > 0 for isotropic materials. The scatter in Poisson's ratio in Figure 3 is much greater than that in the velocities in Figure 2 as is expected from Eqn (10). As an example, Christensen (Reference Christensen1996) notes that an error of over $9 \percnt$ in Poisson's ratio can result from errors of $2 \percnt$ in both V _P and V _S. For this reason, Poisson's ratio calculated from the measurements of Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974) are not shown in Figure 3, since V _P and V _S for this dataset were digitized from the plots in Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974) so that the errors in V _P and V _S are presumably greater than in the original measurements.

Fig. 3.

Poisson's ratio calculated by combining Eqns (1) and (2) and from the empirical relation of Sidler given by Eqn (6) (dash-dotted curve) compared with that calculated from the measurements of Smith (Reference Smith1965) (red points). The dotted curves show the predictions of Eqns (1) and (2) with P- and S-velocities V _P,ice and V _S,ice of ice computed using a density of ice of 915 kg m⁻³ given by Kohnen (Reference Kohnen1972) and bulk modulus K = 8.90 GPa and shear modulus μ = 3.46 GPa of isotropic polycrystalline ice calculated from the single crystal elastic stiffnesses of Gammon and others (Reference Gammon, Kiefte, Clouter and Denner1983). The full curve shows the prediction of Eqn (2) using a reduced shear modulus μ = 3.20 GPa of isotropic polycrystalline ice.

Although the empirical fits of Sidler (Reference Sidler2015), Kohnen (Reference Kohnen1972) and Diez and others (Reference Diez2014) give a reasonable description of the data over the porosity range for which they were calibrated, they do not provide much insight into the underlying mechanisms for such relations. For example, it is not clear how to modify such relations to investigate the role of anisotropy, pore shape, etc. In contrast, a physics-based approach allows variations in anisotropy, pore shape, etc. to be investigated. In the following section, one possible physics-based approach is introduced and compared with the measurements of Smith (Reference Smith1965) and Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974) plotted in Figure 2.

Theoretical approach

To estimate the bulk and shear moduli of the porous ice frame, the differential effective medium scheme will be used. This has been reviewed, for example, by Norris (Reference Norris1985) and Zimmerman (Reference Zimmerman1990), and has been used widely to estimate the elastic properties of porous media with interconnected pore space such as porous sandstones (e.g. Zimmerman, Reference Zimmerman1990; Berge and others, Reference Berge, Berryman and Bonner1993; Berryman, Reference Berryman2004; Li and Zhang, Reference Li and Zhang2010).

Consider a two phase medium with bulk and shear moduli of phase i =1, 2 denoted by K _i and μ _i. In the differential effective medium scheme the effective elastic properties of the medium are modeled by adding inclusions of one phase into the effective medium in a stepwise fashion. The effective medium begins as phase 1, but is altered at each step as an additional increment of phase 2 is added. Assuming that the two phases and the orientation distribution of inclusions are isotropic leads to two coupled differential equations for the effective bulk and shear moduli K _eff and μ _eff (Berryman, Reference Berryman1992):

(11)$$( 1-y) {\rm d\over {\rm d}y}K_{\rm eff}( y) = ( K_2-K_{\rm eff}) P( y) ,\; $$

(12)$$( 1-y) {\rm d\over {\rm d}y}\mu_{\rm eff}( y) = ( \mu_2-\mu_{\rm eff}) Q( y) ,\; $$

with initial conditions K _eff(0) = K ₁, μ _eff(0) = μ₁. The quantity y is the volume fraction of phase 2. The quantities P and Q depend on the shape of the inclusions and, for ellipsoidal inclusions, have been given by Berryman (Reference Berryman1980b).

Application to dry snow and firn

The elastic properties of ice may be estimated using the Hill (Reference Hill1952) average of the Voigt (Reference Voigt1910) and Reuss (Reference Reuss1929) bounds for the bulk and shear moduli of isotropic polycrystalline ice. Ice that occurs naturally on Earth has a hexagonal crystallographic symmetry. The elastic stiffness tensor of a crystal with hexagonal symmetry is transversely isotropic with the axis of rotational symmetry (c-axis) perpendicular to the basal plane. Taking the x ₃ axis to lie along the axis of rotational symmetry, the non-vanishing elastic stiffness coefficients are C ₁₁ = C ₂₂, C ₃₃, C ₄₄ = C ₅₅, C ₆₆, C ₁₂ = C ₁₁ − 2C ₆₆ and C ₁₃ = C ₂₃ in conventional two-index notation (Voigt, Reference Voigt1910; Nye, Reference Nye1985). Table I lists measured single crystal elastic constants for ice compiled by Gusmeroli and others (Reference Gusmeroli, Pettit, Kennedy and Ritz2012), together with the computed Voigt (V), Reuss (R) and Hill (H) estimates K _V, K _R, K _H, μ_V, μ_R, μ_H, M _V, M _R, M _H, ν_V, ν_R and ν_H of the bulk K, shear μ and compressional M = K + 4μ/3 moduli and Poisson's ratio ν of an isotropic polycrystalline aggregate of ice crystals. For a given set of elastic stiffnesses, the bulk modulus is seen to be tightly bound by the Voigt and Reuss averages, but the shear modulus is less constrained due to the greater shear wave anisotropy of single crystals of ice. Figure 2 shows that at high density the bulk modulus K _s = 8.90 GPa of isotropic polycrystalline ice obtained from the single crystal elastic stiffnesses of Gammon and others (Reference Gammon, Kiefte, Clouter and Denner1983) listed in Table I using the Hill (Reference Hill1952) average of the Voigt (Reference Voigt1910) and Reuss (Reference Reuss1929) bounds gives a good description of the P-wave data. However, the shear modulus μ _s = 3.46 GPa calculated from the single crystal elastic stiffnesses of Gammon and others (Reference Gammon, Kiefte, Clouter and Denner1983) is too high and a reduced value μ _s = 3.20 GPa gives a better description of the shear wave data at high density shown in the lower subplot of Figure 2. For this reason, we use the values K _s = 8.90 GPa, μ _s = 3.20 GPa below.

Figure 4 shows the bulk and shear moduli of the ice frame in snow and firn as a function of porosity obtained using the differential effective medium scheme (Norris, Reference Norris1985; Zimmerman, Reference Zimmerman1990). The pores are modeled as oblate spheroids. An oblate spheroid is an ellipsoid obtained by rotating an ellipse about its minor axis, with equatorial dimension a greater than the polar dimension c. The curves in Figure 4 are plotted for pore aspect ratios of 0.1, 0.2, 0.4 and 1. Also shown are the measurements of Smith (Reference Smith1965) and Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974). Porosity is calculated from measured density as follows:

(13)$$\phi = {\rho_{\rm s}-\rho\over \rho_{\rm s}-\rho_{\rm f}}.$$

Here, ρ is the density of snow or firn, ρ _s is the density of ice and ρ _f is the density of the fluid in the pore space, assumed to be air. The density of ice is assumed to be ρ _s = 915 kg m⁻³ given by Kohnen (Reference Kohnen1972). It is seen that at low porosity the predicted bulk and shear moduli are consistent with the pores having a spherical shape (aspect ratio = 1), as observed also in glacial ice. As porosity increases, however, the measurements fall below the curve corresponding to an aspect ratio of 1. This may indicate that the pores have a lower aspect ratio at higher porosity. However, the elastic moduli only become zero in the limit ϕ → 1 in this scheme, while the elastic moduli are expected to vanish at a critical porosity ϕ _c < 1, estimated by Sidler (Reference Sidler2015) to be ϕ _c = 0.8 for snow. Also, pores in firn close and trap interstitial air in airtight bubbles as porosity reduces below a value of ~0.1 referred to as pore close-off (Stauffer and others, Reference Stauffer, Schwander and Oeschger1985; Schaller and others, Reference Schaller, Freitag and Eisen2017; Fourteau and others, Reference Fourteau2020). Above this value of porosity the connectivity of the pore space is important.

Fig. 4.

(a) Bulk and (b) shear moduli (full curves) of the ice frame in dry snow and firn obtained using the differential effective medium scheme (Norris, Reference Norris1985; Zimmerman, Reference Zimmerman1990) for pore aspect ratios of 0.1 (lower curve), 0.2, 0.4 and 1 (upper curve). The symbols show the measurements of Smith (Reference Smith1965) (red points) and Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974) (blue points).

Although widely used, the differential medium scheme has the disadvantage that the elastic moduli of a porous medium are predicted to vanish only in the limit ϕ → 1. For most porous media, however, the elastic moduli only become non-zero as porosity drops below a critical value ϕ _c, referred to as the critical porosity (Nur and others, Reference Nur, Mavko, Dvorkin and Galmudi1998). To account for this, Mukerji and others (Reference Mukerji, Berryman, Mavko and Berge1995) proposed a modification of the differential effective medium scheme in which the properties of the inclusions are taken as those of the porous medium at the critical porosity. The pore space within the inclusions is therefore interconnected, and the pores are not assumed to have any particular shape. Sidler (Reference Sidler2015) estimates ϕ _c as 0.8 for snow, while the porosity of fresh, uncompacted snow is typically of order 0.9–0.95 (National Snow and Ice Data Center, 2021). In the following, the value ϕ _c = 0.9 is used for the purpose of illustration. The approach is not limited to this value, however, and different values may be used as required.

Figure 5 shows the results of a modified differential effective medium scheme that takes into account the vanishing of the elastic moduli at a critical porosity ϕ _c = 0.9 and the interconnectivity of the pore space above the pore close-off porosity, assumed to be 0.1. For porosities <0.1 the original differential medium scheme is used in which discrete pores are added to the effective medium. Above a porosity of 0.1, however, phase 1 is taken to have bulk and shear moduli of the medium in Figure 4 at a porosity of 0.1 and phase 2 is taken to have the properties of the medium at critical porosity ϕ _c = 0.9. Comparison of the predicted bulk and shear moduli with those calculated from the P- and S-velocity measurements of Smith (Reference Smith1965) and Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974) plotted in figure 5 shows that the modified differential effective medium scheme gives a good description of the data. The elastic moduli are found to vanish for porosity greater than the critical porosity ϕ _c = 0.9 without the need to assume small pore aspect ratios at high porosity as in Figure 4.

Fig. 5.

(a) Bulk and (b) shear moduli (full curves) of the ice frame in dry snow and firn obtained using the modified differential effective medium scheme described in this paper assuming a critical porosity ϕ_c = 0.9 for inclusion aspect ratios of 0.1 (lower curve), 0.2, 0.4 and 1 (upper curve). The symbols show the measurements of Smith (Reference Smith1965) (red points) and Yamada and others (Reference Yamada, Hasemi, Izumi and Sata1974) (blue points).