Grounding zone subglacial properties from calibrated active-source seismic methods

The grounding zone of Whillans Ice Stream, West Antarctica, exhibits an abrupt transition in basal properties from the grounded ice to the ocean cavity over distances of less than 0.5–1 km. Active-source seismic methods reveal the downglacier-most grounded portion of the ice stream is underlain by a relatively stiff substrate (relatively high shear wave velocities of 1100±430 m s⁻¹) compared to the deformable till found elsewhere beneath the ice stream. Changes in basal reflectivity in our study area cannot be explained by the stage of the tide. Several kilometres upstream of the grounding zone, layers of subglacial water are detected, as are regions that appear to be water layers but are less than the thickness resolvable by our technique. The presence of stiff subglacial sediment and thin water layers upstream of the grounding zone supports previous studies that have proposed the dewatering of sediment within the grounding zone and the trapping of subglacial water upstream of the ocean cavity. The setting enables calibration of our methodology using returns from the floating ice shelf. This allows a comparison of different techniques used to estimate the sizes of the seismic sources, a constraint essential for the accurate recovery of subglacial properties. We find a strong correlation (coefficient of determination=0.46) between our calibrated method and a commonly used multiple-bounce method, but our results also highlight the incomplete knowledge of other factors affecting the amplitude of seismic sources and reflections in the cryosphere.

1 Introduction

Grounding zones mark the transition from grounded to floating ice, standing sentinel over much of the contribution of glaciers and ice sheets to sea level. Within the grounding zone the location where the ice sheet ceases contact with the bed (the grounding line) is primarily determined by ice thickness, bed elevation, and the stage of the tide. In the Antarctic, tidally induced migration of the grounding line within the grounding zone varies from near zero in the case of abrupt changes in bed elevation and/or ice thickness to up to 10 km in the case of gently sloping ice plains (Brunt et al., 2011; Dawson and Bamber, 2020). Along with grounding line migration, tides correlate with ice velocity changes upstream and downstream of the grounding zone. Observations include daily velocity variability on Bindschadler Ice Stream (Anandakrishnan et al., 2003); twice daily stick–slip displacement on Whillans Ice Plain (Bindschadler et al., 2003; Winberry et al., 2009; Walter et al., 2011); daily and spring–neap velocity variability on the Ronne–Filchner Ice Shelf, Ross Ice Shelf, and Byrd Glacier (Rosier and Gudmundsson, 2020; Brunt et al., 2011; Marsh et al., 2013); and spring–neap tidal velocity variability on Rutford Ice Stream (Gudmundsson, 2007). Observed velocity variability has generally been attributed to tidal changes in the force balance interacting with the underlying till rheology (Bindschadler et al., 2003; Gudmundsson, 2007; Winberry et al., 2009). Subsequent studies have attributed Rutford Ice Stream's spring–neap velocity variability to changes in subglacial pore water pressure (Rosier et al., 2015), while on Rutford and elsewhere others have pointed to contact with ice shelf pinning points and at the grounding zone as the causes of observed velocity changes (Robel et al., 2017; Minchew et al., 2017; Rosier and Gudmundsson, 2020).

Early efforts to model tidal deflection of ice shelves primarily addressed vertical displacement and the associated development of strand cracks and basal crevassing at the grounding zone (Holdsworth, 1969, 1977). These models, termed stiff-bed fixed-grounding-line models by Sayag and Worster (2013), do not allow the grounding line to migrate, nor do they allow the underlying bed to deform. Despite inconsistencies in the retrieved elastic properties, subsequent applications of these models have successfully reproduced surface displacement (e.g. Vaughan, 1995; Schmeltz et al., 2002) with models accounting for basal crevassing (Rosier et al., 2017) and treating the ice as a viscoelastic material (Wild et al., 2017) shown to be more consistent with observations. The importance of grounding line migration for ice dynamics and the sensitivity of ice flow to tidal forcing has prompted renewed examination of the effect of tides on grounding line migration distances and subglacial conditions both within and upstream of the grounding zone. Sayag and Worster (2011) combined laboratory observations and an elastic sheet model in an analysis that allowed the grounding line to migrate over an elastic bed. Their approach was extended to the implications for subglacial water pressure (Sayag and Worster, 2013), showing pressure gradients alternating direction upstream of the grounding zone forming migrating barriers to subglacial water flow. Walker et al. (2013) used a fixed-grounding-line model with no vertical displacement at the grounding line and a viscoelastic ice sheet–shelf overlying an elastic bed. This approach resulted in alternating pressure gradients that may act to draw water from the ocean cavity at low tide and force it upstream at high tide. Tsai and Gudmundsson (2015) applied a novel elastic fracture approach to grounding line migration, which resulted in migration distances significantly different to elastic beam or hydrostatic approaches. Notably, Tsai and Gudmundsson (2015) demonstrated an asymmetry in grounding line migration whereby for a constant surface slope and a constant coastward bed slope, the grounding line migrates upstream as the tide rises from mean sea level much further than it propagates downstream when the tide falls from mean sea level. The subglacial system can also filter forcings, leading to velocity changes that occur at unexpected frequencies (e.g. Rosier et al., 2015). Robel et al. (2017) attributed such behaviour to the viscoelastic response of the ice shelf as it responds to changes in contact and buttressing at the grounding zone and pinning points. Alternatively, Warburton et al. (2020) coupled processes of upstream fluid flow beneath an elastic sheet and drainage through porous till and showed ice streams and ice shelves can respond at a range of frequencies. They also suggested ocean water may be retained in the subglacial system depending on the porosity of the till.

