the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Seismic attenuation in Antarctic firn

### Stefano Picotti

### José M. Carcione

### Mauro Pavan

We estimate the seismic attenuation of P and S waves in the polar firn and underlying ice by spectral analysis of diving, refracted, and reflected waves from active-source three-component seismic signals obtained in 2010 on the Whillans Ice Stream (WIS), a fast-flowing ice stream in West Antarctica. The resulting quality factors are then successfully modeled using a rock-physics theory of wave propagation that combines White's mesoscopic attenuation theory of interlayer flow with that of Biot/squirt flow. The first theory describes an equivalent viscoelastic medium consisting of a stack of two alternating thin porous layers, both of which have thicknesses that are much greater than the pore size but smaller than the wavelength. On the other hand, in the so-called Biot/squirt-flow model, there are two loss mechanisms, namely the global Biot flow and the local flow from fluid-filled microcracks (or grain contacts) to the pore space and back, where the former is dominant over the latter. The fluid saturating the pores is assumed to be fluidized snow, defined as a mixture of snow particles and air, such as powder, with a rigidity modulus of zero.

- Article
(7965 KB) - Full-text XML
- BibTeX
- EndNote

The shallow parts of polar ice sheets and ice streams, commonly known as firn, are subjected to snow densification, a metamorphic process activated by the pressure gradient due to the accumulation of snow and by the temperature gradient at the near surface (Wilkinson, 1988). As a result, the density of firn increases continuously from the surface to the pore close-off depth, where the material has a density of about 840 kg m^{−3} and can be considered glacial ice (Herron and Langway, 1980). Below this depth, the compression of englacial bubbles further increases the ice density until the maximum value is reached. In the cold and dry internal parts of the Antarctic continent, the densification is slow and the firn layer thickness exceeds 100 m (van den Broeke, 2008). On the other hand, in the high-strain shear margins of ice streams, this metamorphic process can be enhanced due to the effect of strain softening (Oraschewski and Grinsted, 2022), leading to exceptionally thin firn layers.

Because of its density structure, polar firn is generally laterally homogeneous, with the seismic velocity increasing nonlinearly with depth. This velocity profile results in the continuous refraction of seismic waves (Greenhalgh and King, 1981), which are also called diving waves. There are many examples in the literature that use diving waves to estimate the firn velocity–depth structure by picking and inverting first-break traveltimes (e.g., Kirchner and Bentley, 1979; King and Jarvis, 2007; Picotti et al., 2015). The ice-fabric characteristics as a function of depth have been obtained by exploiting the P- and S-wave anisotropic velocities inferred from active-seismic surveys conducted in different settings (e.g., Blankenship and Bentley, 1987; Picotti et al., 2015). However, while tomographic methods for estimating seismic velocity in the whole ice column are well established, algorithms for quantifying the depth dependence of attenuation are underdeveloped. In particular, as far as we know, there are no examples of the vertical profiling and modeling of the intrinsic seismic attenuation of P and S waves in the polar firn so far.

Intrinsic loss is often quantified using the inverse quality factor $\mathrm{1}/Q$, which represents the fraction of wave energy lost to heat in each wave period (e.g., Carcione, 2022; Gurevich and Carcione, 2022). P-wave quality factors (*Q*_{P}) in ice have been measured by several authors in various depth ranges. They were found to range from as low as 5 in the temperate ice at the surfaces of mountain glaciers (Gusmeroli et al., 2010) up to 700 within cold polar ice caps (Bentley and Kohnen, 1976). This wide range of values indicates a strong dependence of the quality factor on temperature, as demonstrated by laboratory experiments (Kuroiwa, 1964). This dependence was also verified by Peters et al. (2012) from seismic measurements in Greenland, where *Q*_{P} decreases with depth due to an increase in temperature. In this case, *Q*_{P} was measured within narrow depth ranges using strong basal and englacial reflections. Furthermore, it is common practice to measure the average *Q*_{P} over the entire ice column using the primary and multiple reflection events (e.g., Holland and Anandakrishnan, 2009; Booth et al., 2012). Regarding the S-wave quality factor (*Q*_{S}), Clee et al. (1969) and Carcione et al. (2021) measured *Q*_{S} in warm mountain glacier ice and reported values of about 23 and 12, respectively. To our knowledge, no other measurements of *Q*_{S} in glacial ice or firn have been published in the literature.

While cold polar ice attenuation is low, this is not the case for polar firn. Bentley and Kohnen (1976) obtained a value of the ice *Q*_{P} from seismic refraction measurements at Byrd Station of about 715 at 136 Hz between 100 and 500 m depth and at an average temperature of about −28 ^{∘}C. Because this value is very high, the intrinsic attenuation in polar ice (below the firn) caused only a minor reduction in the signal amplitude in their case. On the other hand, Albert (1998) performed numerical simulations of seismic waves in firn at the South Pole, Antarctica, between 0 and 300 m, where he neglected seismic attenuation on the basis of Bentley and Kohnen's results. However, as we show in the present work, attenuation in firn cannot be neglected, since the *Q*_{P} factor has values between 4 and 50 down to 40 m depth, thus inducing strong energy loss and velocity dispersion near the surface.

Regarding the attenuation mechanisms, Bentley and Kohnen (1976) assume that their results are “consistent with damping in slightly contaminated ice by a combination of two mechanisms: molecular relaxation at temperatures colder than about −20 ^{∘}C and grain boundary viscosity at warmer temperatures”. However, their study does not provide a physical predictive microstructural model for attenuation. Grain boundary relaxation and melting theories based on the Arrhenius equation can be important at temperatures close to the melting point of ice. In the case of rocks, such a model has been applied to describe seismic attenuation at the crust and mantle (Carcione et al., 2018).

Here we consider three attenuation mechanisms based on physical models described in Pride et al. (2004), Carcione and Picotti (2006), and Carcione (2022). We show that wave-induced fluid flow generates enough heat to explain the levels of intrinsic attenuation. This flow occurs at different spatial scales that can broadly be categorized as “macroscopic”, “mesoscopic”, and “microscopic”. The macroscopic flow is the wavelength-scale equilibration that occurs between the peaks and troughs of a P wave. This mechanism was first treated by Biot (1956) and is simply referred to as “Biot loss”. The microscopic flow occurs when a wave squeezes a porous medium that has grain contacts and cracks; these respond with a greater fluid pressure than the main pore space, resulting in a flow from cracks to pores named “squirt flow”. Finally, mesoscopic length scales are those much larger than grain sizes but smaller than wavelengths. Heterogeneity across these scales may be due to frame variations or to patches of different immiscible fluids. In our case, because of the seasonal alternation, the firn can be considered a finely layered medium created through the deposition of two porous media, with one layer being snow-like with high porosity and the other ice-like with low porosity (e.g., van den Broeke, 2008). When a compressional wave squeezes such a medium, the effect is similar to squirt flow, with the more compliant portions of the material responding with a greater fluid pressure than the stiffer portions. What follows is a flow of fluid capable of generating a significant loss in the seismic band. In the case of Antarctic firn, the seismic wave attenuation cannot be explained by adopting a simple porous model consisting of a rigid structure of ice and a porous space filled with air. Such a simple ice–air model underestimates the seismic attenuation by orders of magnitude. Therefore, the fluid saturating the pore space in this case is assumed to be fluidized snow, a mixture of snow particles and air (Mellor, 1974; Maeno and Nishimura, 1979; Nishimura, 1996). In this study, we demonstrate that the replacement of air with pore-fluidized snow leads to increased attenuation and to quality factor values comparable with those obtained from seismic data. To our knowledge, this is the first attempt to use the concept of fluidized snow as a porous fluid in Biot's theory to model the wave attenuation in firn.

