## 1. Introduction

Firn is present across approximately 99% of glaciated Antarctica (van den Broeke, Reference van den Broeke2008; Ligtenberg and others, Reference Ligtenberg, Helsen and van den Broeke2011). Its thermal and mechanical properties are thus critical in moderating external environment influences upon underlying glacial ice, including ice sheets, valley glaciers and ice shelves. High air content and consequent thermal insulation, high albedo, mechanical ductility, ability to shape-shift under the influence of wind and precipitation and thus heal surface disturbances, and a capability to absorb and re-freeze surface meltwater all enhance this protective role. Should the firn layer becomes compromised, underlying ice becomes susceptible to degradation from multiple processes, including ablation, surface melt accumulation, hydrofracture and albedo reduction (Kuipers Munneke and others, Reference Kuipers Munneke, Ligtenberg, van den Broeke and Vaughan2017).

Seismologically, firn structures create strong spatial gradients and/or sharp seismic impedance contrasts. Such features introduce strong and sometimes complex seismic site responses and localized wavefield effects (e.g. Aster and Shearer, Reference Aster and Shearer1991; Vernon and others, Reference Vernon, Pavlis, Owens, McNamara and Anderson1998; Wang and others, Reference Wang2020). In a parallel sense but at a vastly different scale, sedimentary basins are also characterized by low seismic velocity structures with sharp basal and marginal contrasts that produce strong resonances. Initial studies of horizontally polarized shear (SH) wave resonances in basins were classically depicted by Weinstein (Reference Weinstein1969) and were further expanded upon by Rial (Reference Rial1989) to consider elliptical and cosine triaxial geometries. These studies revealed that resonances in a 3D basin may not be accurately estimated from 1D or 2D counterparts. Recent works by Zhu (Reference Zhu2018), Poggi and others (Reference Poggi, Ermert, Burjanek, Michel and Fah2015), and Roten and others (Reference Roten, Fah, Cornou and Giardini2006) have explored methods to recover dominant resonance modes in geometrically realistic basins, and have introduced spectral decomposition methods to identify dominant wavefield components. Basin structures also introduce surface wave amplifications that can manifest as potentially strong wavefield effects that lacks closed form solutions. From an engineering perspective, classic estimations of long period surface wave basin response based on energy flux conservation typically rely on comparisons with hard rock benchmarks (e.g. Tromp and Dahlen, Reference Tromp and Dahlen1992) to produce relative measures of amplification. Although this approach has been shown to be valid for short periods (Bowden and Tsai, Reference Bowden and Tsai2017) and well-defined basin edges (Brissaud and others, Reference Brissaud, Bowden and Tsai2020), more complete predictions of surface wave amplifications, particularly for 3D structures, require full wavefield numerical methods (Komatitsch and others, Reference Komatitsch2004).