Grounding zones have been directly observed in only a few locations around Antarctica. Beneath Langovde Glacier in East Antarctica Sugiyama et al. (2014) reported a substrate of fine sediment with decimetre-scale dropstones, along with an incursion of sea water far beyond the previously mapped grounding line. In the ocean cavity proximal to the grounding line of Mackay Glacier, Powell et al. (1996) imaged a diverse range of glaciomarine lithologies, ranging from soft till to bedrock and dropstone boulders. Approximately 3 km downstream from Whillans Ice Stream's grounding zone, the WISSARD programme (Fricker et al., 2010) observed an ice shelf melt-out deposit with a mixture of soft mud and rock clasts (Scherer et al., 2015). Begeman et al. (2018) reported oceanographic and geophysical observations from the WISSARD borehole where they found a highly stratified water column with basal melt rates of less than 0.1 m a⁻¹. To further investigate the basal properties beneath Whillans Ice Stream's grounding zone, we here revisit the active-source seismic data reported by Horgan et al. (2013b) and apply and extend amplitude analysis methods previously used in studies addressing the basal boundary of glaciers and ice sheets (e.g. Anandakrishnan, 2003b; Smith, 2007; Holland and Anandakrishnan, 2009; Brisbourne et al., 2017; Zechmann et al., 2018; Muto et al., 2019). These methods require source amplitude and path effects to be estimated, which is often challenging due to variability in source and receiver coupling and strong vertical gradients in density and seismic velocity in the firn. Acquiring data over the ocean cavity allows calibration of these methods due to the presence of a known ice–water reflection interface. This allows us to use and expand on the methods of Holland and Anandakrishnan (2009) (hereafter referred to as H&A2009). H&A2009 reviewed active-source seismic methods for the recovery of subglacial properties, outlined best practices for reducing uncertainties, and presented new strategies for source size determination. Our application and extension of their methods enable a robust estimate of elastic properties beneath the ice at a relatively high spatial resolution. Our profile data cover approximately 50 line kilometres. The nominal horizontal resolution of our method is 240 m (based on the spatial footprint of a 100 Hz wave in a 3860 m s⁻¹ medium at a depth of 760 m), and we are able to image the top and bottom of a water layer ≥3.6 m thick ($\mathit{\lambda }/\mathrm{4}$, where λ denotes wavelength, of a 100 Hz wave in a 1440 m s⁻¹ medium). In theory, water layers down to $\mathit{\lambda }/\mathrm{32}$ (0.45 m) can be imaged; however amplitudes from these layers may not be representative of their elastic properties (e.g. Booth et al., 2012). To explore the relationship between the tidal stage and our results, we also present the timing and tidal stage of our experiment and Global Navigation Satellite System (GNSS) repeat transects along two profiles crossing the grounding zone.

2 Data and methods

We performed amplitude analysis of data from four transects that cross the grounding zone of Whillans Ice Stream (Fig. 1). These data were acquired in the austral summer of 2011/12. Acquisition was composed of an explosive seismic source detonated at approximately 27 m depth, with charge sizes of 0.4 (line 1) and 0.8 (lines 2 and 4) and 0.85 kg (line 3) at a nominal shot spacing of 240 m. Each of line 3's 0.85 kg charges was composed of one 0.4 kg charge and three narrower 0.15 kg charges. All other charges were composed of equal-diameter 0.4 kg charges. The time between burial and detonation varied but always exceeded 24 h. Geophones were buried approximately 0.5 m beneath the snow surface at 20 m spacings and consisted of alternating single-string 40 Hz geophones (even channels) and five-element 40 Hz georods (odd channels; Voigt et al., 2013). Acquisition used an asymmetric split spread with near and far shot–receiver offsets of 10 and 1430 m. Seismic imaging and grounding zone determination at Whillans Ice Stream is presented in Horgan et al. (2013b).

Figure 1 Location map showing the seismic profiles (lines 1–4) crossing the grounding zone of Whillans Ice Stream. Radio echo sounding (RES) basal reflectivity from Christianson et al. (2016). Seismic bed reflectivity (R_b) from this study. Background imagery from MODIS MOA (Haran et al., 2005) and grounding line from Bindschadler et al. (2011). Polar stereographic projection (metres) with a true scale at 71^∘ south.

Following H&A2009, the amplitudes reflected off of the base of the ice and recorded at our geophones (A_i, where i denotes the receiver index) are related to our source amplitude (A₀) by

$$\begin{array}{}\text{(1)}& A_i=A_{\mathrm{0}}{\mathit{\gamma }}_{i}R\left(\mathit{\theta }\right)e^{-\mathit{\alpha }{s}_i}\text{(Eq. 1, H\&A2009)},\end{array}$$

where R(θ) denotes the angle-dependent (θ) reflection coefficient at the base of the ice described by the Zoeppritz equations (e.g. Aki and Richards, 1980). During travel along the path length (s_i) from the source to the receiver, amplitudes are modified by path effects (γ_i) and attenuation (α), all of which are discussed below. Both A₀ and γ_i are amplitudes relative to a reference range (typically d₀=1 m; Holland and Anandakrishnan, 2009; Shearer, 2009).

2.1 Seismic velocity model

Tracing seismic ray paths between the source and receivers requires knowledge of the firn and ice column's seismic velocity. To achieve this we estimated a one-dimensional (1D) velocity model using shallow-seismic-refraction techniques. During shallow refraction surveying, a hammer source was recorded at 0.5 m horizontal intervals with near and far offsets of 0.5 and 579 m. A velocity model (Fig. 2) was then calculated using first-break arrival times and the τ–p (intercept time–slowness) method (e.g. Shearer, 2009), which assumes that the velocity monotonically increases with depth. This method estimated a velocity of 3840 m s⁻¹ at 80 m depth. Below this depth our velocity model consists of an extrapolation to a V_p corresponding to −20^∘C (3860 m s⁻¹; Kohnen, 1974), which is kept constant to the ice base. Kohnen (1974) demonstrated a decrease in V_p of 2.3 m s⁻¹ per degree Celsius decrease in temperature, so our velocity is fairly insensitive to our choice of temperature. Also implicit in our use of a 1D velocity model is an assumption that seismic velocity does not vary laterally within the survey area.

Figure 2 One-dimensional compressional wave velocity profile estimated using the τ–p method.

2.2 Amplitude picking