In the present work, we infer the P- and S-wave quality factor vertical profiles of the firn layer of Whillans Ice Stream (WIS – West Antarctica) using the three-component seismic data recorded on the Subglacial Lake Whillans (SLW) during the active-source experiment described in Horgan et al. (2012) and Picotti et al. (2015). First, we describe the data, the attenuation rock-physics theories, and the methodology used to extract the *Q* factor from the seismograms. Then, we analyze the spectral content of the diving-wave first breaks and refracted waves at the firn bottom, and we compute the *Q*_{P} and *Q*_{S} factors using the frequency-shift technique. Finally, we fit the two vertical *Q* profiles using the presented rock-physics theories. In the “Discussion”, we give a detailed explanation of the importance of this model in studies related to the physical properties of the firn and to the characterization of subglacial media using amplitude variations with offset (AVO) analysis (e.g., Peters et al., 2008; Booth et al., 2012). We describe an alternative procedure for calculating the average *Q*_{P} and *Q*_{S} of the ice column below the firn that uses the reflected waves at the bottom of the glacier and the *Q* profile of the firn. This procedure can be useful in cases where conventional methods of amplitude preconditioning for AVO analysis are not applicable.

In the austral summer of 2010–2011, a comprehensive surface geophysics program surveyed a location on the WIS. The main target of the survey (Horgan et al., 2012; Christianson et al., 2012) was the seismic and radar characterization of the SLW (Fig. 1), an active subglacial lake that is the subject of a subglacial access program Tulaczyk et al. (2014). Four active-source seismic lines were acquired during this program: one line parallel to the ice stream flow direction, following the long axis of the lake, and three transverse lines across the lake (Fig. 1).

Two types of seismic surveys were undertaken: one to obtain the characteristics of the firn and another for the imaging (Horgan et al., 2012) and to define the anisotropy of the whole ice column (Picotti et al., 2015). For the analysis of the seismic velocity fields within the firn, a shallow-refraction dataset was generated using stacked hammer blows onto a wooden body embedded in the snow. This experiment was located on the longitudinal line, as shown in Fig. 1. The data were recorded on a 48-channel seismic system using georods consisting of five 40 Hz geophones wired in series (Voigt et al., 2013). The digital sample interval was 0.25 ms, and the maximum offset was 470 m, which is enough to define the velocity profile of the whole firn layer. To acquire the three components of the ground motion, the georods were rotated through the three mutually perpendicular axes. The horizontal longitudinal and transverse components were parallel and orthogonal to the line direction, respectively. In the shallow-refraction survey, the trace spacing was increased with increasing offset to gain better resolution in the shallow part of the firn, which exhibits the strongest velocity gradient. The trace spacing was 1, 10, and 20 m for offsets lower than 10 m, those ranging between 10 and 290 m, and those larger than 290 m, respectively. For the purposes of this work, the generation and acquisition of S waves polarized in the orthogonal direction with respect to the seismic line (i.e., SH waves) is far better, because these are pure S waves, free of P-wave interferences. We also acquired the SV waves, which were polarized in the longitudinal vertical plane, but these signals are more affected by interferences from diving and multiply refracted P waves propagating in firn. The obtained seismograms are shown in Fig. 2, where the P and SH diving-wave first breaks picked for the traveltime inversion and for the spectral analysis are indicated.

All the other data (Fig. 1), which were used to perform the imaging and to define the seismic anisotropy of the entire ice column, were generated using 0.4 kg PETN (pentaerythritol tetranitrate) charges placed at a depth of 27 m by using a hot water drill. Data were recorded on two 48-channel seismic systems, and the sensors consisted of alternating single vertical 28 Hz geophones and georods spaced 20 m apart (Fig. 3a). In addition, multi-component data were acquired using three-component (3C) continuous recording stations (Fig. 3b). The 3C stations were first placed along the longitudinal line and spaced 24 m apart, after which they were moved to the transverse lines, where the distance was 24 or 240 m (Fig. 1). Each of these stations consisted of a Guralp 40-T broadband seismometer with a 40 s corner period that was coupled to a Reftek RT-130 data acquisition system equipped with GPS timing. The maximum offset was 4320 m for the data recorded using 3C stations and 1910 m for the other data. These data are described in more detail by Horgan et al. (2012) and Picotti et al. (2015).

## 3.1 Overview of the theory

The theory of wave propagation in firn combines White's mesoscopic attenuation theory of interlayer flow (White et al., 1975; Carcione and Picotti, 2006; Carcione, 2022) and that of Biot/squirt flow (Gurevich et al., 2010; Carcione and Gurevich, 2011; Carcione, 2022). The first theory describes an equivalent viscoelastic medium consisting of a stack of two thin alternating porous layers of thickness *d*_{1} and *d*_{2}, such that the period of the stratification is $d={d}_{\mathrm{1}}+{d}_{\mathrm{2}}$. The theory gives the complex and frequency-dependent stiffness *E*, equivalent to the P-wave modulus. The complex velocity is given in Appendix A. On the other hand, in the so-called Biot/squirt-flow model, there are two loss mechanisms, namely the Biot global-flow one (Biot, 1956) and local flow from fluid-filled microcracks (or grain contacts) to the pore space and back. Figure 4 shows a sketch of two sandstone grains in contact.

Gurevich et al. (2010) assumed that a compliant pore forms a disk-shaped gap between two grains and its edge opens into a toroidal stiff pore, where *R* is the radius of the contact (or crack) and *h* is its thickness. The model assumes that the material becomes stiffer when the fluid pressure does not have enough time to equilibrate between the stiff and compliant pores (grain contacts and main voids, respectively). The complex velocity is given in Appendix B.

## 3.2 Phase velocity and quality factor

Denoting the complex velocities corresponding to the mesoscopic (subscript 1) and Biot/squirt-flow (subscript 2) loss models as *v*_{1}, *v*_{2P}, and *v*_{2S}, the global P-wave phase velocity and dissipation factor are

and

respectively, where

and

In the above equations, Re and Im denote the real and imaginary parts. Since the dominant mechanism is the Biot/squirt-flow one, we assume as a first approximation that the P-wave phase velocity is that of this loss mechanism. Alternatively, the phase velocity (Eq. 1), *v*_{pP}(*ω*), can be obtained from *Q*_{P}(*ω*) by using the Kramers–Kronig relations as given in Eq. (8) of Carcione et al. (2020) (see also Eq. 2.141 of Carcione, 2022).

Equation (2) can be demonstrated if we consider that the decay factor of a plane wave along a distance *r* due to the effect of the two attenuation mechanisms is

where the *α*_{1} and *α*_{2P} terms are attenuation factors given by

For low-loss media, *v*_{p1} is the mesoscopic-loss phase velocity and *f* is the frequency (see Eq. 2.129 in Carcione, 2022). Substituting Eq. (6) into Eq. (5), we obtain

where we have assumed that *v*_{p1}≈*v*_{p2P}.

On the other hand, the properties of the S wave are solely described by the Biot/squirt-flow model, such that

and

## 3.3 Estimation of the quality factor