Low-velocity near-surface structure can produce complex and sometimes unexpected wave propagation phenomena in the seismic response and excitation spectrum, and careful processing and attention must be paid to wave type when using ambient noise to image or otherwise infer such structures (Priestly and others, Reference Priestly, Orcutt and Brune1980; Savage and others, Reference Savage, Lin and Townend2013; Ma and others, Reference Ma, Clayton and Li2016). For instance, Ma and others (Reference Ma, Clayton and Li2016) noted the presence of a ‘modal hand-off’ from fundamental to first mode Rayleigh waves occurring at an obvious threshold frequency for ambient noise correlations performed in the Los Angeles basin. They concluded that fundamental mode Rayleigh waves within the basin floor showed increased horizontal to vertical energy ratios (*H*/*V*) at a threshold frequency, above which the low velocities of the basin fill allowed for the emergence of higher Rayleigh modes. Whereas long period fundamental mode Rayleigh waves observed at the surface commonly show retrograde motions for most Earth structures, there exist examples of transitions to prograde motion for shallow low-velocity features (Berbellini and others, Reference Berbellini, Morelli and Ferreira2016; Tanimoto and Rivera, Reference Tanimoto and Rivera2005) that have been confirmed through numerical eigenfunction modeling (Tanimoto and Rivera, Reference Tanimoto and Rivera2005; Denolle and others, Reference Denolle, Dunham and Beroza2012). Denolle and others (Reference Denolle, Dunham and Beroza2012) noted that one or more particularly strong velocity gradients appear to be a prerequisite for this behavior.

We analyze analogous basin-like wave propagation at high frequencies (typically >5 Hz) within the Antarctic firn zone as it grades sharply from snow to glacial ice across a depth range of tens of meters. The densification process of snow accumulation, settling and compaction leads to laterally homogeneous and very strong vertical near-surface seismic velocity and density gradients (Reeh, Reference Reeh2008; Reeh and others, Reference Reeh, Fisher, Koerner and Clausen2005; Li and Zwally, Reference Li and Zwally2011) with *P*-wave velocity increasing from the atmospheric sound range (~300 m/s) to near 3000 m/s across ~50 m (e.g. Diez and others, Reference Diez2016). As a result, snow and firn media underlain by glacial ice present a variety of wave propagation effects that can be viewed as higher frequency analogs of basin phenomena. An additional behavior of the firn system is that seismic velocities in the highly porous near-surface snow layers can be so low that *P*-waves in the solid lattice may travel more slowly than associated sound waves propagating in the open pore space (Sidler, Reference Sidler2015). We study the ambient seismic excitation of Antarctic firn due to wind and anthropogenic sources using a passive three-component nodal seismograph array deployed near West Antarctic Ice Sheet (WAIS) Divide camp, Antarctica (elevation ~1800 m). We demonstrate that highly localized fundamental Rayleigh particle motion reversals, higher mode emergence and trade-offs, spatially variable Rayleigh *H*/*V*, and spatially variable absolute wavefield amplifications are pervasively observed at frequencies above 5 Hz. We also propose a revised interpretation of ambient absolute spectral amplifications recently documented by Chaput and others (Reference Chaput2018) characterized in a microbasin resonance framework. Such modes, processed at isolated geographically widespread stations, can constrain regional impacts of near-surface atmospheric wind, melt and temperature forcing on environmentally vulnerable cryospheric media. Increased understanding of these seismological phenomena more generally advances prospects for single and multi-station methods for inferring structures and monitoring temporal changes in cryospheric media.

## 2. Ambient seismic excitation at the WAIS Divide site

### 2.1. Ambient noise interferometry

The broad use of ambient noise methods over the past few decades has sparked a revolution in imaging and temporal monitoring in seismic settings at a multitude of scales. In short, works by Snieder (Reference Snieder2004), Wapenaar and Fokkema (Reference Wapenaar and Fokkema2006), Wapenaar and others (Reference Wapenaar, Draganov, Snieder, Campman and Verdel2010) and Campillo and Paul (Reference Campillo and Paul2003) demonstrated that under certain conditions, the cross-correlation of multiply scattered or otherwise isotropically excited wavefields at two receivers results in an approximation of the impulse response function between those two receivers. Given that at lower frequencies the seismic wavefield is either isotropically excited or at least predictably excited by the presence of a coastline, this relationship has been broadly exploited to reconstruct surface waves (Ritzwoller and others, Reference Ritzwoller, Lin and Shen2011) and body waves (Nakata and others, Reference Nakata, Chang, Lawrence and Boue2015) for a number of imaging applications. If the wavefield is isotropic then the following holds:

Here, *u*(*x* _{A,B}, *t*) denotes the wavefield recorded at receiver A or B, the * denotes a convolution, with the time reversal turning it into a correlation, and *G*(*x* _{A}, *x* _{B}, *t*) represents the Green's function between receivers A and B. *S*(*t*) is the autocorrelation of the noise and essentially represents the wavelet convolved into the Green's function. A given wavepath in the Green's function is correctly reconstructed if there are sources distributed along its stationary zone (Snieder, Reference Snieder2004). For instance, a correctly reconstructed Rayleigh wave between two surface receivers requires surface sources distribute along the end points of the station pair line. We rebuilt Rayleigh waves through anthropogenic noise in the following study, and as such, we generally limited our study to interstation paths that aligned with a known anthropogenic noise source, WAIS Divide camp.

### 2.2. Noise correlations on the TIME array

As a first step to identify the nature of seismic waves supported by the firn and underlying glacial ice, we performed daily average cross-correlations of ambient noise between all interstation component pairs. Cross-correlation tensors were constructed with 1 h windows and 75% overlap (Seats and others, Reference Seats, Lawrence and Prieto2012). Of the ~12 days of passive data, 8 January featured no local human activity in the immediate array vicinity (i.e. no drilling), and thus the strongest noise source at that time was expected to be the WAIS Divide camp. Cross-correlations were computed in sub-bands of ~5 Hz at frequencies 3–45 Hz with frequency steps of 1 Hz. The resulting correlation tensors were used to beamform the velocity and azimuth of incident seismic energy across the array. We did so by taking the envelope of the correlation functions and by assuming plane wave incidence on the array. Employing a grid search of velocity and azimuth, we used the envelope of the correlation tensor diagonal (same component correlations) to quantify the least squares amplitude discrepancy with respect to synthetic waveforms, created with simple Gaussian envelopes with frequency-matched standard deviations. We further isolated Rayleigh and Love modes by separately considering the P-SV and Transverse planes.

Figure 2 shows example beamforming results in correlation gathers on 8 January (Figs 2a,b), during which WAIS Divide camp was a strong noise source (i.e. no local human noise) and on 12 January (Fig. 2c), during which both camp noise and local drilling activity were present (Fig. 2b). Although the dominant azimuth is the same in Figs 2a,b, there is a clear velocity increase at higher frequencies (above ~17.5 Hz). Figures 2d,e show results for the maximum amplitude in the beamforming, denoting group velocity, at different frequencies for all days, which feature multiple directions of dominant noise due to local drilling efforts and the median apparent velocity at sampled sub-bands. A distinct velocity transition can be observed at 15—20 Hz, where the dominant reconstructed arrivals in the correlation functions abruptly decrease in slowness by ~0.5 s/km. We used 8 January as a benchmark to analyze the finer details of this velocity increase and rotated the array into the direction of WAIS Divide camp (so that the radial components point toward the camp).

This velocity step between 15 and 20 Hz is suggestive of a surface wave modal transition in which a higher mode becomes the dominant carrier of energy. Dispersion is also observable, and particle motions are mostly limited to the radial-vertical plane at frequencies below 30 Hz, consistent with Rayleigh waves. Above 30 Hz, the noise wavefield diversifies somewhat in its angular incidence, but the azimuth continues to dominantly point in the direction of WAIS Divide camp. To more clearly visualize mode separation and transitions, we employ a frequency/slowness (FTAN) technique over horizontal ray paths within ± 12 degrees of the camp azimuth for 3—45 Hz energy (Fig. 2f). A second higher mode emerges above 30 Hz, evidenced by a distinct second modal transition. Variations in modal amplitudes could be due to varying contributions of wind noise and/or amplification/attenuation related to medium properties such as porosity.

### 2.3. Multi-mode particle motions

Given the extremely low velocity and layered character of firn, we approached the study of these waves noting prior results such as Ma and others (Reference Ma, Clayton and Li2016), where strong velocity transitions due to the presence of sediments resulted in a modal transition between fundamental and higher mode Rayleigh waves. In Ma and others (Reference Ma, Clayton and Li2016) the higher mode only propagated at frequencies sensitive to the basin fill, and so could be separated by analyzing its frequency-dependent particle motions. For typical seismic media, fundamental mode Rayleigh waves for Earth models have retrograde motions at the surface, whereas first higher mode Rayleigh waves (e.g. Sezawa waves; Sezawa (Reference Sezawa1935)) mostly display prograde motions.

Figure 3 shows examples of particle motions for the reconstructed correlation gather centered at 3 Hz (Fig. 3a) and 7.5 Hz (Fig. 3b) for 8 January, with the array rotated into the camp azimuth. Rayleigh particle motions for all gathers below ~18 Hz are unusually linear and inclined, with the vertical and radial components having approximately equal amplitudes and being close to in-phase. This would generally suggest *P*-wave energy incident at ~45 degrees, or *S* _{V} waves, but the evident dispersion and the velocity jump at 17—22 Hz are inconsistent with this interpretation. Ma and others (Reference Ma, Clayton and Li2016) suggested that inclined ellipticity, as observed here, may indicate the presence of contaminating body waves (i.e. reflected/refracted *P*, *S* _{V} in the firn medium; Albert (e.g. Reference Albert1998)). We briefly recap the results of the Ma and others (Reference Ma, Clayton and Li2016) derivation here. The summation of a single recorded surface wave and a body wave with arbitrary phase *ϕ* and frequency *ω* can be written as

where **M** is a 2 by 2 matrix related to the amplitudes and phases of both waves and **UΣV**^{T} is the singular value decomposition of **M** (see Ma and others (Reference Ma, Clayton and Li2016), Appendix B). As such, if body wave arrivals exist concurrently with surface wave arrivals, we expect circular particle motions defined by [cos(*ωt* + *ϕ*), sin(*ωt* + *ϕ*)]^{T} to be stretched by the **Σ** operator and rotated by **U**. Particle motion linearization is known to happen near natural *S*-wave structural resonance associated with shallow discontinuities, and directional body wave contamination could result in a rotation of those motions. In our case, both the radial and vertical components of motion are of similar scales for most reconstructed wave intervals in 3—18 Hz, which would imply significant body wave contamination combined with the linearization of motions due to structural resonance. This unusual phase and amplitude matching of the radial and vertical components evades any other obvious explanation, and should be the subject of deliberate future studies, should this phenomenon prove to be pervasive in firn media at higher frequencies than typically studied for imaging purposes.

### 2.4. Rayleigh wave H/V ratios

Along with particle motion determinations, we also estimate *H*/*V* for the three observable Rayleigh modes (referred to here as Rayleigh 0, Rayleigh 1, Rayleigh 2; Fig. 3c), calculated by picking the moveouts of obvious radial/vertical waves shown in Figure 4. In general, it is clear that the camp does indeed provide the strongest source of incoming wavefield energy, as cross-correlations are notably one sided. What this means, in terms of Eqn (1), is that the wavefield is not isotropic (i.e. it had a dominant direction), and that we are not recovering both terms in the right-hand side. Rather, we are correctly rebuilding the surface wave component for one of the two Green's function terms and little else. For frequencies between 18 and 30 Hz, there appears to be a secondary source of energy resulting in a partial reconstruction of the second Green's function term (e.g. Fig. 4, 27 Hz). This implies either a band limited source of noise such as wind, or persistently reflected waves from buried or topographic features outside the array. Though interesting, this point is left to future study.

*H*/*V* as a function of frequency is generally calculated as $\sqrt {{( R^2 + T^2) \over Z^2}}$, where *R*, *T*, *Z* are Fourier amplitude for frequency-dependent windows centered around a given moveout and are calculated for interstation paths within ± 12 degrees of the WAIS Divide azimuth (resulting in 29 paths). However, since we are isolating Rayleigh waves for paths satisfying the stationary phase criterion (i.e. pointing toward WAIS Divide camp), Rayleigh waves will be confined the P-SV plane, and corresponding Rayleigh *H*/*V* ratios simply become |*R*/*Z*|.

*H*/*V* shows strong variability across modes (Fig. 3c), with Rayleigh 0 showing an *H*/*V* near unity (e.g. Fig. 3b), Rayleigh 1 being vertically dominated, and Rayleigh 2 showing particularly elevated *H*/*V* that are clearly path dependent and with peak amplitudes that can vary by as much as one order of magnitude. We focus on Rayleigh 2, which features dominantly radial energy, particle motion linearization, and a high *H*/*V* with clear spatial variability, indicative of a relative lack of body wave contamination presumably due to source mechanisms inherent in the noise source at WAIS Divide. Mapping the smoothly interpolated *H*/*V* peaks for Rayleigh 2 to their respective paths yields a distinct pattern of lower peak frequencies in the upper portion of the array. *H*/*V* is well known to be sensitive to shallow-most structural layering (e.g. Pina-Flores (Reference Pina-Flores2016)), and is likely indicative here of firn structure variability across the array. Peak *H*/*V* is generally assumed to decrease in frequency content with a deepening first discontinuity, although isolating the specific Rayleigh modes responsible for the behavior is generally difficult when dealing with single stations. Given that we are able to isolate three distinct Rayleigh modes through cross-correlation, such an analysis is possible, and we describe it in the following section.

For noise definitively originating from WAIS Divide camp (e.g. limited to Rayleigh 0 and 1 on 8 January), we also spatially display particle motion directions and ray paths at various frequencies for WAIS near-azimuth paths. For frequencies between 5 and 15 Hz, where Rayleigh 0 dominates, we observe persistent spatially dependent transitions from retrograde to prograde motion across the array (Fig. 4, panels 2–5), and more complicated transitions for the higher modes as a function of frequency. The presence of both prograde and retrograde Rayleigh 0 waves has been previously observed in basin settings (e.g. Boue and others, Reference Boue, Denolle, Hirata, Nakagawa and Beroza2016), but this is to our knowledge the first time this has been observed to transition across such a small spatial range and at such high frequencies. Both Rayleigh 2 *H*/*V* and Rayleigh 0 particle motions share the same spatial pattern with reversed particle motions, indicative of similar structural influences, and pointing to variable near-surface structure under the array. The more chaotic appearance of particle motions of higher modes may be due to the low samples to wavelength ratio of these modes, and the quasi-linearity of radial/vertical plane motions throughout. That being said, there is a general particle motion reversal between, for example, the spatial pattern at 12 Hz and that at 33 Hz.

As a last observation, we study manifestations of recently documented firn resonance peaks (Chaput and others, Reference Chaput2018) which, unlike *H*/*V* ratios, are not relative amplifications but rather are absolute in nature and currently evade simple explanations.

## 3. Ambient spectral amplifications

Locally and temporally variable absolute amplifications in ambient spectra were observed here above 5 Hz, and we compare them to recently documented narrow band firn resonances first observed on Ross Ice Shelf (RIS) (Chaput and others, Reference Chaput2018; MacAyeal, Reference MacAyeal2018). Continuous passive data recorded at broadband stations deployed on RIS showed narrow band resonances that dynamically responded to environmental forcing effects on scales of minutes to months, including high-temperature anomalies and wind and accumulation-driven surface deposition/erosion effects. Although these observations showed some coherence consistent with regional weather data, stations were too isolated and distant from point weather measurements to analyze more subtle and short-term spatial and temporal variations in the seismic response. Passive data from the TIME array, however, allow more insight into this apparently widespread firn environment phenomenon, and several interesting temporal and spatial trends in peak patterns are recorded on this array.

Given our circular array geometry, we are able to make a number of diagnostic observations concerning the nature of the observed frequency peaks. We first study the spatial variability of the prominent modes as a function of angular location. Figure 5 shows radial component spectrograms for the first 70 h of data at all stations. We notice that adjacent stations often have similar spectral patterns, and that there exists a roughly 2π variation pattern in the peak frequency (red dotted curves in Fig. 5). At some stations, however (e.g. 16, 20, 21, 32), the pattern jumps between adjacent stations (black dotted lines in Fig. 5), with frequency values approximately twice that of neighboring averages (i.e. a factor 2 in terms of peak frequency with respect to its neighbors). These sudden jumps suggest the excitation of a higher harmonic. The smooth spatial variability of the lower pattern, interrupted by the occasional higher jump, suggests a sensitivity to firn structure variability over the array aperture, as also indicated by other observations noted in previous sections. For instance, if these resonances are sensitive to shallow structure (in the *H*/*V* sense), then a 2π variability would imply an approximate 2D dipping plane in terms of structure across the array. We evaluate this hypothesis in the following section. A persistent 60 Hz peak is recorded uniformly across the array, which we attribute to electrical power generation and use at WAIS Divide camp.