Amplitudes were picked on frequency-filtered and amplitude-scaled shot records guided by profiles stacked by common-depth-point. On every shot record we attempted to digitise the direct arrival, primary bed return, and first long-path multiple of the bed return (Fig. 3). The low impedance contrast at the ice–bed interface meant the long-path multiple could not be reliably picked in the grounded part of the profiles. Amplitude picking selected the zero crossing preceding the side lobe of the wavelet. Amplitude extraction was then performed on shot records with only bandpass filtering applied. Amplitudes were extracted within the wavelet encompassing the first side lobe, the central lobe, and the next side lobe. Within this wavelet, peak positive, peak negative, and root-mean-squared (rms) amplitudes were extracted. We avoided picking bed returns where direct arrival energy interferes with the bed wavelet. Our data are from ice thicknesses of approximately 730–790 m, and direct arrivals interfere with the reflection from the base of the ice beyond offsets of approximately 700 m. While the channels with five-element georods showed better signal-to-noise ratios for imaging, we here present an analysis of the single-string geophones as their amplitudes exhibit less channel-to-channel variability, the cause of which we attribute to more variability in coupling when burying the georods. Our analysis also uses the rms amplitudes, with the positive and negative peaks used to define polarity. We tested the use of peak amplitudes and fixed-wavelet-length approaches and found both resulted in a greater distribution of source sizes and less robust estimates of basal reflectivity.

Figure 3 (a) Example shot record from floating portion of line 2 (kilometre 4.8–6.7). (a) Inset shows schematic travel paths for direct (red), primary (purple), and multiple (red) rays. Right-hand panels show wavelets and picks for the direct arrival (b), primary return (c), and multiple return (d).

2.3 Path effects

Path effects (γ_i) modify the source amplitude during its propagation to the receiver. We calculated the total path effects as

$$\begin{array}{}\text{(2)}& {\mathit{\gamma }}_i={\displaystyle \frac{\cos {\mathit{\theta }}_i}{s_i}}\sqrt{{\displaystyle \frac{z_{\mathrm{0}}}{z_{\mathrm{1}}}}},\end{array}$$

where θ_i denotes the angle between the incoming ray and normal incidence, z₀ and z₁ denote the acoustic impedance at the source and receiver respectively, and s_i denotes the path length travelled between the source and receiver. Equation (2) therefore accounts for the angle at which the incoming ray arrives at the vertical-component receivers (cos θ_i), amplitude scaling due to the different acoustic impedance at the source and receiver ($\sqrt{\frac{z_{\mathrm{0}}}{z_{\mathrm{1}}}}$, e.g. Shearer, 2009), and geometric spreading along the ray path $(\mathrm{1}/s_i)$. We estimated all near-field effects using the 1D velocity model (Fig. 2) and the density–compressional-wave-velocity relationship of Kohnen and Bentley (1973). The high vertical gradients in density and velocity in polar firn lead to a cos θ_i correction≈1, as θ_i≈0, and a significant $\sqrt{\frac{z_{\mathrm{0}}}{z_{\mathrm{1}}}}$ correction $(\sim \sqrt{\mathrm{10}})$ due to the different source and receiver burial depths.

H&A2009 noted that placing receivers at a free surface results in a further doubling of recorded amplitudes for normal incidence returns. We tested including free surface amplification but did not apply it to the analysis presented here due to the burial of our receivers, although the shallow burial depth of 0.5 m may justify its inclusion. If included, this additional amplification would have resulted in a halving of the source sizes for two of our methods (the multiple-bounce method and the known reflector method; Sects. 2.4.1 and 2.4.3, Table 1). Including free surface amplification would have had a small effect (<17 %) on the direct-path method source size median values (Sect. 2.4.2, Table 1). Regardless of whether free surface amplification is included or excluded, our choice of preferred method for estimating A₀ would not change. The recovered bed properties also would not change as the same path effects used to calculate source size are later used to estimate bed properties.

Table 1 Source size (A₀) estimates.

2.4 Source size and attenuation

Source size (A₀) is often estimated using the ratio of the primary bed return amplitude (A_i) and the long-path multiple amplitude (A_m,i) (e.g. Röthlisberger, 1972; Smith, 1997; Peters et al., 2008; Brisbourne et al., 2017; Zechmann et al., 2018). This approach, termed the multiple-bounce method by H&A2009, removes the need for an independent estimate of attenuation. However, low impedance contrast at the bed, low signal-to-noise ratios, or closely spaced subglacial reflectors, can all complicate the multiple-bounce method of determining source amplitude. Here we explore this and other methods for determining the source amplitude because more-accurate source amplitude estimates will enable improved investigation of the basal properties resolved by seismic surveys. These methods fall into three categories: (1) multiple-bounce methods, (2) direct arrival methods, and (3) known reflector methods. We present the results for each of the four profiles individually as three different source sizes and configurations were used.

2.4.1 Multiple-bounce methods

Our multiple-bounce methods used the primary–multiple amplitude ratio to estimate A₀ and followed H&A2009. The first method requires near-normal incidence returns but does not require knowledge of attenuation (α):

$$\begin{array}{}\text{(3)}& A_{\mathrm{0},i}={\displaystyle \frac{A_i^{\mathrm{2}}}{A_{\mathrm{m},i}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\gamma }}_i}}\text{(Eq. 6, H\&A2009)},\end{array}$$

and the second method requires close-to-normal incidence returns and an estimate of attenuation:

$$\begin{array}{}\text{(4)}& A_{\mathrm{0},i}={\displaystyle \frac{A_i^{\mathrm{2}}}{A_{\mathrm{m},i}}}{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{m},i}}{{\mathit{\gamma }}_i^{\mathrm{2}}}}{e}^{\mathit{\alpha }(\mathrm{2}{d}_i-d_{\mathrm{m},i}})\text{(Eq. 7, H\&A2009)},\end{array}$$

where d_i and d_m,i and γ_i and γ_m,i denote the path length and path amplitude factor (Eq. 2) for the primary and multiple bed returns respectively. A₀ is then calculated as the average A_0,i for each shot. Equations (3) and (4) give near-identical A₀ estimates with root-mean-squared differences ≈0.1 %. Henceforth for the amplitude ratio method we report only the results from Eq. (4) with an angle cut-off of <10^∘ and assuming an attenuation α=0.27 km⁻¹ (following Horgan et al., 2011). This attenuation corresponds to a seismic quality factor (Q) of 30–300 for 10–100 Hz waves in a 3860 m s⁻¹ medium. H&A2009 noted that Eq. (4) is weakly dependent on uncertainties in α. Long-path multiples from shots in which the primary reflections were from the interface between ice and seismically thick water resulted in 60, 19, 9, and 24 estimates of A₀ for lines 1–4, respectively (left column Fig. 4, A_0MB columns Table 1).

2.4.2 Direct-path methods

