4. Discussion

The stacked theoretical autocorrelograms together with the Bedmap2 ice thickness variations across three profiles are shown in Figs 5 and S5. There is good agreement between the observed and the theoretical arrivals of the reflection responses for most stations. The relative errors of the observed first P wave reflections to the theoretical ones are within 5, 10 and 15% at a total of 19, 27 and 32 stations, respectively. Similarly, the relative errors for S wave reflections for 17, 27 and 29 stations are within 5, 10 and 15% threshold, respectively. In addition, the variations between the Bedmap2 ice thicknesses and the reflection arrivals for stations along the three profiles are closely correlated. Given that the Bedmap2 ice thicknesses adopted to generate theoretical teleseismic waveforms are associated with certain uncertainties at these stations (Fretwell and others, Reference Fretwell2013; Yan and others, Reference Yan2018), we can conclude that the autocorrelation method can effectively image both the P wave and the S wave reflection responses of the ice-bedrock interface.

Fig. 5. Synthetic stacked radial and vertical autocorrelograms for each station. Panel (a) shows Bedmap2 ice thickness variations along the three profiles and the reference thicknesses (red dots) used to build models. Panel (b) displays synthetic vertical autocorrelograms. Yellow and blue crosses denote the arrival times of the observed first and second P wave reflection responses (corresponding to the yellow and blue circles shown in Fig. 4a). Panel (c) is similar to panel (b), but for radial autocorrelograms. Model parameters used to generate theoretical teleseismic waveforms are comprised of Bedmap2 ice thickness, v _p and v _s values. v _p is set to 3800 m s^–1 referring to previous studies and v _s is deduced from the relationship between v _p and v _p/v _s ratio.

Given the lack of ‘ground truth’ ice thickness from direct measurement methods such as drillings, calculated ice thickness is validated by comparing the ZAC estimate with the corresponding Bedmap2 ice thickness at each site. Differences between the ZAC estimates and the Bedmap2 ice thickness are within 100, 200, 300 and 400 m for 15, 24, 29 and 33 stations, respectively. The relative errors of the ZAC estimates to the Bedmap2 ice thickness are the same as those of the observed arrival times to the theoretical arrival times because the same v _p value is used in calculations of ice thicknesses and models to generate theoretical autocorrelograms. Our estimates show a good agreement with results obtained from Pham and Tkalčić (Reference Phạm and Tkalčić2018), which investigated Antarctic ice-sheet properties using an autocorrelation method with an improved scheme of selecting the spectral whiten width. Comparison shows, for the 19, 27 and 33 out of 38 stations that exhibit clear P wave reflections in both studies, the relative errors are within 5, 10 and 15%, respectively (Fig. 6). Note that instead of using ice thickness results estimated by Pham and Tkalčić (Reference Phạm and Tkalčić2018) directly, we adopt their P wave arrival times and use the same P wave velocity value as in our study to calculate ice thicknesses since different v _p values were used in two studies. Though the relative errors of ice thickness at stations P061 and N190 appear over 50% (Fig. 6), we find that the arrival times are shown in autocorrelograms in Pham and Tkalčić (Reference Phạm and Tkalčić2018) and our studies are in fact very close. The differences are caused by their recording mistakes as the arrival time values recorded in the table are not consistent with the ones shown in the autocorrelograms (Pham and Tkalčić, Reference Phạm and Tkalčić2018).