To estimate the quality factor from the seismograms, we adopted the frequency-shift method (e.g., Quan and Harris, 1997; Picotti and Carcione, 2006; Picotti et al., 2007). Let *S*(*f*) and *R*(*f*) be the quasi-Gaussian spectra observed at the source and at a receiver, respectively, located at a mutual distance *d* in a homogeneous and isotropic medium. The frequency-shift approach is based on the property that, as the wave propagates through the medium, the high-frequency part of the spectrum decreases faster than the low-frequency part (Quan and Harris, 1997). This effect may be quantified by measuring the resulting downshift $\mathrm{\Delta}f={f}_{S}-{f}_{R}$, where *f*_{S} and *f*_{R} are the frequency centroids of *S*(*f*) and *R*(*f*), respectively. Then, if we approximate the spectrum *S*(*f*) by a Gaussian with variance ${\mathit{\sigma}}_{S}^{\mathrm{2}}$, we have

where *v* is the P- or S-wave velocity and Δ*t* is the traveltime. Although this method is only applicable under the hypothesis of constant *Q* in the bandwidth of interest (Quan and Harris, 1997; Dasgupta and Clark, 1998), it can also be applied with good reliability to frequency-dependent *Q* media. For example, Picotti et al. (2007) estimated the frequency-dependent *Q* from numerically simulated seismograms by using the spectral-ratio and frequency-shift methods, which gave values in good agreement with the synthetic theoretical models. Moreover, they found that the classical spectral-ratio method is more sensitive to this problem, resulting in noticeable departures of the logarithm of the spectral ratio from the linear trend. On the other hand, the frequency-shift method is less sensitive to this problem when the spectrum of the signal is Gaussian because the quality factor variation is generally weak in the frequency band of active seismic experiments (typically 5–500 Hz). It becomes important only towards the edges of the seismic frequency band, where the integral contributions of the spectral amplitudes to the calculation of the variances and the frequency centroids are negligible (e.g., Picotti and Carcione, 2006).

In the more general case of inhomogeneous media, the frequency shift along a raypath can be expressed as

Quan and Harris (1997), where the integral is taken along the raypath and ${\mathit{\alpha}}_{\mathrm{0}}=\mathit{\pi}/Qv$ is the attenuation coefficient, i.e., the attenuation factor is linearly proportional to the frequency and defined as *α*_{0}*f*. To infer the vertical *Q* profile of firn, we consider the whole firn column as a sequence of homogeneous-*Q* layers. The *Q* of each individual layer is obtained by adopting a layer-stripping method that considers the cumulative attenuation of all the overlying layers and exploits the monotonic increases in the firn's P- and S-wave velocities. The layer-stripping method is a well-established technique in both traveltime and attenuation tomography (e.g., Yilmaz, 2001; Böhm et al., 2006; Rossi et al., 2007). The P- and S-wave velocity fields (Fig. 5a) were previously computed in Picotti et al. (2015) by the traveltime inversion of first arrivals, following the method of Herglotz and Wiechert (Herglotz, 1907; Wiechert, 1910; Nowack, 1990). This technique, applied for the first time to Antarctic seismic data by Kirchner and Bentley (1979), calculates the velocity–depth function from traveltime-offset-picked first breaks. It is well suited to situations where the velocity increases monotonically with depth, and for this reason it was successfully applied to firn by several authors (e.g., King and Jarvis, 2007). The resulting diving-wave paths in the velocity gradient (Fig. 5b) are obtained using an optimized ray-tracing algorithm based on the shooting method, as described in Picotti et al. (2015). The curves displayed in Fig. 5a represent improved versions of those originally presented in Picotti et al. (2015).

Let's consider layer *N* in this succession and assume that the value of *Q* in the overlying *N*−1 layers is already known. In our plane-layered model, the depth of each layer coincides with the depth of maximum penetration of each ray at increasing offsets. The ray corresponding to this layer reaches a receiver at the surface whose recorded first-break frequency centroid is *f*_{R}. Under our hypothesis, and in accordance with Eq. (11), the total frequency shift in the overlying *N*−1 layers is

where *v*_{i} and *d*_{i} are the velocities and the lengths of the ray segments in each layer, respectively, and the factor of 2 accounts for both the downgoing and the upgoing waves. Then, the quality factor of the layer *N* is obtained by Eq. (10), where *d*=*d*_{N}, *v*=*v*_{N}, and the frequency shift in the layer *N* is given by $\mathrm{\Delta}f=\mathrm{\Delta}{f}_{N}={f}_{S}-{f}_{R}-\mathrm{\Delta}{f}_{N-\mathrm{1}}$.

The layer-stripping method requires that the quality factor of the shallowest layer is known a priori. Moreover, the frequency centroid *f*_{S} and variance ${\mathit{\sigma}}_{S}^{\mathrm{2}}$ of the source must be known as well. These two values can be estimated using the first breaks recorded by the receivers placed close to the source. For this purpose, it is necessary that the acquisition geometry is properly designed to have a sufficient number of traces at very short offsets. Then, particular care must be devoted to acquiring first breaks with high signal quality, e.g., the wavelets must be free of interference effects and their amplitudes must not be clipped. If such signals are available, the following considerations allow for the computation of the required spectral properties. It is convenient in this case to estimate the difference in the frequency content of propagating wavelets at different receivers. We consider the closest available signal to the source as reference wavelet, which is compared with the first breaks recorded at increasing offsets. In other words, we compare the spectrum *R*_{1}(*f*) of the reference wavelet to that of another wavelet *R*_{2}(*f*), which propagates for an additional traveltime Δ*t* in the same shallow layer. If the traces are sufficiently close to each other and to the source, it is possible to characterize the quality factor of the uppermost layer as the average of all calculated values, and the variance ${\mathit{\sigma}}_{S}^{\mathrm{2}}$ can be approximated to that of the reference wavelet. Finally, the dominant frequency of the source *f*_{S} follows from Eq. (10) by using the reference first-break traveltime.

## 4.1 Quality factor vertical profiling

As previously described, a requirement for the layer-stripping technique is knowledge of the quality factor of the shallowest layer as well as the dominant frequency and variance of the source. We determined these properties by using the first breaks recorded close to the source: up to a maximum offset of 10 and 6 m, for P and S waves, respectively.

Figure 6 shows the time histories and the corresponding spectra of the considered P-wave first breaks.

Unfortunately, the first breaks recorded at offsets of less than 7 m are contaminated by the SV waves, and the amplitudes very close to the source are clipped. Taking the signal recorded at 7 m offset as the reference wavelet, we computed the frequency shifts and, using Eq. (10), we estimated the average P-wave quality factor (and standard deviation) of the shallowest layer to be *Q*_{P}=4.1 ± 1.3. Calculating the maximum penetration depth of the ray emerging at 10 m offset, we estimated the thickness of this layer to be about 2.7 m. Moreover, the estimated dominant frequency and variance of the source are *f*_{S}=574 Hz and *σ*_{S}=171 Hz, respectively.

Figure 7 shows the wavelets and the corresponding spectra of the considered S-wave first breaks.

Because the amplitudes of the first break recorded at 1 m offset are clipped, this trace cannot be used as the reference. Taking the signal recorded at 2 m offset as the reference wavelet, we computed the frequency shifts and, using Eq. (10), we estimated the average S-wave quality factor (and standard deviation) of the shallowest layer to be *Q*_{S}=1.9 ± 0.4. By calculating the maximum penetration depth of the emerging ray at an offset of 6 m, we estimated the thickness of this layer to be about 1.6 m. Furthermore, the estimated dominant frequency and variance of the source are *f*_{S}=496 and *σ*_{S}=146 Hz, respectively.

