Near-surface seismic anisotropy in Antarctic glacial snow and ice revealed by high-frequency ambient noise

Abstract

Ambient seismic recordings taken at broad locations across Ross Ice Shelf and a dense array near West Antarctic Ice Sheet (WAIS) Divide, Antarctica, show pervasive temporally variable resonance peaks associated with trapped seismic waves in near-surface firn layers. These resonance peaks feature splitting on the horizontal components, here interpreted as frequency-dependent anisotropy in the firn and underlying ice due to several overlapping mechanisms driven by ice flow. Frequency peak splitting magnitudes and fast/slow axes were systematically estimated at single stations using a novel algorithm and compared with good agreement with active source anisotropy measurements at WAIS Divide determined via active sources recorded on a 1 km circular array. The approach was further applied to the broad Ross Ice Shelf (RIS) array, where anisotropy axes were directly compared with visible surface features and ice shelf flow lines. The near-surface firn, depicted by anisotropy above 30 Hz, was shown to exhibit a novel plastic stretching mechanism of anisotropy, whereby the fast direction in snow aligns with accelerating ice shelf flow.

1. Introduction

Flowing glacial ice is widespread in Antarctica, where snow accumulation is progressively being transformed by compaction processes into glacial ice, much of which eventually advects to the ocean and, in places, supports large ice shelves (e.g. Reeh, 2008; Hein and others, 2016; Rignot and others, 2019; Smith and others, 2020). Flow from the continental ice sheet to ice shelves in West Antarctica is governed by thinning and progressive acceleration through multiple outlet glaciers and ice streams, several of which, notably Pine Island Glacier (Joughin and others, 2021) and Thwaites Glacier (Parizek and others, 2013), are considered critical to the long-term stability of the West Antarctic Ice Sheet (WAIS) (Lhermitte and others, 2020). Weakening or collapse of ice shelves promotes and leads to increased inland ice flux to the ocean (Rignot and others, 2019) which, in turn, leads to a global rise in sea level. This effect was well documented following the collapse of the Larsen B ice shelf in increased flow rates from the grounded glaciers that had previously fed it (Rignot and others, 2004; Scambos and others, 2004). Many ice shelves in Antarctica are currently undergoing thinning (Rignot and others, 2013; Paolo and others, 2015), and this has prompted a number of initiatives aimed at improving our understanding of the dynamic response of ice shelves to atmospheric (Chaput and others, 2018), oceanic (Bromirski and others, 2010), geodynamic (Bassis and others, 2008; Scambos and others, 2009; Bromirski and others, 2015) and structural (Walker and others, 2013; Ledoux and others, 2017) forcing mechanisms. Some ice shelves have been particularly well studied in visual imagery and in situ structural interpretation. Studies of Ross Ice Shelf (RIS) for instance have included airborne radar, satellite image-based tracing of surface flow lines and other features (Fahnestock and others, 2000; Glasser and others, 2015), categorization of provinces (Ledoux and others, 2017) and observations of rift or crevasse seismicity (Olinger and others, 2019). These studies provide evidence for the persistence of large-scale crevasse and other anisotropic features. However, in cases where deeper investigations into ice structures are possible (e.g. radar-based studies), resolved sensitivities may reflect only a subset of effects inherent to buried structures. For instance, recent works (Young and others, 2021) at WAIS Divide have shown that polarimetric radar measurements are sensitive to ice crystal orientation, which can be strongly influenced by the historical strain regime of the sampled ice mass, but may not necessarily detect advected crevassing.

Fracture and ice fabric structures related to strain history may be modified by flow and evolution, producing accumulated and transported variations in bulk elastic structure. Spatial and temporal flow variability implies that the vertical seismic profile of an ice parcel may incorporate overlain fracture and/or strain features reflecting variations in snow and ice elastic moduli (e.g. Frolov and Fedyukin, 1998). These inherited features, when advected into a floating ice shelf, may significantly affect rift formation, flow and seismic wave propagation (Jeong and others, 2016; Lipovsky, 2018). Glacial ice in Antarctica furthermore often features thick firn cover (Diez and others, 2014; van den Broeke, 2008) in which surface snow with up to 90% porosity gradually densifies into solid ice over tens of meters. The small vertical scale, high parametric contrastand high sensitivity of the firn to atmospheric forcing on a broad range of time (e.g. minutes to years) and space (meters to 100 s of kms) scales make it both difficult to study broadly and locally and can lead to somewhat exotic wave propagation effects (e.g. Sidler, 2015) such as extremely slow compressional waves and intra-pore air waves. Where ice shelves are concerned, compaction and melting of porous firn layers leading to pond formation and subsequent hydrofracture are likely critical processes in destabilization and were integral to the rapid collapse of Larsen B in 2004 (Scambos and others, 2009; McGrath and others, 2012; Banwell and others, 2013; Kulessa and others, 2014; Alley and others, 2016; Kuipers Munneke and others, 2017). Remote means of monitoring (i.e. satellite imagery) were however unable to observe this process directly, motivating the development of other methods capable of both high-frequency sampling and adequate spatial and depth sensitivities.

Seismic waves offer a means to illuminate hidden ice shelf structures that complements the capabilities of overhead imagery. The ambient seismic wavefield at a seismograph site is generated by a superposition of processes that may tend to be statistically stationary when averaged over suitably long time periods. The relatively well-studied long-period microseismic spectrum between approximately 25 and 8 s period, for instance, is caused by coupling of ocean gravity waves with the solid Earth (Bromirski and Duennebier, 2002; Bromirski and others, 2005; Anthony and others, 2015; Nishida, 2017). The higher frequency ambient wavefield is generally a combination of wind coupling with Earth's surface (e.g. Withers and others, 1996) and with local structures. In Antarctica, this may include appreciable anthropogenic noise when camps or research stations are nearby. Ambient noise seismic interferometry has successfully been applied to study seasonal ice sheet and large glacial basal hydrological variations (e.g. Mordret and others, 2016; Zhan, 2019) using seismic waves with frequencies of up to a few Hz. Efforts to apply time-domain correlation-based seismic interferometry at higher frequencies are still evolving and have shown promise for temporal monitoring and imaging of near-surface dynamic structures (e.g. Picozzi and others, 2009; Larose and others, 2015; Mao and others, 2019), but typically suffer from instability in the high-frequency seismic source (i.e. wind), which in turn destroys coherence in repeated correlation functions. As shown by Chaput and others (2018) and Chaput and others (2022) however, high-frequency (greater than several Hz) seismic wave propagation in the very strong near-surface density and velocity gradient of the firn layer can lead to strong spectral amplifications that are pervasive throughout Antarctic firn seismic stations, persistent in time, and with spectral contents that respond acutely to multiscale environmental forcing. Leveraging these resonances directly in the frequency domain is thus an attractive avenue of research that not only permits single-station analyses, but appears largely insensitive to noise source instability.