To ensure that these observations are indeed part of the temporally varying and environmentally linked phenomena documented by Chaput and others (Reference Chaput2018) on the Ross Ice Shelf, we examine the spectral peak patterns in terms of temporal variability with respect to weather. As in Chaput and others (Reference Chaput2018), we discount instrument resonances as the root cause of these spectral peaks, given that the nodal seismographs are buried and have been tested for such behavior in a laboratory setting (i.e. Ringler and others, Reference Ringler, Anthony, Karplus, Holland and Wilson2018). Two broad observations further support this notion:

1) The seismic nodes were (shallowly) buried at depths of ~0.2 m, and remained so for the duration of recording. Temperatures remained significantly below freezing, precluding node exhumation due to local wind/instrument coupling or freeze/thaw effects (as observed, for instance, by Johnson and others (Reference Johnson, Vernon, Nakata and Ben-Zion2019)).

2) These patterns of resonances and their response to environmental forcing, summarized in Chaput and others (Reference Chaput2018), are observable across a variety of instrumentation configurations and sites, with both nodal and broadband instruments, and with and without surface flagging. Spectral peaks furthermore show spatial coherence on small to intermediate scales (e.g. Fig. 5) that preclude the possibility of individual instrument coupling or resonance effects, which would produce spatially random or spatially identical patterns.