The frequency centroids of the diving P- and S-wave first breaks are obtained from the spectra of the selected waveforms, as described in Picotti and Carcione (2006). The amplitude integrals are calculated, for each wavelet, in the frequency band from zero to the maximum frequency of the signal. Since this high cutoff frequency depends on the signal-to-noise ratio, a statistic is calculated to evaluate the dispersion due to noise. The mean values of the obtained distributions of the frequency centroids for both the P- and S-wave first breaks are shown in Fig. 8. The corresponding standard deviations are less than 3 Hz. Both curves show a rapid decrease in the dominant frequency of the signal as a function of offset, except for an increase in frequency for the P waves between 10 and 30 m offset, which is an effect due to the diving waves traveling in that strong velocity–*Q* gradient.

As described above, the P- and S-wave velocity fields shown in Fig. 5a represent improved versions of those originally computed by Picotti et al. (2015) by the traveltime inversion of first arrivals, following the method of Herglotz and Wiechert. The uncertainties in the velocities were obtained by perturbing the first-break traveltimes according to the dominant frequency of the selected wavelets and then repeating the Herglotz–Wiechert inversion to produce different velocity profiles as a function of depth. The averages and standard deviations of the obtained velocity distributions are displayed in Fig. 5a. The maximum P- and S-wave velocities, verified using the first arrivals refracted at distant offsets on the seismograms acquired using the explosive source, are 3864 ± 30 and 1947 ± 25 m s^{−1} at 60 ± 5 m depth, respectively. At short offsets, errors increase due to the steep velocity gradient near the surface.

The two frequency curves shown in Fig. 8, the velocity profiles, the characteristics of the spectral source, and the *Q*_{P} and *Q*_{S} of the most superficial layer are the inputs for the layer-stripping procedure used to calculate the P- and S-wave quality factor profiles of the entire firn column. The uncertainties in the *Q*_{P} and *Q*_{S} values are obtained by repeating the procedure using the previously calculated frequency centroid and velocity distributions and perturbing the quality factors of the most superficial layers by their standard deviations. Then, from the two *Q*-factor distributions in the output, we derive the corresponding mean and standard deviation values with respect to depth.

Figures 9 and 10 display the resulting *Q*_{P} and *Q*_{S} profiles for the whole firn column, from the surface to the ice, together with the dominant frequency computed at the maximum penetration depth of each diving raypath.

The first plot (Fig. 9) shows a slow increase in *Q*_{P} from the previously computed minimum value of 4.1 ± 1.3 close to the surface to a value of 120 ± 35 at about 50 m depth. Then, it increases sharply to over 300 ± 140 at about 58 m depth. A noticeable increase in uncertainty below 40 m depth is due to a consistent reduction in both the traveltime and the frequency shift in each layer. The second plot (Fig. 10) shows a moderate increase in *Q*_{S} from the previously calculated minimum value of 1.9 ± 0.4 near the surface to an average maximum value of about 250 ± 90, which remains almost constant at depths greater than about 40 m.

Since the offset range is not sufficient to appreciate where the *Q*_{P} factor stabilizes at the maximum value in the ice, we use the dataset acquired with the explosive charges placed at 27 m depth (Fig. 3a) to analyze the refracted P-wave first breaks recorded at large distances up to the maximum available offset of 1910 m. Figure 11 shows an example of selected wavelets and the corresponding spectra together with a reference wavelet recorded at 670 m offset. All these signals travel along similar raypaths in the firn, while their paths in the underlying ice have different lengths. Following the same principle adopted for the computation of the *Q* factor in the shallowest layer, we compare the spectrum of a reference wavelet to those of other wavelets that propagate for additional traveltimes in the deep ice. Because these paths through the ice are long, the differences in spectral centroid can be better appreciated and the quality factors can be calculated more reliably. This procedure should be repeated for a large number of selected wavelets in order to perform a statistical analysis. In this case, the average P-wave quality factor of the ice is *Q*_{P}=380 ± 70 at 60 ± 5 m depth, where the average measured temperature is about −24 ^{∘}C (Engelhardt and Kamb, 1993). Using the same procedure, we also analyzed the refracted S-wave first breaks from 330 to 470 m offset (Figs. 2b and 8), resulting in an average ice S-wave quality factor of *Q*_{S}=260 ± 80 at 60 ± 5 m depth. This value confirms the result provided by the analysis of the diving waves.

## 4.2 Modeling of the seismic properties of the firn

Firn is assumed to be a deposition of two porous media, with one layer being snow-like with high porosity and the other ice-like with low porosity. The grains (ice) have the properties *K*_{s}=10 GPa, *μ*_{s}=5 GPa (shear modulus), and *ρ*_{s}=917 kg m^{−3} in both layers (Gurevich and Carcione, 2022). The fluid saturating the pores is assumed to be fluidized snow, which is defined as a mixture of snow particles and air, like powder, with a rigidity modulus of zero. We consider *K*_{f}=571 MPa, *ρ*_{f}=200 kg m^{−3} and *η*=0.1 Pa s. The properties of this material have been investigated by Mellor (1974), Maeno and Nishimura (1979), and Nishimura (1996).

For this study, it is important to measure the properties of the fluidized phase of the snow. This task can be performed during the coring of firn samples, as indicated in Nishimura (1991). The main apparatus consists of two parts: a fluidized snow feeder and an inclined chute, where it is possible to store the disintegrated snow in fluidized conditions. Measurements are made at a temperature of −15 ^{∘}C to avoid adhesion effects between snow particles. In this context, the bulk density, elastic velocity, and viscosity can be measured.

The physical properties of the firn layer are obtained from the density model using functions of porosity that have been shown to be suitable for snow. The density profile as a function of depth is obtained by using the following empirical relationship (Kohnen, 1972):

where *V*_{P}(*z*) is the vertical P-wave velocity displayed in Fig. 5a, and *V*_{P,ice}=3864 m s^{−1} is the velocity in ice, which we assume equal to the maximum computed P-wave velocity.

The porosity obtained from the experimental density (Eq. 13) is

Figure 12 shows the experimental density and porosity, where the former increases and the latter decreases monotonically with depth, mainly due to compaction.

Then, for each layer, the dry-rock bulk modulus which best fits the data of Johnson (1982) is

The dry-rock shear modulus is

where *ν* is the Poisson ratio. Below the pore close-off depth, approximately 35 m in this case, the medium is mainly ice (*ρ*≥840 kg m^{−3} and *ϕ*≤10 %), and the Poisson ratio is better approximated from the inverted wave velocities as follows:

(Mavko et al., 2009), where *V*_{P} and *V*_{S} are the P- and S-wave velocities displayed in Fig. 5a, respectively.

The permeability is

(Sidler, 2015; Gurevich and Carcione, 2022), where *C*=0.012 kg^{2} m^{−4} and we consider that the porosity *ϕ* is obtained from the experimental velocities and density as given by Eqs. (13) and (14).

For the mesoscopic loss model, we assume that *ϕ*_{1}(*z*)≈*ϕ*(*z*), ${p}_{\mathrm{1}}={p}_{\mathrm{2}}=\mathrm{0.5}$, and *ϕ*_{2}=*γ**ϕ*_{1}, where *γ* is small, i.e., layer 2 (ice) has a much lower porosity than layer 1 (snow). Here we assume that *γ*=0.1. On the other hand, the squirt-flow model has the following values of the parameters (Carcione and Gurevich, 2011): $h/R=\mathrm{0.015}$, *ϕ*_{c}=0.0002, and *K*_{h}=1.38*K*_{m}, where *K*_{m} is the bulk modulus with the grain contacts and cracks open. The attenuation is due to the softer layer with much higher porosity, since the cracks are open and the global-flow loss (Biot) mechanism is effective. The contribution of the stiff layer is negligible.