Complex media with preferentially aligned structures tend to propagate seismic waves at different velocities depending on propagation direction and polarization. For transient (e.g. earthquake) sources, this can be observed in the time domain where incoming shear wavefields are split between orthogonally polarized S-waves, though Rayleigh and Love waves are also sensitive to anisotropy (e.g. Smith and Dahlen, 1973; Deschamps and others, 2008; Diez and others, 2014). Diez and others (2016), for example, documented anisotropic (approximately 4–18 Hz) surface wave dispersion on RIS using a seismic array and ambient surface wave excitation created by activity at a nearby seasonal camp. Observed radial S _V/S _H anisotropy inferred from Rayleigh and Love wave dispersion measurements was interpreted as being due to the layered character of the firn. Shear wave anisotropy is demonstrably sensitive to preferentially aligned structures such as crevasses (Aster, 2019; Lindner and others, 2019; Zhan, 2019; Hollmann and others, 2021), but may also be influenced by small-scale effects such as ice crystal axis alignment (Bennett, 1968; Picotti and others, 2015).

Here, we present a new approach capable of constraining azimuthal anisotropy in Antarctic firn and ice from passive seismic measurements. Ambient seismic measurements recorded at two seismic deployments in West Antarctica with very different apertures and recording durations (Fig. 1) were processed using an ambient noise methodology that estimates anisotropic parameters at individual stations in the frequency domain by observing ambient firn resonances (Chaput and others, 2018).

Fig. 1. Seismic stations used in this study. (a) The circular TIME array at WAIS Divide, consisting of 24 three-component Fairfield Gen2 nodal instruments, with the blue star marking the (central) location of calibration shots. The passive recording duration was approximately 250 h. The red ellipse in the map of Antarctica denotes Ross Ice Shelf. (b) Ross Ice Shelf Moderate Resolution Imaging Spectroradiometer (MODIS) imagery (Ledoux and others, 2017) showing the 34-station RIS/DRIS broadband seismographs deployed during 2014–2017 (Bromirski and others, 2015).

In the first case, we analyzed data from a circular 1 km aperture array near WAIS Divide Camp consisting of 24 nodal seismometers deployed as part of the Thwaites Interdisciplinary Margin Evolution (TIME) project (Barbuzano, 2020). At this site, active source data allowed for direct comparison of azimuthal anisotropy between shots recorded at the array and single-station ambient methods. Data from the ice shelf-wide RIS network consisting of 34 broadband stations deployed during 2014–2017 (Bromirski and others, 2015) were also reprocessed using the single station method, and results were interpreted in the context of satellite imagery, other ice shelf morphology and upstream ice-sheet processes.

1.1. The time-varying spectrum of firn-trapped seismic energy

Chaput and others (2018) noted that high-frequency ambient seismic excitation recorded by near-surface broadband seismometers in Antarctic firn provides stable information about near-surface elastic conditions. In that case, single station recordings of the ambient seismic wavefield on RIS revealed persistent spectral peaks that smoothly varied in frequency and amplitude over time scales of minutes to months for stations distributed across the ice shelf (e.g. Fig. 2). The peaks showed unique features at every station, but sometimes shared common spectral features and evolution at inter-station distances of up to many tens of km. Numerical simulations incorporating clustered surface sources over low velocity near-surface firn media with >80% porosity were shown to qualitatively reproduce some of the simpler spectral behaviors (i.e. those with only a few peaks and where frequency shifts could be accounted for by varying surface sources). Detailed physics of the source and medium interactions responsible for these signals remains somewhat speculative. Chaput and others (2018) hypothesized that these firn vibrations were excited by turbulent coupling between wind and semi-periodic dune and sastrugi features and/or firn-coupled seismic station infrastructure. Importantly, regionally correlated spectral features were found to be highly sensitive to near-surface changes associated with wind and temperature. Systematic spectral changes were also noted to occur during snow accumulation and scouring events during which surface snow topography was rearranged over a few hours. A prominent warming event that generated small amounts of near-surface melt over a large extent of RIS (i.e. 1000 s of km²) (Nicolas and others, 2017) also produced correlated spectral changes over a geographic region consistent with the melt region.

Fig. 2. (a) Moving-window power spectra for all three components across 200 h of continuous seismic velocity data at all three components of station 13 (shallow burial) of the TIME array, showing a strong firn resonance peak that drifts in frequency between approximately 20 and 30 Hz. This peak features prominent doublet behavior, though this is not always obviously solely from the spectrogram. (b) Similar plot for 20 days of continuous three-component data recorded at RIS station DR06, for a seismograph buried at ~2 m depth, showing a more complex pattern of evolving spectral peaks. The peak near 40 Hz shows clear doublet behavior, and most others do as well upon closer inspection. All three components show resonance peaks in all cases.

Recent analyses by Chaput and others (2022) of data from the 1 km aperture circular array of nodal seismometers deployed near WAIS Divide Camp (Fig. 1a) leveraged anthropogenic noise from WAIS Divide camp (5 km away) to reconstruct inter-station cross-correlation functions and study the behavior of up to 4 Rayleigh wave modes propagating across the array. These reconstructed Rayleigh waves were then analyzed in terms of motion ellipticity in the P-SV plane and horizontal to vertical energy ratios (H/V). Both those metrics were mapped to stations on the circular array, and it was found that ellipticity reversals in higher Rayleigh modes occurred in conjunction with lower peak frequency H/V ratios for those same modes. Such Rayleigh ellipticity reversals typically occur in sedimentary basins with severe parameter gradients (Denolle and others, 2012; Ma and others, 2016), and H/V ratios are often used in structural site response studies (e.g. Pina-Flores and others, 2016). Ambient spectral resonances as discussed here in this manuscript were also broadly observed on the TIME array (Fig. 2a), and their frequency content was found to vary spatially in a manner correlated with the other two observations. It was thus concluded that all three observations were sensitive to structural firn variations on the order of 100 s of m to a few km through surface wave excitation. As was the case for RIS data, the spectral peaks recorded on the TIME array exhibited splitting on the horizontals, but both peaks of the doublets generally showed a vertical expression as well (e.g. Fig. 2), further suggesting Rayleigh wave excitation. In this manuscript, the ‘doublet’ refers to the slightly offset pair of frequency peaks that compose what are referred to broadly as ‘spectral resonance peaks’ and ‘splitting’ refers to their offset, most often observed on the horizontal components.

Although environmental forcing such as storms often greatly perturbed the spectral peak patterns, the splitting of peak doublets was not affected by these forcings. When the frequency of doublets drifted, both peaks strongly tracked each other, suggesting separable source/structure effects.

The pervasiveness of these features leads to the natural query: Is spectral peak-doublet splitting associated with azimuthal anisotropy?

An azimuthal anisotropy ground truth is provided by active source determination at the TIME array, and we describe this processing below. Given limitations of the amplitude spectrum in terms of clearly discriminating peak doublets, we propose a practical approach to map these split peaks to both offset magnitudes and corresponding particle motion axes.

2. Methods

This section is organized as follows.

1. (1) We describe active source data processing at the TIME array, where multimode surface wave azimuthal anisotropy is directly resolved.

2. (2) We describe the requirements for an effective peak tracking algorithm to systematically analyze peak splitting, and propose a projection more suitable than time/frequency (spectrogram) for locating subtle peak doublets.

3. (3) A resonance peak splitting algorithm is proposed to systematically track peak doublet offsets at single stations and retrieve information relative to peak particle motions.

2.1. Active source anisotropy determinations at the TIME array

The TIME array data provide observations of passive spectral resonance peak splitting and well-constrained active sources at the center of the array. The analysis of these offset peaks in terms of firn and ice azimuthal anisotropy can thus be directly compared to ground-truth observations in the active source dataset.