Figure 6 shows three-component spectrogram segments compared with the gradient of the wind azimuth *θ* times the square of the wind speed *V* _{w} (to identify intervals of both strong winds and changing wind direction) and wind speed. The following observations are noted:

1) Resonance peaks are in general visible on all three components, but are stronger on the horizontals often by up to an order of magnitude in amplitude (although this is variable).

2) A resonance peak may show ‘splitting’, where peak doublets are observed that vary proportionately in frequency with time and show up to 15% peak frequency differences. In such cases, split peaks are typically recorded on both the vertical and horizontal components, possibly reflecting anisotropic structure (e.g. Fig. 6a).

3) Storm events with blowing snow (noted in field notes by days with very low visibility) were associated with frequency shifts by Chaput and others (Reference Chaput2018). Peaks are also observed to sometimes drift in frequency with changes in wind azimuth (Fig. 6b).

4) Peaks are sensitive to temperature anomalies, and will slowly decrease in both frequency content and amplitude with persistently high surface temperatures and may fluctuate daily even if temperatures approach freezing.

5) Peak patterns also occasionally show evidence of harmonic behavior following storm events (Fig. 6c).

The persistent mapping of energy onto all three components, the addition of splitting on the horizontals, and the increase of amplitude corresponding to higher wind speeds, suggest a *P* − *S* _{V} resonance excitation from a wind source. Although we cannot discount the possibility of *S* _{H} resonances being excited as well, they are not observed conclusively on the TIME array, as indicated by the mapping of peaks onto the vertical component even in cases of peak doublets. The increase in peak frequency content with the passing of storms is intriguing, and Chaput and others (Reference Chaput2018) hypothesized that multiple processes could account for this. Jumps in frequency content occur as winds are dying down and have time scales as short as a few hours. Two mechanisms were proposed, namely a change in surface coupling structure through deposition and stripping of surface material, and the rapid post storm formation of wind slab, a denser, more rigid layer of snow with low porosity several centimeters in thickness. Using 1.5D numerical simulations limited to flat layers with realistic firn velocity models, Chaput and others (Reference Chaput2018) were able to show that both mechanisms could produce an upward shift in frequency content, although no 2D or 3D effects were explored. Surface wave amplification effects were also noted by Dal Moro (Reference Dal Moro2020), who studied the effect of a thin, stiff layer overlaying a low velocity medium. In particular, higher mode Love waves were numerically shown to strongly emerge as harmonics with the presence of a stiff thin surface layer, a behavior potentially observed in Figure 6c, despite its mapping to the vertical component. A thorough numerical exploration of these various effects, particularly concerning deposition and temperature, can be found in Chaput and others (Reference Chaput2018).

