2.1. Instrumentation

The coordinated RIS (Mantle Structure and Dynamics of the Ross Sea from a Passive Seismic Deployment on the Ross Ice Shelf) and DRIS (Dynamic Response of the Ross Ice Shelf to Wave-Induced Vibrations) projects (Figs 1, S1) consisted of 34 polar-engineered broadband seismic stations provided by the Incorporated Research Institutions for Seismology (IRIS) Polar Programs. The RIS/DRIS stations were installed in late 2014 and recorded ~2 years of continuous data. Seismographs were deployed in two main transects: a 1100 km-long ice-front-parallel transect (W–E) and a 425 km long ice-front-perpendicular transect (N–S). The network consisted of (1) a shelf-spanning large aperture array (RS01–RS18) with a median spacing of 83 km and (2) a central medium aperture array (DR01–DR16) with stations spaced at 20–50 km. RS04 is located at the intersection of the two main transects. All stations were sited on floating ice, with the exceptions of RS08 and RS09 on Roosevelt Island, RS11–RS14 in Marie Byrd Land and RS17 on Steershead Ice Rise. .

Most RIS (RS) and DRIS (DR) stations utilized Nanometrics Trillium 120PH posthole sensors direct-buried at depths of 2–3 m below the snow surface; exceptions were RS09, RS11–RS14 and RS17, which were Nanometrics Trillium 120PA sensors installed on phenolic resin pads within shallow vaults. All DR stations and RS04 had a sampling rate of 200 Hz; all other RS stations had a sampling rate of 100 Hz. Stations used solar power during the Antarctic summer and single-use lithium thionyl batteries during the winter. Due to Iridium satellite modem power and bandwidth constraints, only state-of-health information was telemetered, and the network was therefore serviced in 2015 for intermediate data recovery and any other necessary servicing. The signal-to-noise analysis presented here incorporates data from the full deployment period of approximately November 2014 through November 2016. Stations RS10–RS14 remained deployed in Marie Byrd Land until early February 2017 due to logistical issues.

2.2. Teleseismic earthquake signals on the RIS

We perform a signal-to-noise ratio (SNR) analysis for teleseismic earthquake signals (M _w 5.5 or greater) observed by the RIS/DRIS stations at epicentral distances >30°, selecting P-wave, S-wave and surface wave arrivals and their immediate coda, with spectral bandpass filtering and signal times shown in Table 1. The epicentral distribution and the Gutenberg–Richter relation for all events are shown in Figure S2. Predicted arrival times for individual phases for all earthquakes are based on ak135 travel time curves (Kennett and others, Reference Kennett, Engdahl and Buland1995). Bandpass filtering was determined based on the peak periods of the median modified SNR curves.

Table 1. Parameters for teleseismic signals used in this study

$\bar {N}$ indicates the mean number of events across all stations for each season. Summer includes all events between 1 December and 31 March. Winter includes all events between 1 April and 30 November. P- and S-wave signal lengths were chosen based on a survey of observed coda durations and to avoid contamination by subsequent phases. ‘Mean Dist.’ indicates the mean arc distance and standard deviation across all stations and all events. The mean magnitude for all bands, across all stations and all events, was M _w = 6.0 ± 0.5.

To characterize teleseismic P-wave signals, we examine events at epicentral distances between 30° (slowness 0.08 s km⁻¹) and 95° (0.04 s km⁻¹) from each station. The teleseismic S-wave catalog is limited to events between 60° and 95° to prevent interference from surface wave arrivals near 5 km s⁻¹. Teleseismic P- and S-wave arrivals are separated by a minimum of ~300 s (at 30°) and therefore will not mutually interfere. The surface wave analysis uses wavetrains in the 4 to 2 km s⁻¹ arrival window, from earthquakes with epicentral distances between 45° and 100°, to prevent interference from body waves.