Fig. 6. Comparison of the ZAC ice thickness estimates with results obtained using the HVSR and PRF methods in our previous studies (Yan and others, Reference Yan2017, Reference Yan2018). The grey horizontal line in the plot indicates the average ice thickness of the HVSR, PRF and ZAC (this study) estimates for each station. The red circle and its bar represent the Bedmap2 ice thickness and its associated uncertainty for each station (Fretwell and others, Reference Fretwell2013). The ZAC estimates measured by Pham and Tkalčić (ZAC) are also displayed here (yellow stars) (note that the thickness values are a product of P wave arrival times taken from Pham and Tkalčić (Reference Phạm and Tkalčić2018) and v _p adopted in this study; we don't use their ice thickness values as two slightly different v _p values are used in different studies).

As no comparison among ice thickness estimates obtained from the PRF, HVSR and autocorrelation methods has been conducted before, we compared the ZAC estimates with results derived from the PRF and HVSR methods in our previous studies (Yan and others, Reference Yan2017, Reference Yan2018). We have calculated the average ice thickness at the 35 stations that successfully acquire ice thickness estimates using the PRF, HVSR and autocorrelation methods simultaneously. It shows that the ice thickness estimates obtained from the three methods are consistent as the differences of ice thickness values estimated by each method from the average of their ice thickness values are within 200 m at ~30 stations. Table S2 displays the relative errors of the ice thickness estimates with respect to the Bedmap2 ice thickness. We find that the level of relative error in the estimates is also consistent – the relative errors are within 15% for 25 stations, including GM01, P061, BYRD and ST01. This consistency can also be found for the seven stations with relatively large relative errors to the Bedmap2 ice thicknesses (over 15%), such as the N206, GM03 and ST09 stations. Similar disagreements with Bedmap2 ice thicknesses at these stations have also been found by Pham and Tkalčić (Reference Phạm and Tkalčić2018).

It is possible that the relatively large deviations from Bedmap2 ice thicknesses at the seven stations could be attributed to uncertainties in Bedmap2 ice thicknesses. More accurate ice thicknesses could be required to validate these specific results. The deviations could also be caused by the complex subglacial settings beneath these stations. As pointed out by Chaput and others (Reference Chaput2014), it is difficult to model the converted phases from the base of the ice sheet in their PRF study for a few stations including ST04, ST06 and ST09, which could possibly be influenced by anisotropy or basal dip structures. Therefore, we calculated autocorrelograms using synthetic teleseismic waveforms derived from models assuming an ice-sheet layer atop an inclined bedrock layer with dip angles varying from 5 to 30°, to show the effect of nonplanar basal topography on autocorrelograms. Figure S6 shows the theoretical teleseismic waveforms with a range of back azimuths between 0 and 360° generated using the raysum codes (Frederiksen and Bostock, 2000). Figure S7 displays the resulting autocorrelograms obtained using the theoretical teleseismic waveforms with the processing strategies mentioned in section 2. We can find that the autocorrelation reflections from different back azimuths show clear variations in values of amplitudes and arrival times. Such features become more remarkable as the dip angle increases. For example, the arrival time of the vertical autocorrelation reflections shown in the left panel of Fig. S7 decreases from 1.55 to 1.40 s when the dip angle increases from 5 to 15°. In addition to decreasing arrival times, the autocorrelations become more complex when the dip angle exceeds 20°. The same characteristics can also be found in radial autocorrelograms. With the PWS technique applied, we found that the azimuthal effect on the autocorrelograms resulting from the inclined basal structures is suppressed as a consistent autocorrelation reflection from different back azimuths is obtained (Fig. S8). Arrival times identified from the consistent autocorrelograms only reflect an average ice thickness beneath a station. Therefore, basal dips could certainly affect the arrival times of P and S wave reflection responses and consequently lead to errors of ice thickness estimates (and v _p/v _s ratios).

The PRF and the autocorrelation methods share certain similarities in data selection criteria and in the underlying equations in the frequency domain

(5)$$ZAC\lpar w \rpar = P^{^\ast}\lpar \omega \rpar P\lpar \omega \rpar $$

(6)$$PRF\lpar \omega \rpar = \displaystyle{{P^{^\ast}\lpar \omega \rpar S\lpar \omega \rpar } \over {\max \lcub {P^{^\ast}\lpar \omega \rpar P\lpar \omega \rpar ,\; \,c\cdot P_{\max}^{^\ast} P_{\max}} \rcub }}$$

where P(ω) and S(ω) are the vertical and radial component spectra from teleseisms at angular frequency ω, and P*(ω) is the complex conjugation of P(ω) (Sun and Kennett, Reference Sun and Kennett2016). The water level parameter c in Eqn (6) is used to avoid division by very small numbers by reducing the magnitude of the trough in the vertical component spectrum (Clayton and Wiggins, Reference Clayton and Wiggins1976). For small c, Eqn (6) approaches deconvolution, while it approaches a scaled cross-correlation as c becomes large. Considering the similarities of the PRF and autocorrelation methods, we compare their results for ice thickness and v _p/v _s ratios. It turns out that ice thickness estimates obtained from the PRF and autocorrelation methods are closely matched (for the total of 38 stations that are common to both methods, the relative error of the ZAC estimates to the PRF estimates for 22 and 34 station are less than 5 and 10%, respectively). Taking the uncertainties in the estimated v _p/v _s ratios into account, the values of v _p/v _s ratios determined using PRF and autocorrelation methods are generally consistent and ~2.00 for most stations (see Fig 4c). The differences are less than 0.1 for 17 stations, while the differences for stations such as BENN, ST02, ST07 and ST08 exceed 0.35. Further analysis indicates that the quality of PRF waveforms used to conduct H-Kappa stacking (Zhu and Kanamori, Reference Zhu and Kanamori2000) for the four stations is poorer than that of the autocorrelograms (Fig. S4). As is commonly acknowledged, the inherent trade-off between the thickness and shear-wave velocity for the PRF method also contributes to uncertainties in the v _p/v _s estimates.

To clearly show the capacity of the two methods for detecting ice properties, we have built a set of models (Fig. S10) to compare results obtained from both PRF and autocorrelation methods (Anandakrishnan and Winberry, Reference Anandakrishnan and Winberry2004; Chaput and others, Reference Chaput2014). Conclusions can be drawn by comparing results using different pairs of models (Figs 7, 8 and S11–S13). Table S1 displays the theoretical arrival times of P and S wave reflection responses calculated from model parameters (Eqns (7–9)), as well as the observed arrival times derived from the synthetic autocorrelograms using a bootstrapping technique (Fig. S14).

(7)$$t_{{\rm p(s)}}^{{\rm ice}} = \displaystyle{{2H_{{\rm ice}}} \over {v_{{\rm p}({\rm s})}}}\; $$

(8)$$t_{{\rm p(s)}}^{{\rm ice + sed}} = \displaystyle{{2\lpar {H_{{\rm ice}} + H_{{\rm sed}}} \rpar } \over {{\bar{v}}_{{\rm p}({\rm s})}}}$$

(9)$$\bar{v}_{{\rm p}({\rm s})} = \displaystyle{{(H_{{\rm ice}} + H_{{\rm sed}})v_{{\rm p}({\rm s})}^{{\rm ice}} v_{{\rm p(s)}}^{{\rm sed}}} \over {H_{{\rm ice}}v_{{\rm p(s)}}^{{\rm sed}} + H_{{\rm sed}}v_{{\rm p(s)}}^{{\rm ice}}}} $$

where H _ice and H _sed are ice and sediment thicknesses, and $t_{{\rm p(s)}}^{{\rm ice}}$, $t_{{\rm p(s)}}^{{\rm ice + sed}}$ are the theoretical arrival times of P (S) wave reflection responses in ice, and in ice and sediment layers. v _p(s) is P (S) wave velocities in ice, and $\bar{v}_{{\rm p(s)}}$ is the weighted average P (S) wave velocities in ice and sediment layers. Figure 7 shows the results obtained from model 1 and model 2 assuming two ice-sheet layers with different values of v _p/v _s ratios. It is clear that the ice thickness and v _p/v _s ratio estimates obtained from the PRF and H-Kappa methods are closer to the theoretical values than those obtained from the autocorrelation method (Fig. 7). The deviations of the theoretical v _p/v _s ratios from the estimated ones are attributed to the differences between the theoretical and the observed P and S wave arrival times. For example, the theoretical P wave and S wave reflection arrival times of model 1 are 1.55 and 2.94 s, while the observed values are 1.4 and 2.83 s (Table S1). Figures 8 and S11 show that both the autocorrelograms and the PRF waveforms become complicated when sedimentary layers are involved. The H-Kappa technique failed to acquire optimum estimates of ice thickness and v _p/v _s ratios using PRF waveforms as the phases from the sediment-bedrock interface greatly obscure phases from the ice-sediment interface (Figs S11–S13). A comparison of the theoretical and the observed arrival times of P or S wave reflection responses indicates that the autocorrelation method can only identify interfaces associated with strong impedance contrasts (Table S1). For example, Figure S10 shows that the P wave impedance contrasts (z = ρ · v _p, where z is the impedance contrast, and ρ and v _p are density and P wave velocity of the medium) set at the sediment-bedrock interfaces for models 4 and 6 are stronger than those set at the ice-sediment interfaces for models 3 and 5. Consequently, the observed P wave arrival times identified from models 3 (2.32 s) and 5 (2.35 s) are (nearly) equal to the theoretical times that were calculated with both the ice sheet and the sediment layers included (2.35 s for models 3 and 5). By contrast, the observed P wave arrival times of models 4 (1.46 s) and 6 (1.50 s) correspond to the theoretical times calculated with the single ice sheet layer (1.55 s). Similar results can be found for radial components (Table S1). Comparing results for models 3 and 4, we can find that the observed S wave arrival time of model 3 (4.73 s) is close to the theoretical arrival time (4.46 s) calculated including both the ice sheet and the sediment layers, while the observed arrival time of model 4 (2.77 s) approaches the theoretical arrival time (2.94 s) with the ice-sheet layer alone. We can also find that the observed S wave arrival time of model 6 (3.38 s) corresponds to the theoretical arrival time with the single ice-sheet layer (3.45 s) as the v _s value in ice decreases from 2.04 km s^–1 in model 4 to 1.74 km s^–1 in model 6, thus increasing the impedance contrast at the ice-sediment interface. Pham and Tkalčić (Reference Phạm and Tkalčić2018) pointed out that the low velocity of sediment results in low impedance contrast, explaining the relative transparency of sediment below an ice-sheet layer. In this sense, the autocorrelation method should be applied with caution to study low velocity sediment layers in glaciological settings. Alternatively, a nonlinear inversion method that is similar to the PRF waveform inversion approach could be effective to obtain the optimum model that fits with the observations.

Fig. 7. Estimates of v _p/v _s ratios obtained from the PRF and autocorrelation methods using theoretical teleseimic waveforms based on models 1 and 2 (panels (a) and (e)). Panels (b), (c) and (d) show the stacked autocorrelograms, PRF waveforms and H-Kappa results for model 1, and the right panels display results for model 2. The arrival times used to calculate v _p/v _s ratios in panels (b) and (d) are automatically picked using the bootstrapping technique (Koch, Reference Koch1992) shown in Fig. S14. The v _p/v _s ratio estimate obtained from PRF and H-Kappa methods is closer to the theoretical value than that measured using the autocorrelation method.

Fig. 8. Effect of a sediment layer (model 3) on v _p/v _s estimates obtained from the PRF and autocorrelation methods. Panels (b), (c) and (d) present the stacked autocorrelograms, PRF waveforms and the H-Kappa results for model 1, and the right panels show results for model 3. This comparison illustrates that neither the PRF method nor the autocorrelation method can clearly identify the ice-sediment interface (Table S1). The estimated v _p/v _s ratios obtained from the PRF and autocorrelation methods also deviate from the theoretical values.

The modeling results do not suggest that the autocorrelation method is superior to the PRF method, but it is still an effective alternative to the PRF method as it only requires a single component, which expands its application to datasets that were collected using single-component sensor decades ago. Moreover, the autocorrelation method can constrain both P wave and S wave velocity structures independently if both vertical and radial component records are available, while the PRF method can only reveal S wave velocity profiles. It is worth mentioning that there is a promising way to better constrain the thickness, v _p and v _s of an ice-sheet layer as a recent study has proposed a new H-Kappa-v _p stacking algorithm, which exploits independent constraints from the P wave autocorrelation and receiver functions (Delph and others, Reference Delph, Levander and Niu2019).

Unlike the close relationship between the PRF and autocorrelation methods, there are more differences than similarities between the HVSR method and the autocorrelation method. First of all, seismic ambient noise data are usually adopted for HVSR processing, though earthquake records can also be used. Secondly, the spectrum peak observed by the HVSR method represents S wave resonance frequency and only reflects S wave energy trapped in a medium, while the autocorrelation method can reveal P and S wave energy separately. Both the autocorrelation and PRF methods require a certain number of teleseismic earthquakes, and such data usually take a relatively long time to collect in Antarctica. By contrast, the HVSR method applied to ambient noise records only requires a recording time of several hours, making it a valuable complement to the autocorrelation and PRF methods (Yan and others, Reference Yan2018). Furthermore, the HVSR method combined with an inversion algorithm can be used to reveal inner stratification structures of ice sheets (Yan and others, Reference Yan2018), which the autocorrelation method has failed to detect (Pham and Tkalčić, Reference Phạm and Tkalčić2017).

As the v _p/v _s ratio is closely related to Poisson's ratio, we also calculated Poisson's ratio at each station using the relation shown in Eqn (10) (Christensen, Reference Christensen1996). The computed Poisson's ratios range from 0.32 to 0.38 (Table 1) and the average value (0.35) is slightly higher than the fixed value 0.33 used in previous studies (Gudmundsson, Reference Gudmundsson2011; Rosier and others, Reference Rosier, Gudmundsson and Green2015; Ramirez and others, Reference Ramirez2016).

(10)$$\upsilon = \displaystyle{{{(v_{\rm p}/v_{\rm s})}^2-2} \over {2{(v_{\rm p}/v_{\rm s})}^2-2}}$$

Variations in the value of Poisson's ratio at different sites could be attributed to temperature differences in different regions. Thus, rather than using a constant value for Poisson's ratio, it is more reasonable to consider a variable Poisson's ratio for large-scale ice-sheet modeling.