As a further test of spatial variability, we roughly interpolate the lowest spectral peaks in Figure 5 according to the red lines, while making the assumption that the higher peaks noted by the black dotted ellipse are factor 2 harmonics of the peaks noted at nearest neighbors. These inferred peaks are then mapped to their respective stations and are compared with *H*/*V* and Rayleigh 0 particle motions (Fig. 7). Notably, a cluster of lower frequency peaks appears to map to the upper portion of the array in the same manner as lower frequency *H*/*V* peaks (Rayleigh 2) and prograde particle motions (Rayleigh 0), suggesting that surface wave effects may also be responsible for these amplifications. The first two effects, particle motions and *H*/*V* ratios for individual Rayleigh modes, can be numerically studied, and we do so in a 1D manner in the following section. Resonance peaks, however, have not been definitively constrained in terms of source/structure interaction. Thus, we subsequently explore several possibilities accounting for both SV plane polarization and the smooth spatial variation shown in Figure 5.

## 4. 1D numerical simulations

We seek to test the extent to which realistic near-surface velocity models can convincingly simulate the surface wave phenomenology observed here.

For 1D velocity and density profiles, pseudo-analytical methods for rapidly obtaining sensitivity kernels and other seismic properties are provided by several publicly available packages. We made use of RAYLEE (Haney and Tsai, Reference Haney and Tsai2015, Reference Haney and Tsai2017), a finite element model for surface wave simulation and inversion written in MATLAB©.

Firn densification processes on average tend to follow generally predictable velocity and density profiles controlled by their regional environments. For instance, compaction and ablation processes have been broadly studied in Antarctica, and firn profiles on average were modeled as having similar characteristics for RIS and WAIS Divide regions (van den Broeke, Reference van den Broeke2008). Local variations, however, can be severe and diverge significantly away from steady-state assumptions (e.g. Ligtenberg and others, Reference Ligtenberg, Helsen and van den Broeke2011; Reeh, Reference Reeh2008) owing to the formation of ice lens layers, hoar frost and other atmospherically driven structures. Given that a high-resolution seismic model is not yet available for this area, though such preliminary efforts are in progress as per Adam Booth (pers. comm. 2021), we proceed as follows. We begin with a smooth firn velocity model obtained from inversion of anthropogenic (camp noise) 4–18 Hz Rayleigh and Love waves recorded on the Ross Ice Shelf by a small-aperture (5 km nonuniform spacing) short-period seismic array (Diez and others, Reference Diez2016). Although likely not an exact average match to our site, perturbations of this model allow us to investigate variability in Rayleigh wave sensitivity kernels and eigenvalues in such media.

We investigate two main types of perturbations away from our starting firn model:

1) A series of smooth slowing perturbations for the firn profile, and

2) a series of low-pass filtered random Gaussian noise perturbations applied the initial firn profile.

We seek here to model several observations: the reversal of particle motions in Rayleigh 0; the specific pattern of *H*/*V* perturbations for all three modes; the likelihood that a planar deepening of firn structure across the array could account for spatial variations thereof.