The ambient noise field of the RIS in the 0.4–30 s period band has seasonal variations of ~5–20 dB (Baker and others, Reference Baker2019). This variation is bimodal between the austral summer and winter and reflects the attenuation of wind-sea and ocean gravity waves in the Ross Sea by the formation of spatially continuous sea ice during the winter months. The absence of extensive, continuous sea ice determines the onset and termination dates of the ‘summer’ high-noise state analyzed in Baker and others (Reference Baker2019). We similarly subdivide our catalog into ‘summer’ earthquakes occurring during the approximate open-water interval between 1 December and 31 March, and ‘winter’ earthquakes occurring during the remainder of the year during which sea ice is broadly contiguous at the ice shelf front. Note that we henceforth refer to ‘winter’ and ‘summer’ noise conditions as indicating these annual sea ice-determined time periods, consistent with Baker and others (Reference Baker2019), and not the formal austral seasons.

2.3. Spectral characterization of teleseismic signals

The SNR ratios of individual events generally depend on a variety of source and propagation factors (e.g. moment, radiation pattern, depth, distance and mantle heterogeneities) which must be minimized to ensure that any observed temporospatial variations are solely the result of receiver-side features. We use event stacking and source normalization (described below) to approximate a globally-averaged teleseismic wavefield for comparison to station-specific noise fields. Additionally, we focus on the teleseismic excitation of the ice shelf beyond that of initial phase arrivals (in contrast to White-Gaynor and others (Reference White-Gaynor2019)).

For body wave arrivals, the phase arrival and its subsequent coda are extracted in fixed duration time segments (Table 1) regardless of epicentral distance, back azimuth and magnitude. S-wave signals with epicentral distances >~82° are referenced to the SKS phase. For surface wave trains, we analyze 4–2 km s⁻¹ surface waves, with a total signal length that is dependent on epicentral distance. Data segments for all three arrival types include 10 s of pre-arrival noise as a buffer against potential deviations from predicted arrival times and to allow tapering of processing artifacts without attenuation of the actual signal. For all segment types, the ‘noise’ time series is drawn from immediate pre-arrival data using an equal number of seconds as the ‘signal’ analysis. For example, for a P-wave arriving at t ₀ = 0 s, the noise data are comprised of T _n = [t ₀ − 110…t ₀] s, while the signal data spans T _s = [t ₀ − 10…t ₀ + 100] s. All data are de-meaned, de-trended, cosine-tapered by 5 s, decimated to 10 Hz, and rotated to (Z, R, T) (vertical, radial, transverse) coordinates using USGS NEIC Comcat parameters.

SNR spectra are calculated as the ratios of the signal and noise power spectral densities (PSD). PSDs are estimated using Welch's method (Welch, Reference Welch1967), with the number of (Hann-tapered) sub-segments determined by requiring 80% overlap and sub-segment lengths approximately equal to ten times the upper period limits indicated in Table 1. PSD estimates are subsequently smoothed by averaging over 1/8 octave bins in 1/16 octave increments.

Due to the geographic extent of the RIS/DRIS array and the nonuniform distribution of earthquake sources (Fig. 1), not all stations observed the same population of earthquakes given the event selection criteria. These disparities are especially exaggerated during our defined high-noise, low-sea-ice summer months (1 December–31 March), which include roughly a quarter of the number of events as the winter months (1 April–30 November). For the goals of this study, we require an idealized teleseismic source that is uniformly observed by all stations, such that interstation variations in SNR are predominantly controlled by receiver-side phenomena. To approximate such a source, we apply a correction analogous to the relative radiometric normalization used in remote sensing, but operating on the source rather than the receiver. We normalize individual earthquake SNRs to a theoretical earthquake model, with a scaling factor based on USGS seismic moment and losses from geometric spreading. For P- and S-waves, the corrected SNR′ is:

(1)$${SNR}' = {SNR} {M_{0{ref}}\over M_0} \bigg({D\over D_{{ref}}}\bigg)^{2} \ \comma\; $$

where SNR is the uncorrected signal-to-noise ratio, M ₀ and M _0ref are the event and reference seismic moments, respectively and D and D _ref are the depth-dependent, source-to-receiver ray path distances through a spherical ak135 Earth model. For surface waves, the correction is:

(2)$${SNR}' = {SNR} {M_{0{ref}}\over M_0} \bigg({\sin\lpar \Delta\rpar \over \sin\lpar \Delta_{{ref}}\rpar }\bigg)\ \comma\; $$

where Δ and Δ_ref are the great arc distances, and the sine function corrects for two-dimensional spreading (Stein and Wysession, Reference Stein and Wysession2009). We use a reference event with a magnitude of M _w = 6.0 (the approximate mean of all observed events), band-specific great arc distances as listed in Table 1, and a depth of 10 km.

We make no attempt to correct for source function beyond median stacking, nor do we eliminate overlapping events. For our S-wave analysis, the strength of the P-wave coda integrated into the S-wave noise estimate may be dependent on epicentral distance if the P-wave coda decays appreciably within the noise window. That is, S-wave noise estimates for more distant earthquakes may integrate a smaller percentage of the P-wave coda than those for less distant earthquakes. Based on our generally symmetric earthquake catalog (Figs 1, S2), we expect that median stacking will neutralize any spatial biasing relative to the receiver-side spatial variations. Surface wave noise estimates incorporate the entirely of the associated P- and S-wave codas regardless of epicentral distance and are therefore insensitive to this issue.

To evaluate seasonal variations in SNR, we calculate the seasonal median SNR from all appropriate summer and winter events, as defined above and by Baker and others (Reference Baker2019). We evaluate spatial variations across the RIS/DRIS array with respect to the mean of station-specific median seasonal PSDs for the arrival bandpass ranges listed in Table 1.

3.1. Elastic structure

The RIS spans an area of 487, 000 km² and floats on the relatively shallow continental shelf waters of the Ross Sea embayment with only sparse internal pinning points. RIS/DRIS station sites have ice thicknesses of 200–400 m (median 325 m; $\widetilde {h}$) and water column thicknesses of 100–700 m (median 464 m; $\widetilde {H}$) in BEDMAP2 (Fretwell and others, Reference Fretwell2013).

Laterally varying structures in ice shelves are introduced by the convergence of source glaciers near the grounding line and by the subsequent advection of shelf ice toward the RIS calving front, at up to 1 km a⁻¹ for the RIS (Mouginot and others, Reference Mouginot, Rignot and Scheuchl2019). Internal structures include suture zones between tributary glaciers that persist to the terminus of the RIS ice front. Also present are subglacial and surface crevasses and rifts that are generally parallel with the calving front, reflecting a broadly tensional stress state in the seaward flow direction. This tensional stress field is magnified near the free-floating terminus (LeDoux and others, Reference LeDoux, Hulbe, Forbes, Scambos and Alley2017) and is cyclically influenced by tidal tilt (e.g. Olinger and others, Reference Olinger2019). Significant vertical structure includes a meteoritic snow-firn layer that transitions to glacial ice over tens of meters (e.g. Diez and others, Reference Diez2016), and for some shelves, a seasonally modulated basal freeze layer. The sub-shelf seafloor may be comprised of up to several kilometers of low velocity lithified sediment overlying a high velocity basement (e.g. Beaudoin and others, Reference Beaudoin, ten Brink and Stern1992).

We assume bulk elastic properties for the RIS system as listed in Table 2. We also assume that the RIS is laterally and vertically homogeneous and isotropic, with ice/water and water/seafloor interfaces that vary smoothly at scales much longer than the ice thickness. We justify this simplification by noting that the RIS/DRIS stations were deliberately sited in regions of solid glacial ice several kilometers or more distant from crevassed or rifted areas as a matter of safety. For teleseismic body wave arrivals (i.e. Table 1) with near-normal angles of incidence, the maximum ice shelf basal piercing distance would be 230 m for a P-wave with a ray parameter of 0.08 s km⁻¹, assuming a 700 m thick shelf with a P-wave velocity of 3.87 km s⁻¹. For vertical features, our simplification is justified by the resolution limit of the teleseismic waves of interest. Assuming a quarter-wavelength limit, the highest resolving wave would be an SV-wave generated by a teleseismic P-wave at the ice/water interface, with a period of 0.5 s and a seismic velocity of 2.0 km s⁻¹, and a minimum resolution of 250 m. For contrast, the low-velocity firn layer becomes important only in the uppermost 60 m (Kirchner and Bentley, Reference Kirchner and Bentley1979; Beaudoin and others, Reference Beaudoin, ten Brink and Stern1992; Diez and others, Reference Diez2016).