Figure 13 shows the P- and S-wave velocities and dissipation factors as a function of frequency close to the surface, where *ϕ*_{1}=0.7, *ϕ*_{2}=0.07, and the solid and fluid properties of the two layers are the same. We can see that the dominant mechanism is the Biot (global flow) one (Biot, 1956), with a strong relaxation peak (high attenuation) and velocity dispersion. The squirt-flow peak is weaker and located at lower frequencies, while the mesoscopic loss contributes only to the P wave. The P- and S-wave dissipation peaks occur at a frequency of approximately 282 Hz. On the other hand, the model predicts that below a given depth (from 40 to 50 m), the squirt-flow loss predominates over the Biot global one.

Let us now study the dependence of attenuation on depth. We consider a frequency of 282 Hz, corresponding to the dissipation peaks in Fig. 13. Figure 14a shows a comparison between the experimental and theoretical quality factor profiles as a function of depth, while Fig. 14b shows the experimental errors in *Q*_{P} and *Q*_{S}. As expected, the uncertainties increase with increasing quality factor, in agreement with the fact that the attenuation effect weakens with increasing *Q* and can be observed with high accuracy only over large distances. The differences between the theoretical quality factor and the profiled *Q* are within the experimental errors for almost the entire firn column.

Studying the structures of polar ice sheets and basal materials is essential for modeling the responses of ice masses to climate change. The mechanisms of the basal movement of glaciers, which strongly influence the mass balance of Antarctica, are poorly understood. The presence at the ice bottom of enormous quantities of water-saturated sediments, ponds, and subglacial lakes determines the dynamics of the ice streams (e.g., Bindschadler et al., 2003; Winberry et al., 2009), which represent the main drainage conduits of ice from the interior of Antarctica to the ocean. Furthermore, the identification and characterization of subglacial lakes is also important for studying the climate history of Antarctica, including the possible presence of extremophile life (Siegert et al., 2011).

Seismic reflection techniques are powerful tools for mapping the physical properties of subglacial environments; specifically, the roughness at the subglacial boundary, water occurrence and fluid layer thickness, rock and sediment type, sediment porosity and fluid saturation, and the thickness of the subglacial sedimentary layer (e.g., Blankenship et al., 1987). These parameters are all crucial for understanding the basal flow mechanism, which is modulated by the subglacial hydrology and deformation characteristics of the subglacial till. Seismic imaging is the only way to map these properties over large areas, as drilling to the ice bottom is extremely expensive and only provides local information.

The application of advanced AVO techniques has successfully exploited detailed tomographic velocity models and reflection amplitudes at the basal boundary to gather considerable information on the underlying material properties in various contexts (e.g., Blankenship et al., 1987; Anandakrishnan, 2003; Peters et al., 2008; Booth et al., 2012). However, quantitative estimates of these subglacial properties are difficult due to the many challenges that characterize AVO techniques. The two most important limitations are that the source characteristics are often unknown and that the attenuation of P and S waves can be highly variable vertically across the entire ice column (and in some cases also laterally) and poorly constrained. Both of these conditions reduce the ability to quantify bedrock reflection coefficients, affecting uncertainties about the physical properties of subglacial materials.

Our method allows for a better understanding of the seismic properties of the firn, which are useful for estimating the ice-sheet mass balance from satellite observations of ice-sheet elevation (e.g., Wingham, 2000; Alley et al., 2007). As explained, since *Q* is strongly influenced by the porosity (and density) profile of the propagation medium, measuring both velocities and quality factors has the potential to provide additional information on the physical structure of the firn, avoiding costly coring.

Furthermore, the present work provides alternative means to correct for seismic reflection amplitudes at the glacier floor and can be very useful in cases where conventional amplitude preconditioning methods for AVO analysis are impossible. Indeed, conventional methods ignore the complex *Q* structure of the firn and usually rely on the presence of multiple reflections to calculate the source amplitude and average quality factor of the entire ice column at normal incidence (e.g., Holland and Anandakrishnan, 2009; Booth et al., 2012). In addition to the fact that multiple reflections are not always present in seismic data, when they are available, they often show a poor signal-to-noise ratio (e.g., Dow et al., 2013). In the case considered in this work, there is no evidence of multiple reflections in the overall dataset. Therefore, the obtained firn velocity–*Q* profile can help in calculating both the source amplitude and the average quality factor of the ice mass below the firn.

Consider an explosive charge buried deep in the firn (in our case at a depth of about 27 m for most shots) and a ray propagating from the source to a surface receiver. Adopting our layered firn model, the source amplitude can be calculated considering the first breaks and Eq. (7), expressing the decay factor along a distance *r*_{i} in each homogeneous layer. Neglecting the transmission loss between adjacent quasi-layers, and denoting the attenuation factor in layer *i* by ${\mathit{\alpha}}_{i}=\mathit{\pi}f/{Q}_{i}{v}_{i}$, the total amplitude of the source *S*(*f*) is

where *A*_{c}(*f*) is the amplitude correction factor, *f* is the dominant frequency, *N* is the number of crossed layers, *R*(*f*) is the receiver amplitude, and *G* takes into account frequency-independent factors, for example the geometric spreading factor, which is proportional to the total length of the ray. To reduce the uncertainty, Eq. (19) can be applied to a large number of receivers for each shot in order to perform a statistical analysis.

Let us now consider the amplitude *R*(*f*) of a signal reflected from the bottom of the glacier. The reflection coefficient *R*_{c} is obtained by correcting *R*(*f*) for the decay factor as follows:

where the calculation of *A*_{c}(*f*) is now performed along the rays reflected from the bed. Thus, this calculation implies knowledge of the quality factor profile of the entire ice column. In this case, in the absence of englacial and multiple reflections, the average quality factor between the firn bottom and the bed is required. This value can be estimated, for both P and S waves, by following the same principle adopted for the calculation of the *Q* factor in the most superficial layer (direct waves at short offsets) and at the base of the firn (refracted waves at large offsets). This method does not require a priori knowledge of the source characteristics, only the firn's *Q* profile. We compare the spectrum of a reference wavelet with those of other wavelets propagating for further traveltimes in both firn and deep ice. For this purpose, we select a reference wavelet at a short or medium offset (ray A) and another wavelet at a large offset (ray B), as shown in Fig. 15. These signals travel along raypaths with slightly dissimilar lengths in the firn, while their paths in the underlying ice have very different lengths. The greater the difference between the lengths of these paths, the larger the difference in the centroids of signal spectra and the greater the reliability of the calculated quality factor.

Applying Eq. (11) to rays A and B, we obtain

where Δ*f*_{BA} is the difference between the frequency centroids of the reflected wavelets recorded at the receivers *R*_{B} and *R*_{A} (Fig. 15). Separating the contributions of ice and firn in the summations, the average quality factor of the ice mass below the firn is