The first set of perturbations seeks to gauge the impact of average firn structure variations on surface waves, given that the starting model from Diez and others (Reference Diez2016) may, on average, be slightly slower or faster than at our study site, and variable within the aperture of our array. Specifically, particle motion reversals observed in Figure 4 for the Rayleigh 0 mode should be sensitive to these forms of perturbations given the generally more broad nature of fundamental mode kernels. The second set of perturbations seeks to study the effect of severe shallow ‘switch backs’ (i.e. fast/slow/fast and similar transitions) in the top 10 s of meters of the firn, suspected here to be largely responsible for *H*/*V* amplifications of the second Rayleigh mode given its variation in peak frequency over a small spatial scale.

To examine particle motion directions, we adopted the convention of Tanimoto and Rivera (Reference Tanimoto and Rivera2005), whereby the radial eigenfunction is negative and, in the absence of a firn layer and in the majority of Earth media, fundamental Rayleigh particle motions are retrograde. For each perturbed profile, we computed the *ZH* ratio, defined by *ZH* = *U*(*R*)/(1.5 · *V*(*R*)), where *U*(*R*) and *V*(*R*) are the vertical and horizontal eigenfunctions (Takeuchi and Saito, Reference Takeuchi and Saito1972; Tanimoto and Rivera, Reference Tanimoto and Rivera2005), as a proxy for determining frequency ranges over which particle motion reversals can occur. A positive *ZH* implies prograde motions, whereas a negative *ZH* indicates retrograde motions. Similarly, the *H*/*V* ratio for a given Rayleigh mode can be calculated by the absolute value of the ratio of its radial to vertical eigenvalues, and so is closely related to the *ZH* ratio. It is primarily used as a metric sensitive to near-surface discontinuities, and higher mode Rayleigh waves here are obvious candidates for study.

For both types of perturbations, we generate *H*/*V* ratios, *ZH* ratios (indicative of particle motion directions) and depth sensitivity kernels for all modes, which provide additional information concerning potential mechanisms for spectral amplifications (Figs 5 and 6) and allow us to consider which portions of the firn are primarily sampled by surface waves.

Figure 8a shows the results of smoothly slowed firn profiles, starting with the RIS profile computed by Diez and others (Reference Diez2016). It should be noted here that ‘sped up’ versions of this profile did not result in zones of particle motion reversal for Rayleigh 0, and as such are not shown here. As velocity profiles move from faster to slower, a region of prograde Rayleigh 0 motion emerges that narrows and shifts to lower frequencies (Fig. 8a, bottom panels). These results would suggest at first glance that our firn profile at WAIS Divide is necessarily slower than the one on RIS given the emergence of a zone of prograde motions for models 3–10, perhaps due to different firn densification conditions (as per van den Broeke (Reference van den Broeke2008)). We can however note that the corresponding *H*/*V* ratios for this region are not equivalent to those observed in our data (i.e. a near-unity *H*/*V* for Rayleigh 0, a vertically dominated Rayleigh 1 and large values for Rayleigh 2). In short, broadly slowing down the firn profile tends to result in higher *H*/*V* ratios for Rayleigh 0, and this is not observed, though, as discussed previously this could also be due to contamination of Rayleigh 0 by body waves. We anticipate that more random and severe perturbations to these profiles exist, and thus subsequently perturb the profile in a different manner.

Figure 8b shows the results of ten models produced by adding low-pass filtered additive Gaussian perturbations to the initial RIS profile, with decreasing standard deviation as a function of depth, as is expected from densification processes. The amplitude of the perturbations was chosen arbitrarily to qualitatively mirror non-steady state densification models (Ligtenberg and others, Reference Ligtenberg, Helsen and van den Broeke2011). It is immediately clear that a large degree of variability, particularly in *H*/*V* ratios for Rayleigh 2, can be imparted into surface wave eigenvalues by severe shallow variations in structure without slowing the profile on average. *ZH* ratios for Rayleigh 0, although clearly variable over different perturbed models, never become positive for any of the models in Figure 8b, nor do *ZH* ratios for Rayleigh 1 become negative, suggesting that a slight slowing of the profile from the base model is indeed additionally required to account for our observations. Several of the random models (1,3,7 in Fig. 8b, Rayleigh 2 *H*/*V*) show that high values of *H*/*V* can be obtained in the 30–50 Hz range for Rayleigh 2, while maintaining near unity *H*/*V* for Rayleigh 0 and < 1 for Rayleigh 1. These realistic random models thus capture fine scale high-frequency effects related to the propagation of multi-mode Rayleigh waves, but also underscore the inherent difficulty in performing a unique structural inversion for the firn based on these data. A combination of the two perturbations is likely here (i.e. a slower average profile, and a set of random but severe perturbations), but a joint inversion of these quantities is beyond the scope of this paper.