Table 2. Elastic parameters for the RIS used in this study

Firn values are the geometric mean of empirical results published in the listed references.

^aKirchner and Bentley (Reference Kirchner and Bentley1979).

^bFofonoff and Millard (Reference Fofonoff and Millard1983).

^cBeaudoin and others (Reference Beaudoin, ten Brink and Stern1992).

^dKing and Jarvis (Reference King and Jarvis2007).

^eGriggs and Bamber (Reference Griggs and Bamber2011).

^fFretwell and others (Reference Fretwell2013).

^gDiez and others (Reference Diez2016).

3.2. Intralayer resonances

High seismic impedance contrasts in the vertical structure of an ice shelf system – specifically at the ice/water, water/seafloor and sediment/basement interfaces – create strong reverberatory wavefields by multiply reflecting incident body waves. A similar but less severe effect is common in geologic basins, where low velocity sediment overlies high velocity basement rock. The amplitude of this effect is maximized when the wavefield constructively interferes with itself. For any individual layer, the resonance periods, P _R are:

(3)$$P_{\rm R} = {2Z\eta\over n + m} \qquad {\rm for}\ n = 1\comma\; 2\comma\; 3\comma\; \ldots \comma\; $$

where Z is the layer thickness, η is the vertical slowness of the incident primary wavefield, the integers n are harmonic mode orders and m accounts for phase reversals at reflective interfaces (modified from Press and Ewing, Reference Press and Ewing1951; Crary, Reference Crary1954). For resonances spanning multiple layers, the period is simply the series summation of all layers. For the ice shelf system, strong resonance wavefields may be generated for S-waves within the ice shelf and P-waves within the water column.

The layered structure of an ice shelf creates a highly efficient waveguide for plane-media polarizations of shear waves. SH-waves within an ideal ice shelf with perfectly horizontal boundaries undergo lossless, in-phase reflections (m = 0) at the free surface and the ice/water interface; at any resonance period, P _R, an SH-polarized shear wavefield will propagate laterally along the ice shelf as a shear horizontal plate wave (Press and Ewing, Reference Press and Ewing1951; Rose, Reference Rose1999). SV-waves will experience lossless, out-of-phase reflections (m = 1) from either boundary only at the critical angle $\theta _{\rm c} = \sin ^{-1}\lpar \beta _{\rm i} \alpha _{\rm i}^{-1}\rpar$; a critically reflected SV-wavefield (m = 1) at any P _R will propagate laterally along the ice shelf as a Crary wave (Crary, Reference Crary1954). SV-waves incident on the ice/water interface at all angles other than θ _c will lose energy via SV-to-P conversions within the ice and water (Fig. S3).

For the structure specified in Table 2, θ _c = 29.6°, equivalent to a ray parameter of p = 0.263 s km⁻¹. Therefore, we do not expect to observe excitation of Crary waves by teleseismic body waves, for which p = 0.04–0.12 s km ⁻¹. However, the steeply incident ray parameters typical of teleseismic P-waves result in non-negligible P-to-SV conversion coefficients (10–20%) for both the free surface and the ice/water interface, with strong SV reflection coefficients (>80%), and no phase shift (m = 0) (Fig. S3). Furthermore, for the ice shelf values listed in Table 2, and a ray parameter of p = 0.06 s km⁻¹, Eqn (3) yields P _R = 0.33 s, which is generally within the spectral content of the P-wave teleseismic signal. We therefore expect to observe significant SV-resonant energy associated with these arrivals. We will refer to these as Crary resonances to differentiate them from true, lossless Crary waves.

SH-wave energy may potentially be observed coincident with teleseismic P-wave arrivals as a result of scattering from sloping interfaces or due to ice anisotropy. We do not expect to observe Crary resonances associated with teleseismic S-waves, due in part to their longer periods being incompatible with the RIS resonant periods, and for additional reasons to be discussed later in this section.