As done for the refracted waves at the bottom of the firn, we adopted the dataset acquired with the explosive charges placed at a depth of 27 m (Fig. 3) in order to analyze the reflected waves from the glacier bottom. Figures 16 and 17 show two examples of selected wavelets and corresponding spectra up to maximum offset values of 1910 and 2892 m for P and S waves, respectively. These are the largest offsets available in our dataset that show signals with sufficient quality in terms of signal-to-noise ratio. We performed a ray tracing by adopting the algorithm described in Picotti et al. (2015), which allows the calculation of the ice velocity versus the propagation angle. Then, using the wavelet frequency centroids, the firn quality factor profiles, and the computed traveltimes along the ray segments, we obtained the average ice quality factors below the firn (for both P and S waves) using Eq. (22). This procedure should be repeated for a large number of selected wavelets in order to perform a statistical analysis. In this case, the mean P-wave quality factor of the ice (and standard deviation) is *Q*_{P}=320 ± 60. Using the same procedure, we also analyzed the reflected S waves (Fig. 3b), obtaining an average ice quality factor of *Q*_{S}=205 ± 50. These values agree with the result provided by the diving and refracted waves. However, this analysis shows a slight decrease in ice quality factors compared to the maximum values estimated at the bottom of the firn. This decrease reflects the substantial dependence of the quality factor on temperature (e.g., Kuroiwa, 1964; Peters et al., 2012), which increases with respect to depth in polar ice sheets and ice streams. In our case, at the SLW site, the temperature increases from a minimum value of about −25 ^{∘}C (annual mean temperature) near the surface to a maximum value at the bed which is close to its pressure melting point (Engelhardt and Kamb, 1993; Tulaczyk et al., 2014).

A physical explanation of the seismic attenuation (quality factor) in the polar snow and ice masses is essential to gaining insight into the ice sheet and deeper geological formations. This study shows a novel approach, based on a physical model which depends on parameters that can be derived from core sample analyses, to interpret and explain the observed seismic wave attenuation in firn and in the underlying ice. In particular, the proposed model is useful for performing data processing and amplitude variations with offset inversions (AVOs) to extract basal petrophysical properties from active-source multichannel seismic data.

In this work, we estimated the P- and S-wave attenuation profiles in the firn and underlying ice of Whillans Ice Stream (West Antarctica) from the spectral analysis of diving, refracted, and reflected waves of active-source three-component seismic data. The resulting experimental quality factors range from values lower than 5 at the surface to approximately 300 and 250 at about 60 m depth for P and S waves, respectively. Then, the P-wave quality factor further increases up to a maximum value of about 380 in the underlying ice. This attenuation model allowed us to infer the average P- and S-wave quality factors of the entire ice column beneath the firn up to the bed, which is otherwise not possible using standard methods. The estimated average *Q* factors are slightly lower than the maximum values for the ice at the firn bottom, i.e., approximately 320 and 205 for P and S waves, respectively. This decrease is in agreement with the increase in temperature as a function of depth typical of the polar ice caps.

The firn wave propagation model combines White's mesoscopic attenuation theory of interlayer flow with that of Biot/squirt flow, and the fluid saturating the pores is assumed to be fluidized snow. The dominant attenuation mechanism is found to be that of the global Biot flow, with the theory showing good agreement with the experimental values for the whole firn column. The model also predicts that below the pore close-off depth (approximately 35 m in this case), the squirt-flow loss predominates over the global Biot loss. The proposed model enables seismic experiments to be interconnected with mechanical laboratory analyses of firn samples collected on a study site, allowing for both the calibration of seismic surveys and the inversion of firn properties. To achieve this aim, knowledge of both the seismic velocity and the attenuation is also important because they could allow, at least theoretically, the porosity (and density) profile of the polar firn to be obtained by means of surface or borehole seismic experiments, which is useful for estimating the ice-sheet mass balance from satellite observations of ice-sheet elevation. In light of this information, surface seismic experiments or well logging in inexpensive hot-water-drilled holes could replace costly core drills.

The complex P-wave modulus of two periodic thin porous layers is (White et al., 1975; Carcione, 2022)

where

with ${p}_{l}={d}_{l}/({d}_{\mathrm{1}}+{d}_{\mathrm{2}})$, $l=\mathrm{1},\mathrm{2}$,

where *ω* is the angular frequency, and $\mathrm{i}=\sqrt{-\mathrm{1}}$.

For each layer,

is the dry-rock P-wave modulus, where *K*_{m} and *μ*_{m} are, respectively, the bulk and shear moduli,

is Gassmann's P-wave modulus,

is Gassmann's bulk modulus,

where *ϕ*, *K*_{s}, and *K*_{f} denote the porosity and the bulk moduli of the grains and saturant fluid, respectively. The coefficient *α* is known as the effective stress coefficient of the bulk material. Finally, *η* is the fluid viscosity and *κ* is the frame permeability.

Let *ρ*_{s} and *ρ*_{f} denote the mass densities of the grains and fluid, respectively, and let

denote the mass density of the bulk material.

The complex velocity is

This model considers only the P-wave attenuation, since the shear modulus *μ*_{m} is real.

The squirt-flow poroelasticity stiffnesses are the Gassmann bulk and shear moduli,

where

*ϕ* is the porosity, *K*_{m} and *μ*_{m} are the bulk and shear moduli of the drained matrix, and *K*_{s} and *K*_{f} are the solid and fluid bulk moduli, respectively. We explicitly indicate the functional dependence of *α* and *M* on *K*_{m}, since we shall replace this modulus by a modified matrix (or frame) complex modulus $\stackrel{\mathrm{\u203e}}{K}$ which includes the squirt-flow mechanism. In the same manner, *μ*_{m} will be replaced by $\stackrel{\mathrm{\u203e}}{\mathit{\mu}}$. The new moduli are complex valued and frequency dependent.

In the Biot/squirt-flow model, the bulk and shear moduli of the saturated rock at low frequencies are given by Gassmann's equations,

where $\stackrel{\mathrm{\u203e}}{K}$ and $\stackrel{\mathrm{\u203e}}{\mathit{\mu}}$ are the bulk and shear moduli of the modified frame, including the unrelaxation due to the presence of the squirt-flow mechanism, and *α* and *M* are given by Eq. (B2) with *K*_{m} substituted by $\stackrel{\mathrm{\u203e}}{K}$. For simplicity, we keep the same notation for the Gassmann moduli, but now they are complex valued and frequency dependent.

Gurevich et al. (2010) obtained the modified dry moduli in the following form:

where *K*_{m} and *μ*_{m} are the dry-rock bulk and shear moduli at the confining pressure *p*_{c}, *K*_{h} is the dry-rock bulk modulus at a confining pressure where all the compliant pores are closed, i.e., an hypothetical rock without the soft porosity, and *ϕ*_{c} is the compliant porosity. This is so small – nearly 0.001 for most rocks – that the total porosity *ϕ* can be assumed to be equal to the stiff porosity. The key quantity in Eqs. (B4) is the effective bulk modulus of the fluid saturating the soft pores:

where *J*_{0} and *J*_{1} are Bessel functions.

The complex velocities of the P and S waves are

and

respectively.

The data presented in this work are available through the following previous publications: Horgan et al. (2012) and Picotti et al. (2015).

SP mainly conceived the study, processed and analyzed all data, produced most of the figures, and wrote the manuscript. SP and MP designed the experiment and acquired the seismic data in Antarctica. JMC contributed to the basic ideas, proposing the theoretical models, producing some of the figures, and writing several parts of the text. All authors discussed the results and were involved in drafting the manuscript.

The contact author has declared that none of the authors has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

We would like to thank Huw J. Horgan and Sridhar Anandakrishnan for their support.

This work has been supported by the Italian National Program of Antarctic Research (PNRA – WISSLAKE Project) and by the National Science Foundation (grant no. 0944794, 0632198, 0424589).

