Supplementary material for "Passive Seismology for Investigating the Ice-Bedrock Interface Zone of East Antarctica"

Supplementary Material for Chapter 5: The subglacial conditions of the Totten Glacier, East Antarctica, from horizontal-to-vertical spectral ratios (HVSRs) of seismic ambient noise This document contains the supplementary material for Chapter 5, including: •Texts S5.1 and S5.2; •Figures S5.1 to S5.20. Text S5.1: Workflow for geostatistical simulations of local bed topography To infer azimuthal variations of ice thickness (H(i)) within the substation HVSR footprint, we use the GStatSim Python package MacKie et al., 2023 to perform sequential Gaussian simulations of the local bed topography (Figs. S5.2–S5.11). Around each of the 10 stations analyzed in this study (Fig. 5.3), we define a 10 km x 10 km square grid, within which we select point measurements of bed topography from the collection of available RES data (Fig. 5.5a). For the ICECAP flight lines, bed elevation measurements are easily accessed from the Bedmap database (Frémand et al., 2023); for the INGV flight lines, only ice thickness measurements are provided, which we convert to bed elevation values with the use of the Reference Elevation Model of Antarctica (REMA; Howat et al., 2019). Following MacKie et al. (2023), the set of bed elevation measurements are then spatially interpolated across the grid with 50 m resolution. Significant non-stationarity is found in the gridded bed elevation data: we therefore detrend the gridded data using a radial basis function with a smoothing factor of 500 m, matching the coarse resolution of the Bedmap3 gridded product (Pritchard et al., 2025). The resulting grid of residuals are subsequently transformed into a Gaussian distribution, after which an isotropic variogram is calculated and modeled. The isotropic variograms estimated from the gridded residuals at most stations is well behaved, with the semi-variance smoothly increasing with lag towards a sill of ∼1 and clearly described by an exponential variogram model. For the 1 TI3A-C stations however, the semi-variance overshoots the sill and displays some cyclicity with lag, indicative of an underlying spatial trend (Gringarten and Deutsch, 2001). The presence of significant bed topography at the TI3A-C site is further highlighted by computing anisotropic variograms for the gridded residuals within the substation footprints (Fig. S5.13), mapping distinct variogram behavior with azimuth and suggesting a prevailing topographic trend along the 60°–150°axes. We therefore choose an exponential variogram model with anisotropy along the 150°azimuth for the geostatistical simulations at the TI3A-C stations, setting the major and minor ranges to 1000 m and 500 m respectively. For each station, a sequential Gaussian simulation is computed 50 times using the gridded bed elevation residuals and variogram model parameters, where the sampling space is defined by ordinary kriging with a 1000 m search radius and 100 conditioning points. The resulting ensemble comprises 50 realizations of local bed topography at 50 m resolution, which we then use (along with the station elevation from REMA) to estimate the variation in H(i)within the substation HVSR footprint (Figs. 5.5a, S5.12 and S5.20). Text S5.2: Workflow for HVSR inversions of subglacial structure To estimate the average thickness (H(s) S) and S-wave velocity (V(s) S) within the substation HVSR footprints of stations TI7A and TI8, we apply the HVSR inversion workflow developed by Kelly et al. (2025b) (contained in Chapter 4), which performs a misfit of HVSR peak and trough frequencies in a parameter search using the Neighborhood Algorithm (Sambridge, 1999). For each station, we identify HVSR peaks that are clear and consistent across the recording period (Fig. 5.7). The peak frequencies are then used to define frequency ranges for window rejection (Cox et al., 2020) applied to the single HVSR measurements (Fig. S5.16) computed across the available recording at each station (Fig. S5.1). The window rejection algorithm of Cox et al. (2020) iteratively removes windows from the HVSR based upon the lognormal statistics of peaks found within each frequency range, refining the sharpness of the peaks in the median HVSR and excluding windows where recorded seismic energies are not related to subsurface resonance. At TI8A for instance (Fig. S5.16b), the algorithm removes numerous windows with anomalously low and high HVSR amplitudes, likely attributed to local fieldwork noise and small data gaps. The inversion framework of Kelly et al. (2025b) involves evaluating the least-squares misfit (Ψ) between the peak and trough frequencies from the observed HVSR (γ(o) k) against those from the modeled HVSR (γ(m) k), i.e., Ψ = N+M X k=0 hγ(o) k−γ(m) ki2,with N=N′(p′≥P), M =M′(t′≥T),(S5.1) where krepresents the relative order in (ascending) frequency of each peak/trough feature from 1 to N+M, with Nand Mrespectively denoting the number of peaks and troughs selected in the observed HVSR (Fig. S5.16). N′and M′respectively define the number of peaks and troughs in the modeled HVSR, which are filtered according to their prominence p′and t′against set prominence values P= 0.4and T= 0.2to ensure that sufficiently strong features are selected from the modeled HVSR (Kelly et al., 2025b). We generate modeled HVSRs using the body2 wave resonance (BWR) model of Herak (2008), which was determined by Kelly et al. (2025a) (contained in Chapter 3) and Kelly et al. (2025b) to be well suited for representing Antarctic HVSRs. We also set the condition γ(p) k=(γ(p) k′if N′(p′≥P) = Nand M′(t′≥T) = M, ⊥else,(S5.2) where k′represents the relative frequency order of peak and trough features in the modeled HVSR, similar to k. Eq. S5.2 ensures that the observed multi-modal resonance pattern from the crystalline basement is replicated in the modeled HVSR within a specified frequency range (Fig. S5.16), recognizing that prominent resonances from other subsurface discontinuities in the Antarctic IBIZ are not expected (Kelly et al., 2025a). The initial sampling for the Neighborhood Algorithm search (Section 5.3.2) involves evaluating Eqs. S5.1 and S5.2 across the parameter space for the subglacial layer, gridded into 100 m (for H(s) S) and 50 m/s (V(s) S) intervals. The parameter space is consequently reduced to small regions (Fig. S5.17) where the condition given in Eq. S5.2 is satisfied. The corresponding set of samples are then used to initialize the Neighborhood Algorithm search (Sambridge, 1999), instead of the standard approach of randomly sampling the full parameter space. The Neighborhood Algorithm employs Voronoi cells to approximate a misfit surface over the parameter space, which is used to preferentially guide sampling within lower-misfit regions and refine the Voronoi tesselation in the process. In our case, we set arbitrarily high misfit values to the excluded regions of the parameter space to force the search towards viable solutions of Eq. S5.2 (Fig. S5.17). For the Neighborhood Algorithm, we set the search parameters ni= 20 (i.e., the number of iterations), nr= 25 (i.e., number of lowest-misfit Voronoi cells selected for resampling), and ns= 50 (i.e., number of new samples generated), as per Kelly et al. (2025b). The behavior of the parameter search is illustrated in (Fig. S5.18). To incorporate uncertainties in H(i), the initial sampling and parameter search is performed for each realization of simulated bed topography (Text S5.1), where H(i)is estimated for each realization from the difference between the REMA station elevation and the average of simulated bed elevation values within the substation HVSR footprint and allowed to vary according to the associated standard deviation of bed elevation values. Uncertainties in V(i) Sare additionally included by allowing V(i) Sin the parameter search to vary within 1800–2000 m/s (Section 5.2.1). The set of solutions determined from the Neighborhood Algorithm (i.e., the red lines in Fig. 5.9) are used to estimate the mean and standard deviation in H(i),V(i) S,H(s) S, and V(s) S, given in Table 5.1. 3 Table S5.1: Summary of the subsurface parameters used in our HVSR inversion to estimate the subglacial low-velocity zones at stations TI7A and TI8A-C. H= layer thickness, VS= S-wave velocity, VP= P-wave velocity, ρ= density, and QSand QP= Sand P-wave Q-factors. Refer to Kelly et al. (2025b) for an explanation of these parameter choices. Layer H[m] VS[m/s] VP[m/s] ρ[kg/m3]QS, QP Ice sheet Taken from geostatistical simulations 1800–2000 3870 920 100, 100 Subglacial layer 0–5000 200–2800 1330–4720∗1520–2490∗100, 100 Basement n/a 3000 5050∗2540∗100, 100 ∗Values approximated from corresponding VSusing Brocher (2005). 4 TI1A TI3A TI3B TI3C TI4A TI6A TI6B TI6C TI7A TI8A TI8B 24/12 31/12 07/01 14/01 21/01 28/01 Date (day/month) TI8C 2018 2019 17889 17896 17903 17910 17917 17924 Component: 1 Component: 2 Component: Z 5 Figure S5.1: (Previous page.) Overview of the seismic data collected at Totten Glacier during the 2018/2019 austral summer (Fig. 5.3), trimmed to exclude periods during instrument deployment and collection. Component waveforms are differentiated by color, whilst the black arrows above indicate prominent teleseismic earthquakes (i.e., magnitudes >6) recorded by the stations. 6 7 Figure S5.2: (Previous page.) A summary of the geostatistical simulations of local bed topography performed at station TI3A (see Text S5.1 for more details). (a) and (b) display the normal score transformation of the gridded residuals in bed elevation, detrended using a radial basis function. The resulting isotropic variogram and fit to different models is shown in (c) where, for the TI3A-C stations, anisotropy is evident and further investigated (Fig. S5.13). (d) presents one of the 50 realizations of simulated bed topography, with the ensemble mean and standard deviation shown alongside in (e) and (f), respectively. Station TI3A is located by the red triangle, with other nearby stations shown as gray triangles. (g) provides a closer look into (d) within the substation HVSR footprint, which is denoted by the dashed red circle. (h) and (j) illustrate the corresponding bed elevation and uncertainty in Bedmap3. In (d) and (g), the flight lines are colored according to Fig. 5.5a. 8 Figure S5.3: A summary of the geostatistical simulations of local bed topography performed at station TI3B (red triangle). Refer to Fig. S5.2 for plot details. 9 Figure S5.10: A summary of the geostatistical simulations in local bed topography performed at station TI8B (red triangle). Refer to Fig. S5.2 for plot details. 16 Figure S5.11: A summary of the geostatistical simulations in local bed topography performed at station TI8C (red triangle). Refer to Fig. S5.2 for plot details. 17 14501500 Ice thickness H ( i ) [m] 0 30 60 90 120 150 180 Azimuth [°] TI8A (a) 14501500 Ice thickness H ( i ) [m] TI8B 14501500 Ice thickness H ( i ) [m] TI8C 16501700175018001850 Ice thickness H ( i ) [m] 0 30 60 90 120 150 180 Azimuth [°] TI6A (b) 16501700175018001850 Ice thickness H ( i ) [m] TI6B 16501700175018001850 Ice thickness H ( i ) [m] TI6C 21502200225023002350 Ice thickness H ( i ) [m] 0 30 60 90 120 150 180 Azimuth [°] TI3A (c) 21502200225023002350 Ice thickness H ( i ) [m] TI3B 21502200225023002350 Ice thickness H ( i ) [m] TI3C Figure S5.12: Azimuthal variations in substation ice thickness (H(i)) for the (a) TI8A-C, (b) TI6A-C, and (c) TI3A-C stations. Refer to Fig. 5.5 for plot details. 18 0 250 500 750 1000 1250 1500 1750 2000 Lag [m] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Semi-variance TI3A 0 250 500 750 1000 1250 1500 1750 2000 Lag [m] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Semi-variance TI3B 0 250 500 750 1000 1250 1500 1750 2000 Lag [m] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Semi-variance TI3C 0° (Geo. N) 30° 60° 90° 120° 150° Figure S5.13: Anisotropic variograms for the TI3A-C stations, with the azimuth taken at 30° intervals clockwise relative to geographic north. A clear distinction in the range and sill is observed for 60°and 150°compared to other azimuths. Refer to Text S5.1 for details. 19 0.1 1 2 3 Frequency [Hz] 1 2 HVSR TI8A (a) 0.1 1 2 3 Frequency [Hz] TI8B 0.1 1 2 3 Frequency [Hz] TI8C 0.1 1 2 3 Frequency [Hz] 1 2 HVSR TI6A (b) 0.1 1 2 3 Frequency [Hz] TI6B 0.1 1 2 3 Frequency [Hz] TI6C 0.1 1 2 3 Frequency [Hz] 1 2 HVSR TI3A (c) 0.1 1 2 3 Frequency [Hz] TI3B 0.1 1 2 3 Frequency [Hz] TI3C 24/12 31/12 07/01 14/01 21/01 28/01 2018 2019 Figure S5.14: Median HVSR curves from all daily timeseries analyzed at the (a) TI8A-C, (b) TI6A-C, and (c) TI3A-C stations. Refer to Fig. 5.7 for plot details. 20 Figure S5.15: The azimuthal characteristics of the median HVSRs for the (a) TI8A-C, (b) TI6AC, and (c) TI3A-C stations. Refer to Fig. 5.8 for plot details. 21 Figure S5.16: The single HVSR measurements used for the Neighborhood Algorithm search of subglacial low-velocity zones beneath the TI7A and TI8A-C stations, as outlined in Section 5.2.3 and Text S5.2. The red and gray HVSR curves respectively illustrate the rejected and accepted windows from the window rejection algorithm of Cox et al. (2020), based upon the lognormal statistics of HVSR peaks within pre-defined frequency ranges (horizontal orange bars). The orange diamonds and crosses identify the peaks and troughs in the median HVSR evaluated in Eqs. S5.1 and S5.2 within the frequency range delineated by the gray shaded regions. Refer to Fig. 5.6 for other plot details. 22 Figure S5.17: The sampling of the parameter space used to initialize the Neighborhood Algorithm search of subglacial low-velocity zones beneath the TI7A and TI8A-C stations, exemplified for one of the 50 realizations of simulated bed topography to define H(i)(Text S5.2). The gray and colored cells respectively denote regions of the parameter space that are excluded and validated by Eq. S5.2, with the color-mapping illustrating the respective misfit values from Eq. S5.1. The red star locates the solution determined from the parameter search, as in Fig. S5.18. 23 37 250 500 1000 Number of samples 0.1 0.2 0.3 Misfit Ψ TI7A (a) Search 123 VS [km/s] 1.5 2.0 2.5 3.0 Depth [km] (b) Profiles 0.1 1 2 3 Frequency [Hz] 0 1 2 3 4 HVSR (c) HVSRs 0 2 4 H ( s ) [km] (d) Parameter spread 0.15 0.30 Misfit Ψ 1 2 V ( s ) S [km/s] 46 250 500 1000 Number of samples 0.04 0.06 0.08 0.10 0.12 Misfit Ψ TI8A 123 VS [km/s] 1.4 1.5 1.6 1.7 1.8 Depth [km] 0.1 1 2 3 Frequency [Hz] 0 1 2 3 4 HVSR 0 2 4 H ( s ) [km] 0.06 0.12 Misfit Ψ 1 2 V ( s ) S [km/s] 47 250 500 1000 Number of samples 0.05 0.10 0.15 0.20 Misfit Ψ TI8B 123 VS [km/s] 1.4 1.5 1.6 1.7 1.8 Depth [km] 0.1 1 2 3 Frequency [Hz] 0 1 2 3 4 HVSR 0 2 4 H ( s ) [km] 0.1 0.2 Misfit Ψ 1 2 V ( s ) S [km/s] 46 250 500 1000 Number of samples 0.05 0.10 0.15 0.20 Misfit Ψ TI8C 123 VS [km/s] 1.4 1.5 1.6 1.7 1.8 Depth [km] 0.1 1 2 3 Frequency [Hz] 0 1 2 3 4 HVSR 0 2 4 H ( s ) [km] 0.1 0.2 Misfit Ψ 1 2 V ( s ) S [km/s] 24 Figure S5.18: (Previous page.) An overview of the Neighborhood Algorithm search, exemplified with one of the 50 realizations of simulated bed topography (Text S5.2). For each station, four sub-plots are shown: (a) the misfit behavior of Eq. S5.1 for the entire population of samples, including those from the initial sampling (gray squares; Fig. S5.17) and those generated by the parameter search (pink circles), which converge towards an optimal solution (red star); (b) the corresponding subsurface profiles (gray and pink lines), compared against an ice-basement configuration as described in Fig. 5.9; (c) the corresponding forward-modeled HVSRs (gray and pink curves), with other HVSR features as described in Fig. S5.16; and (d) the corresponding parameter spread of H(s)and V(s) Swith the misfit value, with markers as described in (a). 25