As part of source testing for future experiments, several shot configurations with an accurate circular geometry (i.e. high fidelity 1 km diameter array, shots at the center, Fig. 1a) were performed, yielding a direct measurement of azimuthal anisotropy. The subsequent 10 days of passively recorded data then showed clear examples of spectral peaks (e.g. Fig. 2a), though unlike the 2 years of data recorded on RIS (Fig. 2b), the total frequency range of the peaks was more sparse (due to more limited spectral peak drift in time).

Given that we are solely interested in azimuthal anisotropy as seen by Rayleigh waves, we first analyzed the active source gathers and directly evaluated frequency-dependent delays in the surface wave component of the gather. Figure 3a shows a typical shot gather recorded on the TIME array near WAIS Divide. The much larger velocity of the direct P-wave makes calculating delays difficult due to sample rate limitations. However, variations in travel time are clearly observed in the surface wave component (Fig. 3b) and show a period of π (i.e. over one revolution of the array circle, there are two fast/slow cycles in Rayleigh arrival times), consistent with the presence of roughly orthogonal axes of azimuthal anisotropy.

Fig. 3. Depiction of multimode surface wave anisotropy calculations at the TIME array (Fig. 1). (a) 2 s shot gather recorded at the 24 stations of the circular array. The large arrival near 1.5 s showing a 2π directional dependence is the seismically coupled air wave, which displays wind-generated anisotropic variation. (b) Zoom-in of (a), showing an approximately π pattern in the surface wave arrival times (i.e. two full cycles over the array circle in Fig. 1). The red and black arrows show the slow and fast axes clockwise from South shown in (e). (c) Spectrogram of the vertical component of the shot arrivals at station 11, showing four distinct Rayleigh modes. (d) Isolation and tracking of independent modes.

Anisotropy estimated via surface waves is complicated by the wavefield containing multiple surface wave modes overlapping at similar frequencies, with each mode having its own depth sensitivity. If the medium has depth-varying anisotropy, then each mode may exhibit a different depth-dependent magnitude and direction. As such, instead of processing these data in the time domain (i.e. from Fig. 3b), which would involve poor discrimination, we perform a time/frequency operation on the recorded signal (Figs 3c,d) at each station (i.e. magnitude spectrogram), yielding excellent mode separation at frequencies above ~15 Hz. For each pair of stations (i.e. stations on opposite sides of the shot location), we choose one of the four visible modes in the spectrum (as picked in Fig. 3d at station 34) and stretch the arrivals at one station to align them with arrivals at the other station. We omit any inference from mode 3, however, as it is very weak and thus does not lend itself to reliable tracking. At every frequency in Figure 3d, the ‘fast’ arrival pattern arrives first and is thus time-stretched to the ‘slow’ arrival pattern using: $A_{{\rm slow}} = A_{{\rm fast}} + ( A_{{\rm fast}}\;\ast\; S_{\% })$, where A _slow and A _fast are the slow and fast arrivals and $S_{\% }$ is grid searched from 0 to 0.5. For each station pair, each frequency bin and each mode, the stretch value corresponding to the best alignment in arrivals is retained, thus resulting in a series of 24 × 24 matrix of stretch values. The maximum of those matrices allows for direct estimates of fast/slow axes, since the largest discrepancy (and thus largest stretch) in travel time will occur when stations on the fast and slow axes are compared.

2.2. Spectral peak analysis

As a first step in this analysis, we must accurately trace spectral peak doublet offsets through a variety of temporally variable manifestations. In implementing this method, several complications must be considered.

1. (1) Spectral peaks are observed to strongly change in amplitude and/or to disappear and reappear on time scales of hours to weeks.

2. (2) Individual component-offset peaks (i.e. ones displaying splitting, which we define as splitting observed between two linked peaks with approximately orthogonal particle motion orientations) commonly track each other in frequency on various time scales, precluding implementation of simple average spectrum calculations.

3. (3) Amplitudes between orthogonal components may trade off temporally, with one of the components presenting reduced signal to noise during some time periods.

4. (4) The degree of splitting is observed to be frequency dependent but no linearity in this behavior is assumed (i.e. peak offsets to not scale linearly as a function of frequency). Some peak pairs in one band may glide in frequency independently of other peak pairs.

5. (5) Spectral peak populations can be dense, and candidate paired anisotropic peaks may be difficult to discriminate in some situations.

As shown in Figure 2, a spectrogram representation, whereby the amplitude term of a Fourier transform of a short sliding time window data is taken for each component separately, suffers from amplitude trade-offs on the horizontal components and yields no information on particle motions. Point 4 further implies that a simple doublet frequency domain stretching operation applied to the power spectra (similar to what is used in ambient noise temporal monitoring in the time domain (Brenguierand others, 2008) to quantify coda phase changes) is ineffective in this situation, since standard coda doublet and stretching methods both assume that phase offset is a linear function of time (in our case, it would be a linear function of frequency). This is not the case here for our data as the relationship between anisotropy and depth is not assumed to be constant, both in magnitude and azimuth. Indeed, it can easily be observed that pairs of peaks often show frequency-dependent polarizations. We present a workflow designed to address these points.

2.3. Ambient noise covariance matrix decomposition

In seeking to address amplitude variations in peak doublets (points 1 and 3), we first propose to map our data to a displacement eigenvalue ratio rather than a spectrogram. The foundation of what follows can be found in Jurkevics (1988), Vidale (1988), Aster and others (1990) and Diez and others (2016). We cut our continuous three-component data into non-overlapping 1 h segments and express each segment in terms of eigenvalue ratios and particle motion directions (Fig. 4). These are computed for components i,  j = 1,  2,  3, time windows t and center frequencies f. The dot product of 1 h of component i,  j time domain seismic data X _i and X _j then results in the following 3 by 3 covariance matrix

(1)$$S_{i, j}^2 = X_{i}\cdot X_{\,j}^\ast $$

where the * indicates complex conjugation transpose. The eigendecomposition of this matrix is

(2)$$S_{i, j}^2 = Q \Delta Q^{-1}$$

where Q is the square 3 × 3 matrix whose ith column is the eigenvector v _i of $S_{i, j}^2$ and Δ is a diagonal matrix of eigenvalues E _i. The eigenvectors define the three orthogonal axes of motion averaged over every hour at every frequency and the corresponding eigenvalues are magnitudes of those eigenvectors. We define the principal eigenvalue ratio as seen in Figures 4a,c as:

(3)$$P = {E_1\over E_1 + E_2 + E_3}$$

Note that the eigenvalue ratio does not reflect the dominance of any particular component, which is advantageous, since there are often amplitude trade-offs between horizontal component resonance peaks as a function of time. Resonance peaks are demonstrably near-linear (i.e. P ~1, indicating that the vast majority of energy is accounted for by a single eigenvalue/eigenvector, meaning that particle motions are not elliptical, but only oscillate on a single axis). In Figures 4a,c, we thus have a normalized metric of energy independent of polarity that encompasses power spectral doublets seen on raw North and East components (Fig. 4e). From the primary axes defined by the eigenvectors, we further define θ = tan ⁻¹(V ₂/V ₁), the horizontal particle motion angle clockwise from north (Figs 4b,d), where V ₂ and V ₁ are the east and north components of the eigenvector corresponding to the largest eigenvalue.