This paper was edited by Adrian Flores Orozco and reviewed by Rolf Sidler and Matthias Steiner.

Albert, D. G.: Theoretical modeling of seismic noise propagation in firn at the South Pole, Antarctica, Geophys. Res. Lett., 25, 4257–4260, https://doi.org/10.1029/1998GL900155, 1998. a

Alley, R. B., Spencer, M. K., and Anandakrishnan, S.: Ice-sheet mass balance: assessment, attribution and prognosis, Ann. Glaciol., 46, 1–7, https://doi.org/10.3189/172756407782871738, 2007. a

Anandakrishnan, S.: Dilatant till layer near the onset of streaming flow of Ice Stream C, West Antarctica, determined by AVO (amplitude vs. offset) analysis, Ann. Glaciol., 36, 283–286, https://doi.org/10.3189/172756403781816329, 2003. a

Bentley, C. R. and Kohnen, H.: Seismic refraction measurements of internal friction in Antarctic ice, J. Geophys. Res., 81, 1519–1526, https://doi.org/10.1029/JB081i008p01519, 1976. a, b, c

Bindschadler, R. A., King, M. A., Alley, R. B., Anandakrishnan, S., and Padman, L.: Tidally controlled stick-slip discharge of a West Antarctic ice stream, Science, 301, 1087–1089, https://doi.org/10.1126/science.1087231, 2003. a

Biot, M. A.: Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low frequency range, J. Acoust. Soc. Am., 28, 168–178, https://doi.org/10.1121/1.1908239, 1956. a, b, c

Blankenship, D. D. and Bentley, C. R.: The crystalline fabric of polar ice sheets inferred from seismic anisotropy, in: The Physical Basis of Ice Sheet Modelling, Proceedings of the Vancouver Symposium, August 1987, Vancouver, British Columbia, Canada, IAHS Publ. 170, 17–28, ISBN 0947571361, 1987. a

Blankenship, D. D., Bentley, C. R., Rooney, S. T., and Alley, R. B.: Till beneath Ice Stream B. 1. Properties derived from seismic travel times, J. Geophys. Res., 92, 8903–8911, https://doi.org/10.1029/JB092iB09p08903, 1987. a, b

Böhm, G., Accaino, F., Rossi, G., and Tinivella, U.: Tomographic joint inversion of first arrivals in a real case from Saudi Arabia, Geophys. Prospect., 54, 721–730, https://doi.org/10.1111/j.1365-2478.2006.00563.x, 2006. a

Booth, A. D., Clark, R. A., Kulessa, B., Murray, T., Carter, J., Doyle, S., and Hubbard, A.: Thin-layer effects in glaciological seismic amplitude-versus-angle (AVA) analysis: implications for characterising a subglacial till unit, Russell Glacier, West Greenland, The Cryosphere, 6, 909–922, https://doi.org/10.5194/tc-6-909-2012, 2012. a, b, c, d

Carcione, J. M.: Wave Fields in Real Media. Theory and numerical simulation of wave propagation in anisotropic, anelastic, porous and electromagnetic media, 4th extended and revised edn., Elsevier, ISBN 9780323983433, 2022. a, b, c, d, e, f, g

Carcione, J. M. and Gurevich, B.: Differential form and numerical implementation of Biot's poroelasticity equations with squirt dissipation, Geophysics, 76, N55–N64, https://doi.org/10.1190/geo2010-0169.1, 2011. a, b

Carcione, J. M. and Picotti, S.: P-wave seismic attenuation by slow-wave diffusion: Effects of inhomogeneous rock properties, Geophysics, 71, O1-O8, https://doi.org/10.1190/1.2194512, 2006. a, b

Carcione, J. M., Poletto, F., Farina, B., and Bellezza, C.: 3D seismic modeling in geothermal reservoirs with a distribution of steam patch sizes, permeabilities and saturations, including ductility of the rock frame, Phys. Earth Planet. In., 279, 67–78, https://doi.org/10.1016/j.pepi.2018.03.004, 2018. a

Carcione, J. M., Gei, D., Santos, J. E., Fu, L.-Y., and Ba, J.: Canonical analytical solutions of wave-induced thermoelastic attenuation, Geophys. J. Int., 221, 835–842, https://doi.org/10.1093/gji/ggaa033, 2020. a

Carcione, J. M., Gei, D., Picotti, S., and Botelho, M. A. B.: On the instantaneous frequency and quality factor, Geophys. J. Int., 227, 735–745, https://doi.org/10.1093/gji/ggab250, 2021. a

Christianson, K., Jacobel, R. W., Horgan, H. J., Anandakrishnan, S., and Alley, R. B.: Subglacial Lake Whillans – Ice-penetrating radar and GPS observations of a shallow active reservoir beneath a West Antarctic ice stream, Earth Planet. Sc. Lett., 331–332, 237–245, https://doi.org/10.1016/j.epsl.2012.03.013, 2012. a

Clee, T. E., Savage, J. C., and Neave, K. G.: Internal friction in ice near its melting point, J. Geophys. Res., 74, 973–980, https://doi.org/10.1029/JB074i004p00973, 1969. a

Dasgupta, R. and Clark, R. A.: Estimation of *Q* from surface seismic reflection data, Geophysics, 63, 2120–2128, https://doi.org/10.1190/1.1444505, 1998. a

Dow, C., Hubbard, A., Booth, A., Doyle, S., Gusmeroli, A., and Kulessa, B.: Seismic evidence of mechanically weak sediments underlying Russell Glacier, West Greenland, Ann. Glaciol., 54, 135–141, https://doi.org/10.3189/2013AoG64A032, 2013. a

Engelhardt, H. and Kamb, B.: Vertical temperature profile of Ice Stream B, Antarct. J. US, 28, 63–66, 1993. a, b

Greenhalgh, S. A. and King, D. W.: Curved raypath interpretation of seismic refraction data, Geophys. Prospect., 29, 853–882, https://doi.org/10.1111/j.1365-2478.1981.tb01031.x, 1981. a

Gurevich, B. and Carcione, J. M.: Attenuation and dispersion of elastic waves in porous rocks: Mechanisms and models, SEG Publication, 26, https://doi.org/10.1190/1.9781560803911.fm, 2022. a, b, c

Gurevich, B., Makarynska, D., de Paula, O., and Pervukhina, M.: A simple model for squirt-flow dispersion and attenuation in fluid-saturated granular rocks, Geophysics, 75, N109–N120, https://doi.org/10.1190/1.3509782, 2010. a, b, c

Gusmeroli, A., Clark, R. A., Murray, T., Booth, A. D., Kulessa, B., and Barrett, B. E.: Seismic wave attenuation in the uppermost glacier ice of Storglaciären, Sweden, J. Glaciol., 56, 249–256, https://doi.org/10.3189/002214310791968485, 2010. a

Herglotz, G.: Über das Benndorfsche Problem der Fortpflanzungsgeschwindigkeit der Erdbebenstrahlen, Phys. Z., 8, 145–147, 1907. a

Herron, M. M. and Langway, C. C.: Firn densification, an empirical model, J. Glaciol., 25, 373–385, https://doi.org/10.3189/S0022143000015239, 1980. a

Holland, C. and Anandakrishnan, S.: Subglacial seismic reflection strategies when source amplitude and medium attenuation are poorly known, J. Glaciol., 55, 931–937, https://doi.org/10.3189/002214309790152528, 2009. a, b