Two methods were used to estimate source amplitude from the direct arrival amplitudes (B_i). Direct arrivals have successfully been used to determine source size (Muto et al., 2019) and to normalise shot records (Brisbourne et al., 2017). Following H&A2009,

$$\begin{array}{}\text{(5)}& B_i=A_{\mathrm{0}}{\mathit{\gamma }}_{\mathrm{d},i}{e}^{-\mathit{\alpha }{s}_{\mathrm{d},i}}\text{(Eq. 8, H\&A2009)},\end{array}$$

where B_i denotes direct arrival amplitude at receiver index i, and s_d,i and γ_d,i are the direct arrival path lengths and path amplitude factors. We first estimated A₀ using the direct-path pair method of H&A2009 (H&A2009, Eq. 9). This method uses receiver pairs where the ratio of path lengths $s_{\mathrm{2}}/s_{\mathrm{1}}=\mathrm{2}$ and where the offset is sufficient that depth-averaged attenuation can be assumed the same. This negates the need for an independent attenuation estimate. Our acquisition geometry did not result in pairs where $s_{\mathrm{2}}/s_{\mathrm{1}}=\mathrm{2}$ exactly, so an acceptance distance (x₁) was set such that pairs were used if $s_{\mathrm{2}}\ge \mathrm{2}{s}_{\mathrm{1}}-x_{\mathrm{1}}\wedge s_{\mathrm{2}}\le \mathrm{2}{s}_{\mathrm{1}}+x_{\mathrm{1}}$. We set x₁=14 m through trial and error, looking for the minimum x₁ that resulted in multiple estimates of A₀ for all shots. This resulted in a mode of eight pairs per shot (mean of 7.7, standard deviation of 3.7). A₀ direct pair estimates are shown in Fig. 4 (centre left column) and Table 1 (A_0DP columns).

We also investigated A₀ estimation using all direct arrival amplitudes by fitting the observed B_i values to Eq. (5) and minimising the misfit to determine optimal A₀ and α values. We refer to this method as the direct-path linear intercept method, because

$$\ln {\displaystyle \frac{B_i}{{\mathit{\gamma }}_{\mathrm{d},i}}}=-\mathit{\alpha }{s}_{\mathrm{d},i}+\ln A_{\mathrm{0}}$$

shows that in −s_i versus $\ln \frac{Bi}{{\mathit{\gamma }}_i}$ space every shot record should exhibit a common gradient (α) and independent y intercepts representing ln A₀. Despite this linear form, we solved for best-fitting parameters directly from Eq. (5) using non-linear regression. We restricted our direct arrival analysis to returns from offsets greater than 450 m, and testing up to an offset limit of >800 m did not result in significantly different A₀ and α estimates. For both direct-path methods, path effects (γ_d,i) were estimated both using Eq. (2) and estimating wavefront energy using ray theory (Sect. 6.2 of Shearer, 2009, modified to account for different outgoing and incoming angles). The wavefront energy approach did not result in better A₀ estimates, with a larger distribution and poorer correlation with the known reflector method. We therefore present results using Eq. (2), consistent with our other source size estimates. A₀ direct linear intercept estimates are shown in Fig. 4 (centre right column) and Table 1 (A_0LI columns).

2.4.3 Known reflector methods

Reflections from a known impedance contrast, in this case the floating ice shelf overlying the ocean cavity, allow another method of determining A₀. We estimated a best-fitting A₀ for each ice shelf shot by non-linear regression of Eq. (6) (Eq. 10, H&A2009).

$$\begin{array}{}\text{(6)}& R\left(\mathit{\theta }\right)={\displaystyle \frac{A_i}{A_{\mathrm{0}}}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\gamma }}_i}}{e}^{\mathit{\alpha }{s}_{\mathrm{d},i}}\text{(Eq. 10, H\&A2009)}\end{array}$$

We determined the optimal A₀ for each floating shot by minimising the root-mean-squared misfit between the reflection amplitudes resulting from the Zoeppritz equations for the seismic properties in Table 2 and the observed bed reflection amplitudes (A_i), Eq. (6). To account for the possibility that englacial debris may be present in the basal ice, we also optimised the seismic properties of the ice used in the Zoeppritz equations while keeping the underlying water properties constant. We allowed the basal ice to vary within a range encompassing debris contents of 0–20 % by volume. The range of seismic velocities for this basal ice was estimated using a Bruggeman mixing model following Röthlisberger (1972). We refer to this method as the known reflector method, and the resulting A₀ estimates are shown in Fig. 4 (right column) and Table 1 (A_0KR columns). The method resulted in the same number of A₀ estimates as the multiple-bounce method, and each line's average basal ice properties estimated during optimisation are shown in Table 3. The known reflector method requires an estimate of path effects but is insensitive to our assumption that α=0.27 km⁻¹ as the same α used to determine A₀ is later used in Eq. (6) to determine the basal reflection coefficient. The known reflector method has similarities to the technique used by Smith et al. (2018) in their study of the lithology beneath Subglacial Lake Ellsworth, although here we explicitly estimated A₀, allowed the basal ice properties to vary, and used amplitude versus offset techniques.

Table 2 Range of seismic properties assumed for the lower ice shelf. ν denotes Poisson's ratio.

Table 3 Seismic properties estimated in the lower ice shelf.

Figure 4 A₀ source size estimates for Whillans grounding zone lines 1–4 (rows) using four methods (columns). Left column: A₀ estimates using the primary–multiple amplitude ratio method. Centre left column: A₀ direct pair estimates. Centre right column: A₀ linear intercept estimates. Right column: A₀ estimates from known reflection coefficient method assuming ice overlying water. (See Fig. 1 for line locations.)

2.5 Choosing the best A₀

The known reflector method provided our best estimate of A₀ as judged by its potential to recover accurate estimates of basal reflectivity (e.g. ice–water reflection coefficient where the ice is known to be floating) and its narrow normal distribution (Fig. 4, Table 1). The narrow distribution indicates low source size variability, consistent with a uniform firn–ice profile, a consistent drilling depth and geophone placement, back filling all shots, and allowing at least 24 h before detonation.