Fig. 4. (a)–(b) Eigenvalue ratio and corresponding particle motion directions for the TIME power spectral data shown in Figure 2a for stations 13. (c)–(d) Similar to (a,b), but for RIS station DR06 for the spectral data shown in Figure 2b. (e) Peaks tracked across the north and east components of the data in Figure 2b at DR06, showing pervasively offset spectral peaks. Both peaks in a given pair are clearly evident in (c), and the corresponding splitting directions are constrained in (d).

2.4. Peak splitting methodology

Assuming resonances are dictated by both source and structure, an anisotropic variation in velocity should result in a shift of the resonance frequency, with the lower of two peaks in frequency representing the slow direction and the higher representing the fast direction, according simply to v = λf. If these resonances are caused by surface wave interactions with firn structures, then this assumption does not account for possible dispersion-related effects in velocity–frequency relationships. That noted, we are not inferring azimuthal anisotropy with depth in this study, but only anisotropy with respect to frequency. Particle motion directions corresponding to the slow and fast directions are defined as θ sampled at the time/frequency location of the picked eigenvalue peak.

As a first step, we transform the eigenvalue image into a sparse representation by picking peaks, where only peaks of a specified prominence (i.e. the minimum drop to the nearest trough) are picked and retained for each time bin of the image. For our data, prominence was set to 0.15 through trial and error aiming to balance picked peak numbers with noise (e.g. Fig. 4e). Given that resonance peaks have very high eigenvalue ratios (i.e. >0.8) and are well represented in this fashion, this step simplifies and greatly speeds up the calculation of offsets during the subsequent stretching process described below. For each trace, we start at the lowest identified peak and label it as a candidate slow direction. Upon flagging such a potential ‘slow’ peak, the following conditions are evaluated:

1. (1) The associated ‘fast’ peak is within 20% of a frequency stretch defined by $P_{\rm f} = P_{\rm s} + ( P_{\rm s}\;\ast\; S_{\% })$, where P _s and P _f are the candidate slow and fast peaks and $S_{\% }$ ranges from 0 to 0.2 in this study. This operation is similar to the arrival time stretching approach proposed in the previous section for active source data except that in this case the lower frequency peak of a given doublet becomes the ‘slow’ peak and the higher frequency peak is the ‘fast’ peak owing to v = λf arguments.

2. (2) Corresponding particle motions for P _s and P _f are suitably offset. Given that noise can contaminate split directions, we arbitrarily search for doublets with motions that split within $90\pm 25^\circ$ of each other (e.g. Fig. 5). Lowering this value results in a much smoother image, but progressively negates particle motions as a constraint. Making the range too tight in turn decreases the number of triggered detections, some of which are undoubtedly affected by noise contamination in the particle motions.

3. (3) Noisy signals can randomly result in perpendicular particle motions from one frequency bin to the next. Candidate peaks are thus computed in a moving window approach, where hour time bins ±4 of the current bin are used to compute the variance on both $S_{\% }$ and the fast/slow particle motions angles θ. A candidate peak doublet is accepted as anisotropic if its variance is smaller than $10^\circ$ for particle motions and 0.01 for $S_{\% }$. Once a pair of peaks are flagged as anisotropic, they are removed from the candidate pool.

Fig. 5. Example eigenvalue ratio trace showing three candidate peaks in the eigenvalue ratio with associated particle motion azimuth θ (angle clockwise from north). Each peak is ‘stretched’ within a range of 0–20$\percnt$ anisotropy, and its shifted frequency is compared to higher candidate peaks. If the variance over nine adjoining time bins of θ is suitably small and the mean falls between 90$\pm 25^\circ$, the pair is flagged as anisotropic and removed from the pool. In this example, the first two are flagged as anisotropic, but the third falls outside the stretching search range. θ _slow, θ _fast and the anisotropy magnitude are then summed into their respective images (Fig. 4) through a 2-D Gaussian estimated from the average statistics of peak widths.

An accepted pair of peaks thus yield three values: The magnitude of the peak offset assigned to the average of their frequency content and values of θ for the particle motions of the slow (lower peak) and fast (upper peak) directions. Each time a peak pair is accepted, a unit value, convolved with a 2D Gaussian representative of the average spectral peak shape (such as the peaks in Fig. 5) is added to the three separate matrices as shown in Figure 6a. As such, the images in Figure 6a are progressively ‘painted’ as an increased number of peaks and frequencies are sampled, thus building a map of fast and slow directions using inferred values of θ and splitting magnitude $S_\percnt$ at each station. Obvious maximums are manually (and smoothly where possible) picked in these images and are subsequently used as the basis for histogram-like representations of fast/slow axes.

Fig. 6. (a) Example anisotropy mapping for a 2-year period at DR06 through resonance splitting. In most cases, the lower frequency spectral peaks tend to be more stable over time, resulting in a higher degree of sampling in the images (brighter colors). When clustered peak sequences are temporally persistent, artifacts may occur in these images, given that sequential doublets may interfere with each other. In general, these effects are averaged out in the temporal stacking, but some may persist and require manual verification. (b) Fast/slow axes of anisotropy and corresponding percent magnitude at WAIS Divide.

Due in part to dominantly P-SV plane manifestations, Chaput and others (2022) proposed that resonance features are the product of Rayleigh wave amplifications. Despite persistently weak mapping of peaks onto the vertical component and strongly on in-plane horizontals, we cannot fully dismiss the possibility that some split peaks do not map to the vertical component, which may instead suggest S _H or Love wave resonances. If we are to interpret peak splitting in terms of azimuthal anisotropy, then S _V/S _V (with S _V or Rayleigh wave resonances) and S _H/S _H peak doublets (S _H or Love resonances) would be sensitive to azimuthal anisotropy, while S _V/S _H splitting, as noted in Diez and others (2016) for time-domain Rayleigh/Love inversions, would likely be sensitive to radial anisotropy inherent in the strong vertical gradient of the firn. A pervasive first-order feature of these resonance peaks largely precludes the second possibility, as the vast majority of observable peak doublets map to the vertical component as well as the horizontals. However, to avoid any possible misinterpretations, we take the additional step, described below, of ensuring that any two candidate peaks in a doublet also map, however weakly, to the vertical component. Under this interpretation, we verify that the first two variance tensor eigenvectors represent P-SV plane motion, with a vertical component 4 standard deviations above the noise level of surrounding spectral bands (which, as seen from Fig. 2, is the vast majority of the time). This somewhat arbitrary constraint undoubtedly caused the rejection of more data than necessary, but ensured that splitting was interpreted in the context of azimuthal anisotropy instead of (potentially) radial anisotropy.

Artifacts may however occur in situations where multiple peaks crowd the spectrum and overlap, which can result in 90$^\circ$ oscillatory uncertainty in inferred anisotropy at adjacent frequencies. A manual check was required to discriminate between consistent anisotropy directions (ensuring that peaks track each other over time) and those pairings resulting in false detections (artifacts). Efforts are currently underway to refine and automate this process, but for this study, frequency-dependent anisotropy curves were picked manually.