Horgan, H. J., Anandakrishnan, S., Christianson, K., Jacobel, R. W., Alley, R. B., Heeszel, D. S., Picotti, S., and Jacob I. W.: Subglacial Lake Whillans – Seismic observations of a shallow active reservoir beneath a West Antarctic ice stream, Earth Planet. Sc. Lett., 331–332, 201–209, https://doi.org/10.1016/j.epsl.2012.02.023, 2012. a, b, c, d, e

Johnson, J. B.: On the application of Biot's theory to acoustic wave propagation in snow, Cold Reg. Sci. Technol., 6, 49–60, https://doi.org/10.1016/0165-232X(82)90044-1, 1982. a

King, E. C. and Jarvis, E. P.: Use of shear waves to measure Poisson's ratio in polar firn, J. Environ. Eng. Geoph., 12, 15–21, https://doi.org/10.2113/JEEG12.1.15, 2007. a, b

Kirchner, J. F. and Bentley, C. R.: Seismic short refraction studies on the Ross Ice Shelf, Antarctica, J. Glaciol., 24, 313–319, https://doi.org/10.3189/S0022143000014830, 1979. a, b

Kohnen, H.: On the relation between seismic velocities and density in firn and ice, Z. Geophys., 38, 925–935, 1972. a, b, c

Kuroiwa, D.: Internal friction of ice. III: The internal friction of natural glacier ice, Contributions from the Institute of Low Temperature Science, 18, 49–62, 1964. a, b

Maeno, N. and Nishimura, K.: Fluidization of snow, Cold Reg. Sci. Technol., 1, 109–120, https://doi.org/10.1016/0165-232X(79)90004-1, 1979. a, b

Mavko, G., Mukerji, T., and Dvorkin, J.: The rock physics handbook: tools for seismic analysis in porous media, Cambridge Univ. Press, https://doi.org/10.1017/CBO9780511626753, 2009. a

Mellor, M.: A review of basic snow mechanics, U. S. Army Cold Regions Research and Engineering Laboratory, 1974. a, b

Nishimura, K.: Studies on the fluidized snow dynamics, Doctor thesis, Hokkaido University, Contributions from the Institute of Low Temperature Science, Ser. A, no. 37, 1–55, 1991. a

Nishimura, K.: Viscosity of fluidized snow, Cold Reg. Sci. Technol., 24, 117–127, https://doi.org/10.1016/0165-232X(95)00023-5, 1996. a, b

Nowack, R. L.: Tomography and the Herglotz-Wiechert inverse formulation, Pure Appl. Geophys., 133, 305–315, https://doi.org/10.1007/BF00877165, 1990. a

Oraschewski, F. M. and Grinsted, A.: Modeling enhanced firn densification due to strain softening, The Cryosphere, 16, 2683–2700, https://doi.org/10.5194/tc-16-2683-2022, 2022. a

Peters, L. E., Anandakrishnan, S., Holland, C. W., Horgan, H. J., Blankenship, D. D., and Voigt, D. E.: Seismic detection of a subglacial lake near the South Pole, Antarctica, Geophys. Res. Lett., 35, L23501, https://doi.org/10.1029/2008GL035704, 2008. a, b

Peters, L. E., Anandakrishnan, S., Alley, R. B., and Voigt, D. E.: Seismic attenuation in glacial ice: A proxy for englacial temperature, J. Geophys. Res., 117, F02008, https://doi.org/10.1029/2011JF002201, 2012. a, b

Picotti, S. and Carcione, J. M.: Estimating seismic attenuation (*Q*) in the presence of random noise, J. Seismic Explor., 15, 165–181, 2006. a, b, c

Picotti, S., Carcione, J. M., Rubino, J. G., and Santos, J. E.: P-wave seismic attenuation by slow-wave diffusion: Numerical experiments in partially saturated rocks, Geophysics, 72, N11–N21, https://doi.org/10.1190/1.2740666, 2007. a, b

Picotti, S., Vuan, A., Carcione, J. M., Horgan, H. J., and Anandakrishnan, S.: Anisotropy and crystalline fabric of Whillans Ice Stream (West Antarctica) inferred from multicomponent seismic data, J. Geophys. Res.-Sol. Ea., 120, 4237–4262, https://doi.org/10.1002/2014JB011591, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

Pride, S. R., Berryman, J. G., and Harris, J. M.: Seismic attenuation due to wave-induced flow, J. Geophys. Res., 109, B01201, https://doi.org/10.1029/2003JB002639, 2004. a

Quan, Y. and Harris, J. M.: Seismic attenuation tomography using the frequency shift method, Geophysics, 62, 895–905, https://doi.org/10.1190/1.1444197, 1997. a, b, c, d

Rossi, G., Gei, D., Böhm, G., Madrussani, G., and Carcione, J. M.: Attenuation tomography: An application to gas-hydrate and free-gas detection, Geophys. Prospect., 55, 655–669, https://doi.org/10.1111/j.1365-2478.2007.00646.x, 2007. a

Sidler, R.: A porosity-based Biot model for acoustic waves in snow, J. Glaciol., 61, 789–797, https://doi.org/10.3189/2015JoG15J040, 2015. a

Siegert, M., Kennicutt, M., and Bindschadler, R.: Antarctic Subglacial Aquatic Environments, Geophys. Monogr. Ser., vol. 192, https://doi.org/10.1029/GM192, 2011. a

Tulaczyk, S., Mikucki, J., Siegfried, M., Priscu, J., Barcheck, C., Beem, L., and The Wissard Science Team: WISSARD at Subglacial Lake Whillans, West Antarctica: Scientific operations and initial observations, Ann. Glaciol., 55, 51–58, https://doi.org/10.3189/2014AoG65A009, 2014. a, b

van den Broeke, M.: Depth and Density of the Antarctic Firn Layer, Arct. Antarct. Alp. Res., 40, 432–438, https://doi.org/10.1657/1523-0430(07-021)[BROEKE]2.0.CO;2, 2008. a, b

Voigt, D. E., Peters, L. E., and Anandakrishnan, S.: “Georods”: The development of a four-element geophone for improved seismic imaging of glaciers and ice sheets, Ann. Glaciol., 54, 142–148, https://doi.org/10.3189/2013AoG64A432, 2013. a

White, J. E., Mikhaylova, N. G., and Lyakhovitskiy, F. M.: Low-frequency seismic waves in fluid saturated layered rocks, Physics of the Solid Earth, 11, 654–659, 1975. a, b

Wiechert, E.: Bestimmung des Weges von Erdbebenwellen, I. Theoretisches, Phys. Z., 11, 294–304, 1910. a

Wilkinson, D. S.: A pressure-sintering model for the densification of polar firn and glacier ice, J. Glaciol., 34, 40–45, https://doi.org/10.3189/S0022143000009047, 1988. a

Winberry, J. P., Anandakrishnan, S., Alley, R. B., Bindschadler, R. A., and King, M. A.: Basal mechanics of ice streams: Insights from the stick-slip motion of Whillans Ice Stream, West Antarctica, J. Geophys. Res., 114, F01016, https://doi.org/10.1029/2008JF001035, 2009. a

Wingham, D. J.: Small fluctuations in the density and thickness of a dry firn column, J. Glaciol., 46, 399–411, https://doi.org/10.3189/172756500781833089, 2000. a

Yilmaz, O.: Seismic data analysis: processing, inversion and interpretation of seismic data, SEG Series: Investigation in Geophysics, Tulsa, https://doi.org/10.1190/1.9781560801580, 2001. a