Both our direct-path methods resulted in large standard deviations (Table 1) and correlate poorly with our known reflector estimates (r² (coefficient of determination) of 0.09 for the direct pair method and 0.04 for the linear intercept method, Fig. 5). The linear intercept method resulted in an average $\mathit{\alpha }=\mathrm{1.4}\pm \mathrm{0.5}$ km⁻¹ (mean and 1 standard deviation of the combined results for all four lines). Individual line average values range from 1.0–1.6 km⁻¹. These α estimates are an order of magnitude greater than commonly used published estimates and are not used in our analysis. The multiple-bounce method correlates well with the known reflector method (r²=0.46, Fig. 5). Linear regression of the known reflector estimates with the multiple-bounce estimates results in a best-fitting gradient of 2.2 with an intercept of 180. However, this relationship is dependent on our estimate of α and our γ estimates and will be discussed in Sect. 4.

Figure 5 A₀ estimate comparisons. (a) A₀ estimates from known reflector method against A₀ estimates from multiple-bounce method (coefficient of determination (r²) of linear regression=0.46). (b) A₀ estimates from known reflectivity method against A₀ estimates from the direct pair method (r²=0.09). (c): A₀ estimates from known reflector method against A₀ estimates from linear intercept method (r²=0.04).

2.6 Estimating subglacial properties

Using each line's A₀ values from the known reflector method (Table 1, Fig. 4 right column) we calculated the angle-dependent bed reflection coefficients for each shot gather (R(θ), Eq. 1). Our angle coverage typically extends up to 25^∘, with some shots extending to 30^∘. We present R(θ) as average values within 10^∘ of normal incidence (R_b) (Figs. 6a and 7a) to allow comparison with normal incidence methods reported elsewhere (e.g. Muto et al., 2019). We then calculated the optimal combination of subglacial seismic velocities (V_p,V_s) and density (ρ) (Figs. 6 and 7b–d) by fitting each shot's entire R(θ) to the Zoeppritz equations while imposing reasonable bounds for the subglacial material following Zechmann et al. (2018), expanded to allow for an ice–water interface (Table 4). During optimisation we imposed the additional constraint that the optimal V_p and V_s must result in a realistic Poisson's ratio (ν) of 0.25–0.5 (Hamilton, 1979). Optimisation minimised the root-mean-squared misfit between the observed amplitudes for each shot and those modelled by the Zoeppritz equations using the fmincon algorithm in MATLAB^®. This optimisation uses a trust region approach resulting in rapid convergence. We set the basal ice's seismic properties to those obtained for each line during our A₀ known reflector method in Sect. 2.4.3 (Table 3). We repeat our R_b estimates and the optimisation of V_p, V_s, and ρ values using R(θ) values estimated using all estimates of A₀ for each line, resulting in the same number of estimates of basal properties per shot as there are estimates of known reflector source size per line. In some cases our inversion repeatedly converged on the same solution, implying a misleadingly high precision. To account for this we also estimated our uncertainties by examining the retrieved basal properties from the floating portions of our survey. For all floating portions of the survey, misfit between the recovered properties and theoretical properties resulted in 1 standard deviation uncertainties for R_b of ±0.09, V_p of ±140 m s⁻¹, V_s of ±430 m s⁻¹, and ρ of ±30 kg m⁻³. Uncertainty estimates for each line are shown in Table 5.

Figure 6 Line 1 (a) seismic basal reflectivity at normal incidence estimated from the average value within 10^∘ (R_b). The red line shows radar basal reflectivity from Christianson et al. (2016). (b–d) Box plots of V_p, V_s, and ρ estimated using Zoeppritz fitting and all estimated source sizes. Blue boxes show the 25th and 75th percentiles, whiskers extend to cover data points, and outliers are plotted as black points. Solutions using the mean source size are overlain as black crosses. All estimates use source sizes obtained using the known reflector method. (e) Stacked active-source seismic reflection profile with ice flow from left (grounded ice stream) to right (floating ice shelf). Shot ghost denotes the short-path multiple generated by the ray path from the source up to the ice–air interface then down. For profile location see Fig. 1.

Figure 7 Lines 2 (left), 3 (middle), and 4 (right). (a) Seismic basal reflectivity at normal incidence estimated from the average value within 10^∘ (R_b). The red line shows radar basal reflectivity from Christianson et al. (2016). (b–d) Box plots of V_p, V_s, and ρ estimated using Zoeppritz fitting and all estimated source sizes. Blue boxes show the 25th and 75th percentiles, whiskers extend to cover data points, and outliers are plotted as black points. Solutions using the mean source size are overlain as black crosses. All estimates use source sizes obtained using the known reflector method. (e) Stacked active-source seismic reflection profile. Line 2 is plotted flowing from grounded (left) to floating (right). Lines 3 and 4 are plotted with flow into the page. Shot ghost denotes the short-path multiple generated by the ray path from the source to the ice–air interface and then down. O.c denotes the ocean cavity. For locations see Fig. 1.

Table 4 Seismic velocity (V_p, V_s), density (ρ), and Poisson's ratio (ν) bounds used for Zoeppritz fitting.

3 Results

3.1 Reflection coefficients and basal properties

Line 1 exhibits generally slowly varying R_b values upstream of the grounding zone, before an abrupt change at the grounding zone (Fig. 6). This change occurs over less than 500 m at approximately kilometre 9. V_p, V_s, and ρ values retrieved from Zoeppritz fitting exhibit a similarly abrupt change at the grounding zone. Upstream of the grounding zone binned-mode V_p and V_s values equal 2000 and 1100 m s⁻¹ respectively, and mode ρ values equal 1800 kg m⁻³. Kilometres 3–4 of line 1 exhibit retrieved V_s and ρ values similar to those expected for water, but R_b and V_p estimates suggest otherwise. In the floating portion of the profile most retrieved properties are equal to those expected for water (Table 5). Estimates of V_s are more spatially variable with larger distributions both upstream and downstream of the grounding zone. Line 2 exhibits similar patterns in R_b and retrieved seismic properties to line 1. An abrupt transition is observed at the grounding zone (kilometre 3.6), and the grounded and floating portions are dominated by distinct seismic properties (Fig. 7, left panel, Table 5). Upstream of the grounding zone two retrieved estimates exhibit properties similar to those of water (kilometre 0–0.5); however, neither are unambiguous. V_s estimates are again more variable than other parameters, with most floating shots exhibiting V_s values typical of water. Line 3 (Fig. 7, middle panel) shows both rapid and gradual changes in basal properties along the profile. Rapid changes are observed either side of kilometres 7–8 where a narrow bed feature exhibits V_p and ρ estimates typical of subglacial water. Kilometres 2–4 display a gradual change in R_b while the associated transition in V_p and ρ occurs abruptly over <500 m. V_s estimates are variable along the profile and exhibit scatter within regions thought to be both grounded (kilometres 0–3) and floating (kilometres 3.5–6). Line 4 (Fig. 7, right panel) is dominated by R_b, V_p, and ρ estimates typical of ice over water (kilometres 0–7). The transition from these values occurs over a distance of <1 km beginning at kilometre 7. As with the other profiles the estimates of V_s are variable, but most often the floating portion of the profile (kilometres 0–7) exhibits V_s estimates typical of water (Table 5).