Variances in the 2D Gaussian function were set according to the spectral peak widths (effectively resulting in Gaussian smoothing of the image). More consistently sampled frequencies (i.e. time-persistent spectral resonances with low wander) were more strongly represented in the final images, along with the corresponding particle motion polarization for the fast (second peak in the doublet) and slow (first peak in the doublet) directions. In Figure 6a, artifacts in the vicinity of 25–30 Hz for the fast and slow axes are observable, indicating the presence of multiple clustered peaks and partial mis-mappings of fast/slow directions. Bright spots near ± 90$^\circ$ map to the same splitting direction.

3. Results

We discuss peak splitting estimates first at the TIME array near WAIS Divide as compared with active source estimates and then on the dense portion of the RIS array.

3.1. Anisotropy at WAIS Divide

As a first step, we process the active source data at WAIS Divide as described in the previous section to directly infer azimuthal anisotropy as a function of frequency. Slow/fast axes inferred with shot data are shown in Figure 6b, left panel, and display some variability with a modest change in azimuth around 50 Hz. Figure 6b, right panel, shows the resolved magnitude of anisotropy for three of the four resolved modes in their respective frequency ranges, starting at 2% at 15 Hz and increasing as a function of frequency. The abrupt increase in magnitude above 80 Hz coincides with a reduction of observable surface wave modes at those frequencies but this behavior is only observable at the highest mode (mode 4). That being said, the highest mode is quite visible at frequencies up to 150 Hz, thus suggesting that anisotropy may increase at very shallow depths.

At the lower end of resolvable frequencies in the shot gather (i.e. below 20 Hz, see Figs 3c,d), mode separation is difficult and thus the resolved fast direction is more or less the average of all modes and likely dominated by the lowest most energetic resolvable mode (mode 1 Fig. 3d).

We first apply the peak splitting algorithm to the circular array dataset at WAIS Divide and compare it to active source results. Peak drift and variability occur on multiple time/frequency scales as a function of atmospheric forcing, so arrays recording for longer periods (like the RIS array) tend to sample the total frequency range more thoroughly. The short recording period for the WAIS Divide array resulted in a narrow sampling of the possible frequency range compared to the RIS array, and this was further exacerbated by the sparser manifestation of peaks at WAIS stations.

For reference, Figure 2a shows an example of a spectrogram for the ~200 h dataset at station 13, and its corresponding eigenvalue ratios and particle motion polarizations can be seen in Figures 4a,b. Station 13 features a clear example of a peak doublet between 25 and 30 Hz with roughly orthogonal particle motion axes. The upward shift in frequency content in the middle of the deployment period is due to storm activity, the passing of which has been shown to have such an effect on peaks recorded on RIS as well (Chaput and others, 2018). For each station, we perform the procedure described in the previous section. The very short recording period and the small array aperture thus motivate the summation of contributions from all stations to form a single estimate of peak splitting at the array. Figure 7 shows the splitting magnitude and polarization as a function of frequency averaged over the WAIS array. Peaks are most robustly sampled between 18 and 30 Hz, with much poorer resolution at higher frequencies. Splitting magnitudes for the 18–30 Hz band fall in the 4–7% range (Fig. 7a, top panel), with a dominant particle motion that glides from ~175 to 220$^\circ$ from South (Fig. 7a, bottom panel). This would imply either depth-dependent anisotropy or some degree of unknown error causing this frequency-dependent behavior. This trend is also somewhat more tenuously observed by active source observations (i.e. a ~45$^\circ$ variation between 15 and 25 Hz) (Fig. 6b), though this could also be due to noise and instability inherent in the spectral stretching process. The poorly resolved lower frequencies (<15 Hz) in the active source gather at the small aperture TIME array (Figs 3c,d) do not however permit direct comparison with the lowermost of the resonance peak results.

Fig. 7. (a) Magnitude and fast peak direction average computed through peak splitting for the TIME array. (b) Radial histogram of fast and slow directions over all frequencies computed by peak splitting and by active source data.

We are primarily interested in evaluating how well splitting directions inferred from resonance peaks correlate to fast and slow axes of azimuthal anisotropy resolved by shots. As such, in Figure 7b we display our results as radial histograms for both the shot data (20–150 Hz, Fig. 7b lower panel) and the peak splitting (20–30 Hz, Fig. 7b upper panel). Comparison of their distributions indicates that the fast and slow directions of anisotropy are on average well represented by the first and second peaks in the doublets, though peak splitting clearly features more variability. We next present results for the significantly longer dataset gathered from RIS and test whether observable surface features such as crevasses and flow lines correlate with resonance peak splitting.

3.2. Observations on Ross Ice Shelf

The RIS array (Fig. 1b) consisted of RS (E-W) and DR (N-S) transects. RS was deployed as part of an effort to complement other large continental-scale oriented projects such as POLENET and TAMSEIS in terms of mapping crust and mantle structures beneath West Antarctica, while the motivation for DR was to study the dynamic response of RIS to gravity wave-induced vibrations. Both transects revealed a host of unusual signals (e.g. Chen and others, 2018; Zhao and others, 2019; Aster and others, 2021; Jenkins and others, 2021), including the firn resonances (Chaput and others, 2018). The large aperture RIS array featured a nominal station spacing of 80 km, but also included a denser sub-array of DR stations near a prominent flow convergence suture. The smaller sub-array station spacing potentially provided for more statistically robust processing of spectral peaks. Thus, we primarily focus our analysis on this portion of the array, though individual RS stations are presented as well in the supplement.

For each RIS station, we follow the processing that results in the mappings shown, for example, in Figure 6a. Manually tracking the obvious maxima results in corresponding magnitudes and fast/slow axes as a function of frequency, and these traces are shown in Figure 8 for the dense portion of the array. The magnitude of peak splitting in the vicinity of the denser DR sub-array reaches a peak between 4 and 6$\percnt$ on average in the 5–12 Hz range before largely tapering off at higher frequencies at a value of $\sim 4\percnt$. Fast/slow axes display complex directional behaviors at frequencies below 30 Hz, but abruptly stabilize at frequencies above that range.

Fig. 8. Fast and slow peak splitting directions for the dense DR array as a function of frequency (upper panel) and splitting magnitude (lower panel). The circles represent manually picked points well sampled by resonances in images such as Figure 6, and thus display a non-regular sampling of frequencies. The black trace in the lower panel is the average magnitude of anisotropy over the dense portion of the DR array.

4. Discussion

We discuss the following points:

1. (1) There is a strong qualitative alignment of fast/slow directions between active source and resonance approaches, and splitting results are compared with satellite imagery and Global Positioning System (GPS) motion datasets (Conway and Rasmussen, 2009; Matsuoka and others, 2011). The fast axis of anisotropy is overwhelmingly aligned with the direction of ice flow at frequencies above 30 Hz. There is no known mechanism in ice and snow that accounts for this observation, and as such we propose a qualitative explanation based on plastic deformation of a porous medium under uni-axial strain.

2. (2) Peak splitting obtained through our method clearly undergoes rapid transitions in terms of fast/slow axes as a function of frequency. We address this point by considering the multi-modal behavior of Rayleigh waves in realistic firn media and analyze their depth sensitivity through numerical analysis.

4.1. Comparisons with external datasets