Extending the previous resonance analysis to the water column, we expect teleseismic P-wave codas to excite acoustic resonances (m = 0) with periods near P _R = 0.61 s. This resonance is markedly less efficient than the Crary resonance and will leak P- and SV-wave energy into the ice and seafloor. In Section 4.1 we show that the spectral signatures of these Crary and acoustic resonances may be exploited to accurately estimate the thicknesses of the RIS and the sub-shelf water column.

3.3. Coupling of seismic, acoustic and gravity waves

For elastic waves within an open water column of finite depth H (i.e. bounded by a free surface and a solid elastic Earth), there exists an acoustic cutoff period, $P_{\rm c} = 4H \alpha _{\rm w}^{-1}$, above which any vertical component of displacement ceases to propagate and becomes evanescent (i.e. decays exponentially in the vertical direction) (Ewing and others, Reference Ewing, Jardetzky and Press1957). For oscillations at periods greater than P _c, gravity becomes increasingly relevant to the vertical restoring force, and completely dominates for periods greater than the Brunt-Väisälä period (Apel, Reference Apel1987). For the period band bounded by the acoustic cutoff period and the Brunt-Väisälä period, the propagation of oscillatory energy through the water column lies in the acoustic-gravity wave domain. Traer and Gerstoft (Reference Traer and Gerstoft2014) and Ardhuin and Herbers (Reference Ardhuin and Herbers2013) provide derivations of this wave type as an extension of the nonlinear wave–wave interaction model first proposed by Hasselmann (Reference Hasselmann1966). From the perspective of continuum mechanics, this acoustic-gravity regime may be understood as a transition between compressible and incompressible fluid behavior with increasing period. Useful analytical representations for the response of an ocean above vertically displaced ocean floor are given by Yamamoto (Reference Yamamoto1982). Literature on acoustic-gravity waves is extensive but has focused on wave mechanics confined to the water layer. Recent studies have begun to address acoustic-gravity waves in the presence of an ice layer, but have thus far been restricted to inelastic ice caps (Kadri, Reference Kadri2016) or thin (h ≪ H) elastic sea ice (Abdolali and others, Reference Abdolali, Kadri, Parsons and Kirby2018) and thus do not address coupling between acoustic-gravity and seismic waves.

To obtain an estimate for the cutoff period for the combined ice shelf and water column at the RIS, we use the derivation for the similar case of a three-layered liquid half space (Ewing and others, Reference Ewing, Jardetzky and Press1957). We justify this simplification by noting that teleseismic P-waves propagate within the water column at an angle of incidence of 3.3–6.6° (0.04–0.08 s km⁻¹) and therefore lose <5% energy at any interface via conversion to SV-waves (Figs S3, S4). This approximation does not account for the flexural rigidity of the ice shelf nor the overlying firn layer. The cutoff period, P _c, in this case may be determined using:

(4)$$\tan{m_{\rm i}h} = {m_{\rm i}\over m_{\rm w} \delta_1} \bigg({ \delta_2 m_{\rm s} \tan{m_{\rm w} H} - m_{\rm w} \over \delta_2 m_{\rm s} + m_{\rm w} \tan{m_{\rm w}H} } \bigg)\ \comma\; $$

where m _i,w,s are vertical wave numbers within the ice, water and sediment layers, respectively, and δ ₁ = ρ _i/ρ _w, δ ₂ = ρ _w/ρ _s (Ewing and others, Reference Ewing, Jardetzky and Press1957). For the values listed in Table 2, P _c ≈ 2.25 s and varies within ±0.6 s for the range of sub-station ice and water layer thicknesses.