Table 5 Binned-mode seismic properties estimated using normal incidence methods (R_b) and Zoeppritz fitting (V_p, V_s, and ρ) for the grounded and floating portions of each line. Bin sizes are shown in square brackets. The 1 standard deviation uncertainties were obtained from the misfit in the floating portion of each line.

3.2 Experiment timing and tidal elevation

Seismic shooting occurred at different stages of the tide, resulting in the potential for different tidal heights along profile. Shot and receiver elevations were not directly observed at the time of shooting, so instead we present tidal heights estimated at the floating end of the profile using Erofeeva et al. (2020) (Fig. 8). Figure 8a shows that kilometres 6–12.5 of line 1 were acquired on the falling tide when the tidal elevation varied from +0.1 to −0.6 m. The pronounced change in basal reflectivity that occurs at approximately kilometre 9 on line 1 (Fig. 6) does not coincide with a step in the tidal elevation. Other step changes in tidal elevation along line 1 also do not coincide with changes in basal reflectivity (e.g. kilometre 1, Fig. 6). Lines 2–4 all took less than a day to acquire and for the most part have no major step changes in tidal elevation along the profiles. An exception to this occurs on line 2 where the onset of high basal reflectivity (kilometre 3.6–4.1, Fig. 7, left panel) occurs in proximity to an offshore 0.3 m change in tidal elevation.

Figure 8 Shot timing and tidal elevation from Erofeeva et al. (2020). (a) Line 1. Top subplot shows the timing of shots (blue bars) overlain on the tidal elevation anomaly. Bottom subplot shows vertical tidal anomalies (Erofeeva et al., 2020) at the time of shooting as a function of distance along the profile. (b, c) Same as (a) but for lines 2, 3, and 4. Latitude (lat) and longitude (long) for each tide model time series are shown in each top subplot.

3.3 Repeat elevation profiles across the grounding zone

Repeat kinematic GNSS elevation profiles were acquired along lines 1 and 2 (Fig. 9) and have previously been used to validate the seismically imaged grounding line location (Horgan et al., 2013b). We locate the grounding zone using the standard deviation of elevation observations in 50 m spatial bins after the removal of a single best-fitting spline from each profile. Upstream of the grounding zone we expect this value to represent the method uncertainty, which comes from both the GNSS observations and our ability to repeat a track precisely, combined with a measure of the roughness of the surface. Downstream these combine with the displacement of the ice surface due to the tide. The grounding line is determined to be the point at which the standard deviation changes from values representative of grounded upstream values to those representative of floating values. The pick is subject to some interpretation as roughness and the ability to repeat a track can vary spatially and can correlate with surface slope (e.g. van der Veen et al., 2009).