Figure 9a shows the average slow and fast axes (referred to as shot slow and shot fast) inferred by active sources and the first and second peaks of the doublets (peak slow and peak fast) compared with the GPS-observed flow motions in the region (Matsuoka and others, 2011). Chaput and others (2022) noted that there exist structural variations within the aperture of the TIME array (for instance, the H/V ratio of the highest passively resolved Rayleigh mode is higher in the northeast quadrant of the array, suggesting a variation in site structure), but these would intuitively not manifest as such a variation in travel times (i.e. a period of π over travel time azimuths, Fig. 3b). As such, we assume that whatever layering effects were responsible for variations in Rayleigh particle motions and H/V in Chaput and others (2022) were second-order effects, and that these variations primarily result from azimuthal anisotropy. The average fast directions in both resonance peak splitting and active source surface wave studies are approximately aligned with flow directions away from the divide at frequencies above 18 Hz. Fast directions are best aligned with ice flow at equal Northing to the TIME array, suggesting that the observed resonances are more sensitive to local ice fabric effects than distant ice fabric effects. Interestingly, the fast directions resolved in Figure 7b are very well matched to the dominant eigenvector of crystal alignment inferred by recent polarimetric radar studies near WAIS Divide (i.e. ~19 degrees from North, Young and others, 2021). This fast direction alignment with ice flow runs counter to conventional wisdom of anisotropic mechanisms in moving ice.

Fig. 9. (a) Fast and slow directions for the shot data and the peak splitting analysis at the TIME array compared with the regional ice flow measurements inferred via GPS. Purple arrows represent GPS-inferred ice flow motions in the region (Matsuoka and others, 2011). The background images are MODIS surface elevation data. (b) Radial histogram depiction of fast peak splitting directions for the lower (red bars: 5–17 Hz) and higher (black bars: 30–50 Hz) frequency bands. Results from intermediate frequencies are omitted here to avoid clutter and are generally more weakly excited at the RIS array than the high and low bands. The green lines roughly depict advected crevasses as identified by Ledoux and others (2017).

Fahnestock and others (2000) and Ledoux and others (2017) meticulously documented surface flow lines observed by satellite to infer ice flow history and strain provinces within the Ross Ice Shelf, in general agreement with GPS observations and modeled flow in Klein and others (2020). The large-scale structure of the ice shelf is a product of multiple merging ice streams and glaciers merging at various velocities and directions, thus leading to a complex superposition of structural inheritance both laterally and with depth. The DR portion of the RIS seismic network used in this study spans two different strain provinces of RIS, thus allowing us to interpret our measurements with respect to a significant lateral change in ice shelf structure.

Figure 9b shows snapshots of fast splitting directions in various bandwidths for the dense DR sub-array alongside corresponding flow lines and other features detailed by Ledoux and others (2017). Although variations in splitting directions can occur at any point in the spectrum and are somewhat station dependent, we choose two representative frequency bands, 5–17 and 30–50 Hz, and display their fast directions in terms of radial histograms, with larger amplitudes representing more persistently sampled frequencies. We omit the intermediate frequencies (between ~18–30 Hz), given that this bandwidth frequency often features transitions in resonance peak fast/slow axes and thus histogram bins tend to be broad.

At longer wavelengths, azimuthal anisotropy should be sensitive to larger surface features such as exposed rifts and crevasses such as those visually identified by Ledoux and others (2017) from satellite images that likely propagate to greater depths within the ice shelf (Das and others, 2020). In particular, the dense sub-array is deployed near a ridge caused by the joining of distinct flows around the Crary Ice Rise (near RS17; Fig. 1b). Anisotropy directly on the ridge can be thus complex, where very local structures beyond flow directionality and obvious crevassing become important. West of the ridge, observable structures are largely dominated by crevasses that are perpendicular to the flow direction, while to the east, older advected crevasse fields are plentiful before large rifts gain prominence. At lower frequencies for DR stations, this appears to be reflected in the computed splitting directions, with multiple stations showing fast directions parallel to large visible rifts (RS03, DR14) or advected crevasses mapped by Ledoux and others (2017) (DR15,RS16). Stations deployed over the ridge show complex fast directions at low frequencies (5–17 Hz) that cannot be correlated to simple surface structures, while several stations away from the ridge show fast directions aligned with flow at all frequencies. Except for a few stations deployed on the ridge, nearly all stations robustly display fast directions that are aligned with flow direction at high frequencies (30–50 Hz).

4.2. Numerical modeling of surface waves in firn media

The strong layering, high-density gradient and low velocities of Antarctic firn allow for several strong Rayleigh modes to propagate between 5 and 50 Hz. Chaput and others (2022) detected up to threesuch modes in that frequency band with spatially varying particle motions and H/V ratios. Given that resonance peak frequency distributions were spatially correlated with other Rayleigh wave metrics, it stands to reason that abrupt changes in anisotropic axes as a function of frequency could be due to the depth sensitivities of the various Rayleigh modes being amplified, some of which may overlap significantly in frequency content. As such, we focus on studying perturbations in a smooth velocity model of firn determined via surface wave inversion on RIS (Diez and others, 2016) to see how depth sensitivities might vary as a function of frequency for multi-modal Rayleigh waves. This model is however only a rough reference, and the perturbations considered here are severe enough to encompass a range of realistic firn structures.

Surface wave depth sensitivity kernels in particular should provide insight into the various transitions in fast/slow axes as a function of frequency and show how various perturbations away from smooth, steady-state firn models as depicted by Ligtenberg and others (2014) could cause significant wavefield focusing around specific frequencies.

Rayleigh wave velocity is dictated by the depth-weighted average of its sensitivity kernel (e.g. Aki and Richards, 2002). If the subsurface presents particularly strongly layered structures, the depth sensitivity kernels of surface waves with respect to frequency can become highly focused (i.e. sharply peaked), or exhibit other complex behavior. Given what is known of firn structure from core profiling (Salamatin and Lisenkov, 2008), snow pits (Chaput and others, 2018) and ground-penetrating radar (GPR) (Arcone and others, 2005, 2016), we expect that strong shallow velocity and density discontinuities exist and are further overlain by theoretical firn densification transitions that take into account the progressive rearrangement of snow grains during compaction and individual grain plasticity effects (Alley, 1987; Reeh, 2008; Salamatin and Lisenkov, 2008; van den Broeke, 2008). Meticulous studies of elastic properties of snow (Frolov and Fedyukin, 1998) have revealed that both velocity and elastic moduli do not scale linearly with density, but rather vary according to piece-wise linear descriptions with prominent inflections.

Unlike longer period crust/mantle studies where global average velocity models are reasonably good approximations of true structure, directly relating a given surface wave frequency to its approximate depth sensitivity is a daunting task for firn. Given the multi-mode surface wavefield hosted by firn in the frequency range of interest (Chaput and others, 2022), it is important to study variations in depth sensitivity in these modes as a result of random realistic perturbations of assumed smooth firn velocity/density structure with depth, such as the one resolved by Diez and others (2016).