We now turn our attention to the spatial variability inherent in our observations and in particular, the implication of a shift in *H*/*V* for Rayleigh 2 over the small aperture of our array. GPR studies at other Antarctic firn sites have revealed the presence of modulations in firn layering with scales of 100 s of m to 10 s of km (e.g. Sinisalo and others, Reference Sinisalo, Grinsted, Moore, Karkas and Pettersson2003; Arcone and others, Reference Arcone, Spikes and Hamilton2005a, Reference Arcone, Spikes and Hamilton2005b). Whether these observations are pervasive throughout Antarctic firn is unclear, but we anticipate that semi-planar modulations of firn layering over scales of 100 s to 1000s of meters are possible at our study site. We may not necessarily expect to be able to concurrently model a spatial reversal in *ZH* for Rayleigh 0 while also ensuring that its *H*/*V* remains near unity, given the inherent relationship between the two metrics (e.g. a particle motion reversal implies a zero crossing in *ZH* ratios, which results in either amplified or near-zero *H*/*V*) and the likely contamination from body waves as discussed in previous sections. If body wave contamination is indeed present, as we suspect, then a particle motion reversal as observed in Figure 4 is possible without a resulting *H*/*V* pole or zero. From the three random models that qualitatively produce Rayleigh 2 *H*/*V* in the right range, we choose model 7 (Fig. 8) as a base, and shift the pattern of random perturbations downward along initial firn profile for the RIS. Figure 9 shows the impact of such a structural downward shift on *H*/*V* and *ZH* for all three modes. The *H*/*V* peak at around 40 Hz for Rayleigh 2, for instance, can be moved downward significantly with a small linear shift in structure without affecting Rayleigh 0 or Rayleigh 1. The *ZH* ratio also shows a modest shift for Rayleigh 0 toward a zero crossing, suggesting that such a shift, if combined with a slightly slower firn model on average, could be just enough to cause a particle motion reversal. As such, there is strong evidence that a locally planar dipping firn structure on the order of a few meters is present under our array, deepening from bottom to top in Figure 7 (i.e. west to east), which would account for the 2π variations in *H*/*V* and particle motions noted in Figure 7b,c.

The final point to address here is more speculative and relates to the absolute amplifications observed in Figures 5 and 6. Unlike the observations we have modeled thus far, these spectral amplifications do not currently have a well-described source model despite initial attempts by Chaput and others (Reference Chaput2018). There are two modes of variability for these peaks: spatial (structural) and temporal (source and structure). Given that source studies would likely require high-resolution atmospheric and surface data (such as time lapse photogrammetry to study wind coupling locales), we focus here on structural variability and explore possible excitation mechanisms. Wave resonance mechanisms must account for the P-SV polarization of peaks and their similar sensitivity to spatial variation in firn structure as Rayleigh waves (Fig. 7). Given the novelty of these absolute spectral amplifications, we discuss them in terms of two possible mechanisms, and weigh them with respect to the evidence: *S* _{V} resonances and relative Rayleigh wave amplification. Given that we favor the latter, discussion concerning *S* _{V} resonances can be found in the online Supplementary material.

### 4.1. Surface wave resonances

Rayleigh and Love wave amplifications due to interactions with basin structures do not have closed form solutions beyond simple 1D approximations. Studies targeting surface wave amplification without resorting to strong motion seismograms generally follow a relative calculation, with one nearby station deployed over bedrock with a known velocity profile serving as a benchmark. Initial works by Tromp and Dahlen (Reference Tromp and Dahlen1992), expanded upon by Brissaud and others (Reference Brissaud, Bowden and Tsai2020), Bowden and Tsai (Reference Bowden and Tsai2017), and Feng and Ritzwoller (Reference Feng and Ritzwoller2017) have made use of broadly available ambient noise velocity maps of the continental US to predict and study amplifications in known basin structures while accounting for basin edge effects. Tromp and Dahlen (Reference Tromp and Dahlen1992) initially formulated relative amplification *A* _{n} with respect to a reference site $A_n^R$ as:

Here, *u* _{n}(0) and $u_n^R( 0)$ are the displacement eigenfunctions at the surface for the reference and basin sites, *n* = 1, 2, 3, and *U* is group velocity. *l* _{0} is an integral over eigenfunctions and density ρ and depends on the wave type. For Rayleigh:

and for Love:

We should note that this approach is useful, but is limited to 1D relative effects, and information pertaining to absolute amplification cannot be resolved (but can serve as a relative indicator of amplification in some cases). In our case, we consider the 20 models discussed in the previous section, with a ‘base model’ consisting of RIS profile. The first tenmodels consist of homogeneous slowing functions, and the last ten models consist of the RIS model with random additive Gaussian noise. Figure 10 shows the resulting relative amplifications for the first three Rayleigh modes observed in our data. Note that the frequency range that can be resolved is largely determined by the reference model (i.e. a much faster reference model may yield larger relative amplifications, but higher modes will not propagate at lower frequencies). Two primary effects are noted:

1) Large amplifications are possible with slower firn profiles (Fig. 10a), as expected, and Rayleigh 0 and 1 have similar amplifications (i.e.

*Vertical*~*Radial*). Rayleigh 2, however, exhibits disjointed radial and vertical amplifications, with a large radial dominance near 20 Hz and a large vertical dominance near 50 Hz. Large*H*/*V*ratios are indeed a feature of spectral peaks in Figure 5, which is intriguing, and although we do see radial-only amplifications at many frequencies, the models also predict vertically dominated modes at some frequencies, and this point is not observed in the data.2) Figure 10b shows examples of firn profiles containing the random noise perturbations. Amplifications relative to the initial profile are more modest, but display a wide range of behaviors, particularly in the higher modes. Interestingly, stronger discontinuities in themselves do not guarantee higher relative peaks, as extreme shallow gradients are present in the noisy models relative to the base, but we are limited in terms of interpretation by the relative aspect of these calculations.

Under this source interpretation, we may as before study the impact of shifting a random structure downward to see if a given amplification can also translate in frequency. As such, we choose an example model that has resulted in a peak in the − 20 to 30 Hz range (we arbitrarily choose model 8), which has a modest peak at ~ 25 Hz, and study the effect of a 2 m downward shift in the profile perturbations. As exemplified by Figure 10c, the resonance shifts slightly to lower frequencies. Although it agrees with *H*/*V* response, we do not consider this to be generally diagnostic, given the complexity of firn resonance phenomenology and the source-independent nature of these computations. Rayleigh wave amplifications can be obtained this way with a variety of plane polarizations, and it is thus clear that source excitation (as also inferred by Chaput and others (Reference Chaput2018)) plays a critical role in the generation of these features. As such, although this formulation is only a relative approximation that certainly does not capture the complexities of features shown in Figures 5 and 6, it may offer a glimpse into the structural aspect of firn that allows these resonances to occur.