Figure 9 Repeat kinematic profiling along lines 1 (a, b) and 2 (c, d). (a and c) The elevation (top), residual elevation after removal of a best-fitting spline (middle), and standard deviation of residual elevation in 50 m spatial bins (bottom). (b and d) The timing of the GNSS profile data collection (vertical overlain on the vertical elevation anomaly of Erofeeva et al. (2020).

Repeat elevation profiles for lines 1 and 2 were acquired on the rising tide. The tidal range for line 1 at the time we observed was approximately 1.5 m, while line 2 was observed during a range of approximately 0.35 m. Both profiles exhibit a region of relatively high surface slope that begins upstream of the onset of vertical tidal displacement. We pick the line 1 grounding line at kilometre 9.6 and at kilometre 3.6 for line 2. Well upstream of the grounding zone, our repeat tracks typically all fall within 0.1 m vertically of each other. At the resolution of our data we do not observe migration of the grounding line in the GNSS data.

4 Discussion

4.1 Subglacial properties beneath Whillans Ice Stream's grounding zone

Subglacial material beneath the grounded ice stream exhibits ρ and V_p values in the range of dilatant till but with most V_s values typical of those observed in dewatered tills (Figs. 6 and 7, Table 5) (Zechmann et al., 2018). Our estimates of V_p and ρ for all lines are close to those estimated by Luthra et al. (2016) in their active-source seismic study of a major sticky spot beneath Whillans Ice Plain. V_s estimates from the grounding zone are greater than those estimated by Luthra et al. (2016), although they overlap within uncertainties. When compared with estimates from upstream on Whillans, where Blankenship et al. (1986) measured V_s of 150±10 m s⁻¹, our results indicate significantly stiffer till beneath the grounding zone. Basal shear stress is already known to vary spatially beneath Whillans Ice Stream. Inversion of surface elevation, ice thickness, and remotely sensed velocity observations has resolved spatially variable basal shear stress (Joughin et al., 2004b), and spatially variable rates of change of basal shear stress during the ice stream's deceleration (Beem et al., 2014). Joughin et al. (2004a) estimated low basal shear stress near the grounding zone, similar to that observed elsewhere beneath the majority of the ice plain. Lipovsky and Dunham (2017) introduced spatially variable bed properties in their rate and state friction model to better reproduce the timing and distribution of stick–slip displacement on Whillans Ice Plain. Passive seismic and geodetic observations of Whillans Ice Stream's stick–slip motion have been used to locate asperities beneath the central portion of the ice stream (Walter et al., 2011) and at its grounding zone (Pratt et al., 2014).

The transition in basal properties at the grounding zone of Whillans Ice Stream is abrupt in both longitudinal lines (lines 1–2), occurring over distances of less than 500 m. This is less than the ice thickness of 730–790 m. The transverse lines (lines 3–4) exhibit less abrupt transitions but still show change over distances of less than 1 km. The rapid transition in basal properties suggests that even in the case of a fast-flowing ice stream with low basal shear stress such as Whillans, it is necessary to solve the full Stokes equations if the ice flow velocity field is to be accurately modelled across the grounding line (Pattyn et al., 2013). The radio echo sounding (RES) results of Christianson et al. (2016) provide additional insights (Figs. 6a and 7a). Lines 1, 3, and 4, which all sample the embayment in the grounding zone to the grid north (Fig. 1), all exhibit a drop in RES basal reflectivity of approximately 3–5 dB as the grounding zone is crossed from the grounded ice stream to the floating ice shelf. This change occurs over similar length scales to the seismically detected transition. In contrast, line 2, which crosses the peninsula to the grid south exhibits a gradual increase in RES basal reflectivity of approximately 10 dB after the ice goes afloat, over a distance of approximately 3 km. Christianson et al. (2016) attribute the differences in the RES-detected transitions to the presence of basal roughness (fluting, modelled with a 20 m wavelength and 4 m root-mean-squared heights) and entrained debris in the ice shelf in the embayment and a basal interface that is becoming smoother and losing the basal debris zone due to basal melt at the peninsula. The percentage of entrained debris we obtained during source size estimation is similar across all four lines (6 %–7 %), indicating differing debris content is unlikely to be the cause of the differences in RES basal reflectivity. MacGregor et al. (2011) reported low-frequency (2 MHz) RES bed reflectivity from elsewhere on Whillans Ice Stream and the adjacent Kamb Ice Stream and found negligible change in RES reflectivity when crossing the grounding zone. One possibility discussed by MacGregor et al. (2011) was the presence of brackish water upstream of the grounding line, smoothing the RES-imaged transition from grounded to floating ice.

4.2 Water upstream of the grounding zone

Upstream of the grounding zone several regions (e.g. line 1, kilometres 3–4; line 2, kilometre 0–0.5; line 3 kilometres 7–8) exhibit properties that indicate the presence of subglacial water, although not without ambiguity. This ambiguity likely results from water column thicknesses that are less than one-quarter the dominant seismic wavelength for our data ($\mathit{\lambda }/\mathrm{4}\approx \mathrm{3.6}$ m). Visual inspection of shot records shows that in these regions the thin-layer effects detailed by Booth et al. (2012) result in constructive and destructive interference of our basal wavelet, leading to best-fitting parameter combinations that are not representative of the contrast in properties. A similar phenomenon likely results in the anomalous estimated values at the grounding zone of line 1 (Fig. 6, kilometre 9) and for kilometre 7–7.5 of line 3 (Fig. 7). However, no similar attribution is possible for the V_s outliers in the floating portions of all lines, which instead appear to correspond to low signal-to-noise ratios apparent in visual inspection of the shot records.

Our seismic methods are insensitive to whether the subglacial water is sourced from beneath the ice stream or from the ocean cavity. The WISSARD field site was initially selected as it lay on the subglacial drainage path from Subglacial Lake Whillans to the ocean cavity (Fricker et al., 2010; Carter and Fricker, 2012). Oversnow geophysical surveying, including the data presented here and in Christianson et al. (2013), has shown the potential for estuarine flow across the grounding zone (Horgan et al., 2013a). Shot times, tidal stage, and bed reflectivity lack correlation between changes in tidal height and imaged bed properties. One exception to this occurs on line 2 where the change in bed properties at kilometre 3.6 (Fig. 6, left column) occurs in proximity to a 0.3 m change in tidal height at kilometres 3.8–4.2 (Fig. 8b). We consider this correlation coincidental as line 2's grounding line position appears pinned at kilometre 3.6 by an approximately 6 m change in bed elevation. Also, repeat GNSS profiling (Fig. 9c) indicates vertical change at line 2's grounding line is likely to be much less than that estimated offshore, and even a 0.3 m change in water column thickness would be insufficient to cause the pronounced change in reflectivity observed. Line 1's repeat GNSS profiling (Fig. 9a) locates the onset of vertical tidal deflection 0.6 km downstream of the seismically resolved change in subglacial properties. This indicates the presence of water upstream of the of the GNSS-picked grounding line, but the subjective nature of the GNSS method makes this conclusion tentative. Line 1's repeat GNSS profiling also suggests the region between kilometres 9.6–12 is a zone of ephemeral grounding, resulting in a smaller distribution of elevations over the observed portion of the tidal cycle (Fig. 9a, bottom subplot). Our experiment was not designed to study changing bed properties over a tidal cycle, which would be better examined using tilt metres or fixed GNSS stations and a fixed geophone deployment with a source repeating at the same location.

While our methods are not able to determine the process of stiffening at the grounding zone and ponding upstream, our observations are broadly consistent with the findings of several previous modelling studies. In the nomenclature of Sayag and Worster (2013), our study location appears to be a fixed-grounding-line, stiff-bedded system, although the zone of ephemeral grounding and the 0.6 km difference between our seismically determined grounding zone and that located by our repeat GNSS profiling shows some grounding line migration may be occurring on line 1. Our seismic properties indicate a stiff bed over thicknesses of at least approximately 5 m ($\mathit{\lambda }/\mathrm{4}=\mathrm{5}$ m for a 100 Hz wave in a 2000 m s⁻¹ medium). Estimated seismic velocities and densities imply Young's modulus (E) of 3.1–6.2 GPa in the subglacial material with lines 1, 2, and 4 all exhibiting E=5.2–6.2 GPa. Our observations at this location are not able to identify the asymmetric grounding line migration outlined by Tsai and Gudmundsson (2015). Local variations in bed and surface slope and ice thickness are likely to contribute to this; however the resolution of our GNSS method and our temporal sampling of basal properties also contribute to a lack of fidelity. Stiff till beneath the grounding zone and localised bodies of water upstream of the grounding zone are in keeping with the compression and dewatering of subglacial till due to ice flexure modelled by Walker et al. (2013). Stiffening of the till was also invoked by Christianson et al. (2013) as the cause of the enhanced internal deformation evident in radio echo sounding profiling across the grounding zone. The presence of isolated water bodies also aligns with the alternating pressure gradients causing barriers to water flow upstream of the grounding proposed by Sayag and Worster (2013) and the movement of water upstream of the grounding line modelled by Warburton et al. (2020). Warburton et al. (2020) show that low subglacial permeability should lead to filtering of the response of ice flow to tidal forcing. If this is true for Whillans Ice Stream, then the combination of the low till permeability suggested by our findings and the tidally modulated twice-daily stick–slip motion of the ice stream indicates its response to tides is not controlled by fluid connectivity through the grounding zone till.

4.3 Estimating source size (A₀)

Our preferred method of estimating source size is only possible when a portion of the survey area contains a known reflection interface. The interface need not be known exactly, as demonstrated by our retrieval of basal ice properties alongside estimating source size, provided the shape of the R(θ) response varies with changing properties along with the absolute level of reflectivity. Comparison with other methods used to estimate A₀ demonstrates the efficacy of the commonly employed multiple-bounce method (Fig. 5). A₀ estimated using the multiple-bounce method was, however, approximately twice that estimated using our known reflector method (Fig. 5). This difference can be reduced by a more thorough treatment of the path amplitude factor (γ_i). For instance, applying the geometric loss estimated by Margrave (2003) results in a best-fitting gradient of 1.6. The remaining difference can be accounted for by varying α in our known reflectivity A₀ calculation, with an α=6.0 km⁻¹ resulting in a 1:1 relationship between the multiple-bounce and known reflector methods, albeit with a linear intercept of approximately 100. Instead of using path amplitude factors from Margrave (2003) and adjusting our α estimate, we have chosen the inverse pathlength approach of Eq. (2) and a published α estimate for clarity and to better enable repeatability. The discrepancy between the methods indicates that attenuation (α) and path amplitude factors (γ) remain areas of uncertainty, overcome here by our use of a known reflector. In the absence of reliable A₀ estimates, other attributes of the amplitude reflection curve such as the angle of phase change (e.g. Anandakrishnan, 2003a) can be effective predictors of subglacial geology. Direct-path methods for A₀ estimation have been successfully employed elsewhere (Muto et al., 2019), and greatly simplify R(θ) recovery. Muto et al. (2019) presented data where the sources were buried at 40–50 m depth, compared to our 27 m, and their signal-to-noise ratios are high, as evident in their imaging of englacial seismic reflectivity. The poor correlation between our known-reflector and direct-path A₀ estimates (Fig. 4) shows that further investigation of direct-path methods is warranted. Both the direct-path methods we present would benefit from a greater offset distribution, and the direct pair method would benefit from a greater number of path combinations where $s_{\mathrm{2}}/s_{\mathrm{1}}=\mathrm{2}$ than was available to us. Trace interpolation could also be used here as the direct arrival energy is unlikely to change rapidly. Also, the path effects (γ_d,i) experienced by the direct ray are likely to be inadequately captured by our approach due to the possibility of unaccounted for energy loss and more complex travel paths than those predicted within the firn.

Our Zoeppritz fitting methodology is skilled at recovering both V_p and ρ as demonstrated in the floating portions of all lines where the recovered values are those expected for water (see Table 5 floating estimates). The methodology is less skilful at recovering V_s, likely due to the weaker dependence of the shape of the R(θ) curve on V_s for the angles we observe. Using average source sizes and the known reflector method, we recover the near-zero V_s typical of water for 73 of the 112 floating shots in our survey. Estimating V_s along with ρ allows the shear modulus to be estimated, which can be used to calculate the effective pressure in the till (Luthra et al., 2016). This provides a more direct link between seismic observations and till properties than is otherwise possible from estimates of normal incidence reflectivity (R_b) alone. An acquisition geometry that covered greater angles would improve our ability to estimate V_s; however, limitations due to interference from direct arrivals would still exist. These limitations could be overcome by observing much greater offsets, where direct arrivals no longer interfere with the bed return, or surveying in regions of greater ice thickness. Using multiple charge sizes and configurations also highlights the importance of source configuration. Line 3, which consisted of the largest charges by weight (0.85 kg), resulted in the lowest A₀ estimates calculated from both the known reflector method and the multiple-bounce method. The charges for line 3 were made up of a stack of a single 0.4 kg charge and three narrower 0.15 kg charges. These narrower charges were likely less well coupled with the shot hole wall, and the longer linear configuration resulted in a less effective source. A shorter interval between shot loading and detonation may have also been a factor here as line 3 was shot within 1–2 d of loading.

5 Conclusions

Subglacial material beneath Whillans Ice Stream's grounding zone is relatively stiff, exhibiting V_s≈1100 m s⁻¹ and Young's moduli of 3.2–6.2 GPa, making it more similar to a subglacial sticky spot than to deforming till. The transition from this stiff subglacial sediment to the ocean cavity is abrupt, occurring over distances of 500–1000 m. This seismically imaged transition differs from that imaged using RES, which detects both an abrupt transition and a gradual one at the embayment and promontory respectively (Christianson et al., 2016). Upstream of the grounding line we detect thin, apparently isolated, bodies of water. These findings are consistent with models that compact till within the grounding zone and those that isolate water upstream of the grounding zone, although we cannot detect whether the subglacial water is sourced from the ocean cavity or subglacially. Our comparison of methods used to determine source size (A₀) shows that the commonly employed multiple-bounce method correlates well with the known reflector method available to us. However, our comparison also highlights that path effects (γ_i) are incompletely modelled by the methods employed here and elsewhere. Our findings also reinforce the need for consistency in source placement, configuration, and time between burial and detonation. Overall our methods are skilled at retrieving basal properties at relatively high spatial resolution where the thickness of the subglacial material is sufficient to prevent thin film effects ($>\mathit{\lambda }/\mathrm{4}$). Both V_p and ρ are reliably retrieved, while V_s is recovered less consistently. While we are currently unable to accurately recover all seismic properties for what appear to be thin water layers, our methods do show promise here. These thin layers are pertinent for ice flow, and techniques such as full waveform inversion are likely to prove useful here. These methods, which invert not just for a single amplitude of the basal return but also for the full time series, have been successfully applied to other environments where thin layers with large contrasts in seismic properties have been investigated (e.g. Pecher et al., 1996).