Short-period teleseismic arrivals associated with the P-wave coda (0.5–2.0 s) are therefore predicted to be observable on all three channels at floating stations. However, longer-period (e.g. 10–15 s) teleseismic compressional energy (e.g. S-to-P converted phases originating at the sediment/basement interface) should be observable solely on the vertical channel, reflecting the water column's predominately incompressible response to longer-period vertical seafloor displacement (e.g. Ewing and others, Reference Ewing, Jardetzky and Press1957). Observations of Rayleigh waves (>15 s) are similarly expected to be vertically dominant (as previously noted by Okal and MacAyeal (Reference Okal and MacAyeal2006)). Love wave arrivals should not be observable by floating stations above a planar seafloor due to the zero traction horizontal boundary condition.

3.4. Coupling of seismic and plate waves

A floating ice layer is an elastic plate system and therefore supports a variety of elastic modes that are not observed in solid-Earth or ocean-bottom seismology. On an ice shelf, plate modes are expected to be excited by tractions imposed upon the ice by ocean gravity waves at the ice front and within the sub-shelf water cavity. At ultra-long periods (>3 h), ocean waves may generate normal modes across the entirety of the RIS, which may, in turn, couple into acoustic-gravity waves within the atmosphere (Godin and Zabotin, Reference Godin and Zabotin2016). At periods of 50–100 s, infragravity and tsunami waves originating up to thousands of kilometers away generate flexural-gravity waves (i.e. buoyancy-coupled asymmetric mode (A₀) Lamb waves) which propagate several hundred kilometers into the RIS interior (Bromirski and others, Reference Bromirski2017). In the 1–50 s period band, ocean swell impacts at the ice front generate zeroth-order, symmetric mode (S₀) Lamb waves, in addition to flexural-gravity modes (Chen and others, Reference Chen2018; Aster and others, Reference Aster2019). On the RIS, flexural-gravity wave motion dominates the ambient wavefield at distances <120 km from the ice front. At greater distances, extensional wave motion (i.e. S₀ Lamb) dominates and is observable at least 450 km from the ice front (Chen and others, Reference Chen2018; Baker and others, Reference Baker2019).

Of interest to the present study are the symmetric mode Lamb waves (Lamb, Reference Lamb1917). Particle motion for these waves is retrograde elliptical in the vertical/radial plane, with a high eccentricity oriented parallel to the radial axis (Fig. 2). Lamb waves may be generated within a plate by compressional forces applied normal to the plate end face, or by shear forces applied in traction parallel to the basal plate surface. The theory and application of these methods have been rigorously documented in a number of publications on ultrasonics, where Lamb waves are often used in nondestructive testing (e.g. Viktorov, Reference Viktorov1967; Rose, Reference Rose1999). For this study, we concern ourselves with Lamb wave generation via normal end face tractions, such that:

(5)$$E_{\rm L} \propto A^{2} \cos^{2}\phi \ \comma\; $$

where E _L is the energy of the resulting Lamb wave, A is the amplitude of the source signal and ϕ is the angle between the source traction and the normal to the plate face (Rose, Reference Rose1999). For simplicity, we have restricted the source signal to oscillating within the same plane as the plate surface (i.e. horizontal). For the generation of ocean-coupled Lamb waves, the RIS is oriented such that ocean gravity waves propagating across the Ross Sea generally impact the ice shelf face at oblique-to-normal angles of incidence, maximizing the transfer of horizontal impulse between the incoming waves and the shelf.

The grounding lines of ice shelves are also plate boundaries at which Lamb wave energy may be generated. A relevant situation for this study is the case of teleseismic S-waves arriving near the grounding lines at nearly vertical angles of incidence (Fig. 2). Both SV- and SH-waves have the potential to generate Lamb wave energy (E _L) according to the general relationships:

(6)$$E_{\rm L} \propto {A_{\rm sv}}^{2} \cos^{2}\varphi \comma$$

(7)$$E_{\rm L} \propto {A_{\rm sh}}^{2} \sin^{2}\varphi \comma$$

(8)$$\varphi = 90^{\circ}-\theta_{\rm gl}-\theta_{\rm eq} \comma$$

where A _sv and A _sh are the amplitudes of the SV- and SH-waves, respectively, θ _gl is the strike of the grounding line (with the RIS defined as down-dip using the right-hand-rule convention) and θ _eq is the back azimuth to the source earthquake.