We studied this problem through direct numerical modeling via the finite element package RAYLEE (Haney and Tsai, 2017). The finite element approach provided by RAYLEE was preferred due to its flexible handling of water layers, its direct implementation in MATLAB and its robust representation of very shallow media at high frequencies. RAYLEE does not account for poroelastic effects (i.e. Biot theory) and indirect impacts of poroelasticity are described through extremely low densities and velocities. We began by considering the depth sensitivity of fundamental, firstand second mode Rayleigh waves propagating in a medium defined by the smooth velocity profile provided by Diez and others (2016) near station RS04 (Fig. 1b) and progressively slowed the profiles in subsequent models. Sensitivity kernel results for Rayleigh modes 0, 1 and 2 are shown in Figure 10. For the initial smooth profile (Fig. 10, top row), depth sensitivities of the first three Rayleigh modes also vary smoothly, with increased sensitivity to the near surface as frequencies increased. This effect is greatly magnified in slower profiles, until all modes become almost uniquely sensitive to the near-surface at all frequencies, with the exception of a highly focused locale at the elbow of the lowest frequency range for each mode. These models are extreme, and we do not anticipate that firn profiles on RIS or WAIS Divide could yield such generally low velocities. That being said, firn profiles in areas with less precipitation and lower temperatures (for instance South Pole) have shown evidence of significantly lower velocities in the near surface coupled with a deeper firn in general (van den Broeke, 2008), so this class of behaviors is to a degree plausible.

Fig. 10. Finite element simulations of Rayleigh wave sensitivity kernels for progressively slowed versions of the uppermost RIS velocity model of Diez and others (2016). The top row is the unperturbed model. The hard cutoffs at lower frequencies are numerical cutoffs rather than physical.

As mentioned above, GPR, borehole, snow pit studies and non-steady-state firn densification models have indicated that abrupt back and forth transitions in firn parameters may exist at any depth, but are particularly prevalent in the shallow regions, where settling and melt processes dominate and ice lens and hoar frost layers have been broadly observed (Reeh and others, 2005; Ligtenberg and others, 2014; Alley and others, 2018). Furthermore, snow features several empirical transitions in the slope of the velocity/density relationship (Frolov and Fedyukin, 1998) in the vicinity of 150, 340, 550, 720 and 830 kg m⁻³, representative of hypothesized transitional steps within the process of firn densification.

In order to study the potential depth sensitivity effects of shallow firn structures, we created five realistic models featuring large local excursions dictated by adding low-pass Gaussian noise to the smooth RIS model reported by Diez and others (2016), reproducing the results of Chaput and others (2022) (Fig. 11). In effect, these perturbations aim to qualitatively emulate the presence of strong near-surface layering with large swings in parameters, as predicted by non-steady-state firn densification models (Reeh, 2008; Ligtenberg and others, 2014). Where smooth profile shows corresponding smooth kernels for the first three modes (Fig. 10), the addition of strong noise perturbations in Figure 11 creates complex focusing patterns in the near-surface. We note that certain modes focus at very narrow depths, but large bandwidths, which may explain why spectral peak splitting axes often remain stable over a certain frequency range before suddenly changing (e.g. Fig. 8). This implies that a simple frequency mapping to depth is inherently difficult without additional information such as ice coring or GPR data, or which specific wave mode is being excited.

Fig. 11. Example finite element simulations of Rayleigh wave sensitivity kernels for noise perturbed firn models that introduce random layering, using the upper RIS velocity model of Diez and others (2016) as a reference (Fig. 10).

Upon inspection, these perturbations result in greatly increased focusing of sensitivity kernels to zones in the velocity model with rapid transitions between high positive and high negative velocity anomalies. Such zones can represent transitions from normal firn to hoarfrost to ice, or any similar permutation. The large variability in sensitivity of all three modeled modes may explain the apparently random distribution of peaks, if they are indeed due to the Rayleigh wave amplifications suggested by Chaput and others (2022).

This set of behaviors would explain several peculiar observations in our data. First, there is often a band gap between 12 and 20 Hz observed in our dataset (particularly at RIS stations), which could represent, on average, the falling off the fundamental Rayleigh and the emergence of the first higher Rayleigh mode for slightly slower average firn profiles (Fig. 10). Second, there is a sharp transition in the azimuthal direction of dominant anisotropy in those two frequency bands (5–12 and 25–50 Hz), as mentioned in the previous section, given the abrupt change in depth sensitivity caused by near-surface parameter ‘switch-backs’.

The sharp transition in Figure 8 near 30 Hz may be indicative of such a hand-off in sensitivity, with the consequence that frequencies in the 30–50 Hz range are uniquely sensitive to the top few meters of firn through the dominance of a focused sensitivity kernel. Naturally, without direct knowledge of the velocity profile, it is difficult to map these determinations of anisotropy to depth directly. These numerical models also do not however offer any insight as to why our inferred fast directions, for both active source and peak splitting methods, tend to align with the direction of ice flow at high frequencies. We discuss this below.

4.3. Mechanisms for anisotropy in snow and ice

Glacial ice can exhibit anisotropy through three primary mechanisms: crystal alignment due to plastic accommodation of stress (e.g. Diez and others, 2014), directional fracturing due to non-isotropic stress fields (e.g. Zhan, 2019), and fine layering inherent in the firn densification process (Diez and others, 2016). Early work by Bentley (1972) using active source data characterized deep ice in the West Antarctic Ice Sheet as having a depth-dependent c-axis with a mostly vertical orientation and a tilt normal to flow lines in the area. The fast direction in this interpretation corresponds to the c-axis, i.e. the symmetric alignment of the ice crystals in the direction of maximum compression. Laterally confined flow, e.g. in the case of two convergent ice streams (e.g. Smith and others, 2020), would result locally in the rotation of c-axes toward the compression strain axis (i.e. perpendicular to flow). This is potentially observed at low frequencies for stations RS04 and DR11 deployed over the ridge (Fig. 9b), while DR05 and DR04, also deployed on the ridge, show a mix of behaviors. Note that inherited fabrics can however cause significant deviations from this assumption (e.g. Budd and Jacka, 1989). Ice shelf accumulation and seaward flow (which is ~940 m  a⁻¹ at DR10 and may exceed 1 km a⁻¹ in the seaward regions of RIS, Klein and others (2020)) record features that reflect a rich ice history and may thus create multiple local depth-dependent anisotropic structures that can manifest as azimuthal in nature.

The addition of firn layering complicates matters significantly. C-axis alignment in ice crystals does not explain along-flow fast directions at higher frequencies for the majority of RIS and WAIS Divide stations, since the top several 10 s of meters of the ice shelf do not undergo significant gravitational or other compression, and the accelerating flow of ice toward the coast would lead to flow-perpendicular and not flow-parallel maximum compression (i.e. extension along flow results in compression perpendicular to flow). We explore the peculiar physical properties of porous firn and snow as an explanation.

Whether brittle crevasse formation or ductile deformation dominates anisotropy at a given depth relates directly to the broad field of fracture formation. Under simple extension, ice is known to reach the ductile to brittle transition (d/b) at strain rates in the vicinity of 10⁻⁶ s⁻¹ (Kirchner and others, 2008). Snow, however, reaches this transition in a much more complex manner, with dependencies on temperature, density and structure. Studies by Narita (1984) and Kirchner and others (2008) provide the only systematic laboratory investigations of the (d/b) transition for tensile failure of snow at several densities and temperature, noting values generally ranging from 5 × 10⁻⁴ s⁻¹ to 9 × 10⁻³ through their search range. Specimens were averaged between 250–350 and 350–450 kg m⁻³ and thus, for RIS and WAIS Divide, are only valid as approximations for the near-surface of the ice shelf. Average strain rates for RIS are estimated to be on the order of ~1 × 10⁻¹¹ s⁻¹ (i.e. for unidirectional spreading of 1 km a ⁻¹ over a 600 km ice shelf), and as such, purely from tensile stress inherent in the acceleration of ice shelf flow, brittle failure causing perpendicular fracturing far away from the calving line should not occur in isotropic ice and much less so in snow and firn. Fracture formation thus depends on a plethora of other factors, the broad study of which can yield promising estimates of fracture probabilities (Emetc and others, 2018). Whereas crevasses may heal in ice under compression, those caused by ice shelf flow over asperities persist, advected or otherwise, for the entire flow length of the ice shelf (Ledoux and others, 2017). At densities <550 kg m⁻³, the main densification process in firn is settling, whereby the pore space is progressively closed through rearrangement of snow grains (van den Broeke, 2008). Persistent fractures are thus unlikely to persist for long in snow at the aforementioned densities under tensile stress, which, for RIS and WAIS Divide, occur at depths of under 8–15 m (van den Broeke, 2008).