We view Rayleigh wave amplification as an intuitive candidate here to describe much of the firn ambient seismic wavefield, given the likely *P* − *S* _{V} dominant excitation of wind noise, the relationship between wind and peak drift and the complex behavior of anthropogenically excited Rayleigh waves noted in this study. Although a more detailed nature of these amplifications is difficult to model and interpret, the *P* − *S* _{V} plane constraint on these resonances presents intriguing opportunities such as potential sensitivity to azimuthal anisotropy at single stations, which, in general, is not possible with ambient noise measurements without time domain reconstructions. Finally, the spatial mapping of Rayleigh 0 particle motions (and conversely, the approximate reverse mapping for Rayleigh 2) combined with a similar trend in *H*/*V* and resonance peak estimates (Fig. 7) all suggest sensitivity to similar structures, which reinforces the notion of surface wave amplification. In particular, *H*/*V* peak frequency, largely acknowledged to be a result of Rayleigh wave motion linearization, also maps to decreasing resonance peak frequency, suggesting a direct link between the two processes. We provide further discussion of alternate *S* _{V} explanations and the potential impact of multi-mode depth sensitivity kernels in the online Supplementary material.

## 5. Conclusions

We studied high-frequency (> 5 Hz) ambient seismic wavefield phenomena in Antarctic firn recorded by a 1 km circular array of 22 nodal seismometers deployed in the vicinity of WAIS Divide camp in January 2019. Three primary observations were used to detect shallow firn variability across the array:

1) Anthropogenic sources from the 5 km distant camp and wind forcing generated three distinct Rayleigh wave modes that showed clear modal hand offs. Particle motion directions for the fundamental mode (denoted Rayleigh 0) and conversely for the second higher mode (Rayleigh 2) showed a long period spatial transition from retrograde to prograde moving from west to east across the array (and from prograde to retrograde for Rayleigh 2).

2)

*H*/*V*for the highest detected mode (Rayleigh 2) showed the same spatial transition pattern, moving from higher to lower frequency content from West to East.3) Firn resonance modes identified by Chaput and others (Reference Chaput2018) and observed here show decreased frequency content from west to east and are spatially correlated with the other two Rayleigh wave metrics.

These observations are interpreted as indicating 10–1000 s of m-scale lateral velocity profile changes. Wavefield observations from the array are consistent with Antarctic firn supporting shallow basin wavefields within its characteristic low near-surface and high vertical gradient velocity structure. 1D numerical simulations using firn profiles constrained by seismic inversion at other Antarctic locales also confirms the feasibility of prograde fundamental Rayleigh wave motions in 5—15 Hz with realistic perturbations of the velocity profile. Modeled *H*/*V* ratios for the three isolated Rayleigh modes were suitably replicated for some firn models featuring random Gaussian noise perturbations, and the variability in Rayleigh 2 *H*/*V* and particle motions was accounted for by a small (~ 2 m) downward shift of the perturbation patterns. As such, we conclude that there is a small roughly planar dip in the structure of the firn from West to East, resulting in variable Rayleigh 2 *H*/*V* ratios, a reversal in particle motions for Rayleigh 0, and a similar decreasing pattern for absolute spectral amplifications.

Spectral amplifications in the *P* − *S* _{V} plane were proposed to result from multimode Rayleigh wave interactions with extreme low-velocity firn structures, though the exact generative mechanism remains an open question. The influence of storm activity observed here, which macroscopically and microscopically rearranges surface snow, along with spectral peak changes linked to wind direction, corroborates conclusions from Chaput and others (Reference Chaput2018) that resonances arise from interaction between very low velocity near-surface firn with lateral gradients modulated by layering and waves generated by wind coupling on snowforms (e.g. via dune and sastrugi spatial quasiperiodicity). A 3D array or the presence of rotational instruments would be especially valuable in further understanding the source mechanism of these resonances. Constraining their sensitivities could improve abilities to noninvasively monitor spatial and temporal firn state with simple instrumentation (e.g. small arrays or even single seismographs, augmented with GPR). Observed peak doublets in narrow-band resonances suggest sensitivity to azimuthal anisotropy, which may be exploited in future studies to constrain ice and firn fabric.

This work illustrates that the shallow Antarctic firn cryosphere provides pervasive manifestations of ambient seismic excitation with surface waves and other behaviors that are analogs to those encountered in basin seismology, but with unique phenomenologies. Future ambient noise imaging and temporal monitoring in firn media, particularly those that exploit assumptions of isotropic excitation, should exercise care in the analysis of these data. We expect Antarctica to host a range of firn structures under the influences of temperature, accumulation rates, precipitation and other factors, and it is unclear to what extent the few locales so far investigated may be broadly representative. More high-frequency analyses of ambient wavefields across Antarctica collected specifically or incidentally will no doubt shed light on these questions.

## Acknowledgements

The data and seismographic instrumentation are from the Thwaites Interdisciplinary Margin Evolution (TIME) project and the International Thwaites Glacier Collaboration and will be archived at the Incorporated Research Institutions for Seismology (IRIS) Data Management Center. We acknowledge support by the National Science Foundation (OPP-1739027 and OPP-1744852). The facilities of the IRIS Consortium are supported by the National Science Foundation under Cooperative Agreement EAR-1261681 and the DOE National Nuclear Security Administration. We thank Galen Kaip for his field support. We also thank the US Antarctic Program and the WAIS Divide camp and support staff during the field work. This work was also partially funded by University of Texas, El Paso startup funds (JC).

## Supplementary material

The supplementary material for this article can be found at https://doi.org/10.1017/jog.2021.135.