Thus, given the previous calculations, excepting large open crevasses, we do not expect anisotropy in the top 10 m to be dominated by flow-perpendicular or advected fractures.

This conclusion is supported by the following arguments:

1. (1) Given that the top 15 m of firn consists primarily of an open pore space with a very low bulk modulus (Frolov and Fedyukin, 1998), that the critical strain rate for the d/b transition in snow is several orders of magnitude larger than ice and many orders of magnitude larger than typical RIS strain rates (Kirchner and others, 2008), and that near-surface densification occurs primarily through settling, stresses generated by the acceleration of ice shelf flow toward the ice shelf front must be accommodated plastically. Although deeper ice may feature vertically aligned c-axes in ice crystals (Bentley, 1972), we do not expect this mechanism to dominate in surface snow unless there is significant lateral compression.

2. (2) Figure 8 clearly shows evidence of shifting behaviors between low- and high-frequency spectral splitting (5–17, 30–50 Hz) for the dense sub-array. This strongly suggests a relatively abrupt transition in depth sensitivity for different frequencies and an abrupt change in the way anisotropy is expressed with respect to depth. Thus, there would have to be a typical depth at which perpendicular fracturing is preserved once created and thus dominates anisotropy. Above that depth, firn densification via settling results in purely plastic stress accommodation. The likely multi-modal character of the firn resonances however suggests that this transition in sensitivity could also be strongly tied to surface wave modal hand-offs (Chaput and others, 2022), since each mode can present very different depth sensitivity kernels at a given frequency.

A natural explanation for these observations lies in the properties of porous media (Biot, 1962; Sidler, 2015) and metrics relating wave velocity to tortuosity, which is a measure quantifying the ratio between the actual flow path length between two points and the Euclidean distance between those points, in an open foam under strain (Kirchner and others, 2008). For bi-laterally constrained growth of an open-celled foam (i.e. growth only occurs along one axis), Melon and others (1995) and Melon and others (2008) noted that there exist significant differences in axial and longitudinal tortuosity and this results in large increases in phase velocity for waves propagating along the growth axis. Due to the stress regime, the pore space and solid lattice become aligned with the growth direction. As such, the effective tortuosity of the pore space is reduced in the direction of elongated pores because the corresponding solid lattices are more linear in that same direction, and the distance along the solid lattice between any two points is shorter in the elongated direction than it is in the perpendicular one. Similarly, open-celled bone media feature longitudinally aligned porous structures (referred to as trabecular structures) that also result in the same form of axial anisotropy, with waves propagating parallel to the elongated pores featuring much higher velocities as a consequence of reduced tortuosity in the pore space (Hosokawa and Otani, 1998; Lee and others, 2007).

Although we could not find a direct laboratory example of anisotropic wave propagation for high porosity snow under strain, presumably due to the relative infancy of crysoseismology (Podolskiy and Walter, 2016) and the difficulty inherent in studying snow in the laboratory setting while emulating its field conditions, it is conceivable that this mechanism applies here. As mentioned above, snow deforms plastically under estimated strain rates on RIS, and brittle fracture is not expected to be a dominant mode of stress accommodation for the shallow-most firn. Given that snow grains sinter over time scales of hours to days (Peinke and others, 2019), the inherent cohesion of snow even at the top most layers would thus imply gradual plastic extension of the pore space in the direction of strain. van den Broeke (2008) estimated the age to pore close-off in the firn on Ross Ice Shelf and the general WAIS Divide region to be ~200–500 years. Densification processes imply that moduli for snow are expected to increase very rapidly with depth. Thus, down to pore close-off, the firn can have accumulated several hundred years of plastic strain along flow.

The degree to which progressive pore space reduction would affect tortuosity is unclear, but it stands to reason that there should be a specific depth range in the firn that features maximal anisotropy, i.e. where strain accumulation has resulted in maximal pore space lengthening without losing too much of the pore space to reduce tortuosity. In Figure 8, the lowest frequencies sampled by firn resonances fall between 4 and 8 Hz, and there is indeed a maximum splitting offset shortly above those frequencies and a gradual decrease to $\sim 4\percnt$ between 7 and 16 Hz. Although this explanation is attractive, it is however complicated since the fast direction of splitting at lower frequencies tends to be more variable and often aligns with visible surface structures. Therefore, as described above, anisotropy at depths >10–15 m might be dictated by up to three separate sources (c-axis alignment, fractures/crevasses and flow-generated tortuosity variations), and the available RIS measurements do not allow for systematic discrimination of mechanisms other than qualitatively.

5. Conclusion

In summary, we have presented a novel situation where frequency-dependent local anisotropy can be determined from single station recordings of ambient noise without resorting to time domain reconstructions of surface waves. This unusual situation is made possible by the presence of pervasive spectral peaks recorded in continuous data that are trapped in the relative near surface and exhibit strong azimuthal splitting visible in the three-component seismic recordings. Direct tracking of frequency-dependent splitting is performed for 2 years of data at 34 broadband stations spanning RIS, and the magnitude and fast/slow directions are systematically mapped and interpreted. On average, the magnitude of anisotropy appears to reach a maximum in the 5–12 Hz range. Anisotropic directionality at these frequencies can be complex, but generally reflects large-scale features with surface expressions, such as fractures perpendicular to current flow directions, or rotated and translated fractures from past flow directions. Fast directions of anisotropy in those cases are generally parallel to these features. The magnitude of anisotropy generally falls to $\sim 4\percnt$ at frequencies above 30 Hz, with a fast direction aligned with ice flow. This change in behavior can be attributed to a different mechanism, namely the plastic accommodation of strain in the highly porous firn that produces elongated pore spaces and decreased tortuosity in the flow direction by direct stretching of a porous structure. The transition, in frequency, of these two behaviors is however not smooth in the terms of depth sensitivity. Numerical models of Rayleigh wave sensitivity in realistic firn models with finite element methods support this contention, as large parameter swings in non-steady-state densification models can create highly localized surface wave sensitivity kernels. To our knowledge, no anisotropy analyses of firn core samplings exist due to the difficulty of retaining core integrity, and this should be attempted in future experiments.

Given the agreement between large-scale features observed by satellite imagery, the inferred flow directions and both low- and high-frequency computations of anisotropy, there is strong evidence that these novel observations can be useful to infer strain regimes of at least the top 50 m of firn that are capable of hosting pervasive spectral resonances of this type.