TY - JOUR
AU - Nyblade,, A.
AB - Abstract Reverberations of teleseismic compressional (P-) waves within a glacier or ice sheet may mask signals associated with crustal structure beneath the ice. We remove the signal associated with the ice from teleseismic P-waves using a wavefield downward continuation and decomposition technique that depends on known ice layer properties such as ice thickness, velocity, and attenuation. We test the method using data from nine stations in Antarctica and one station in Greenland. We deconvolve the downward-continued seismic wave vectors to create P-wave receiver functions that minimize the ice-layer reverberations in order to better measure signals from deeper structures. The subsurface P-wave receiver functions have similar sensitivities to crustal structure as those calculated from stations installed on bedrock. Synthetic experiments indicate subsurface P-wave receiver functions can constrain crustal structure more tightly than surface P-wave receiver functions when ice layer properties are known. We model the subsurface P-wave receiver functions using a Markov chain Monte Carlo inversion and constrain the product of crustal thickness and the column-average crustal-slowness beneath the stations. Our subglacial shear speed and thickness estimates are consistent with previous investigations at most stations. At station SUMG in south-central Greenland, our results suggest a thicker crust than from previous estimates. Antarctica, Arctic region, Waveform inversion, Crustal imaging, Wave propagation 1 INTRODUCTION Polar glacier acceleration and melt have contributed to recent sea level rise and are projected to dominate melt processes in the coming decades to centuries (Huybrechts & de Wolde 1999; Thomas et al. 2004; Gregory & Huybrechts 2006). Two critical parameters needed by numerical models of ice sheet flow are heat flow at the base of glaciers and ice sheets and the material properties of the substrate directly beneath the ice (Schoof & Hewitt 2013). Improved estimates of the subglacial seismic structure can reduce uncertainties in ice flow modelling, which in turn can improve the accuracy of sea level prediction. As described below, both active and passive seismic investigations have been performed to better constrain subglacial seismic structure and to provide insight into crustal and upper mantle structure as well as mechanical properties near the ice-basement interface. The geology and tectonic history beneath Antarctica and Greenland are not well known due to their thick ice cover and remoteness. Early geophysical investigations noted that West Antarctica is more active than the old stable East Antarctica (Behrendt 1999) and associated much of West Antarctica with a rift system that has been undergoing crustal thinning since the Cretaceous (e.g. Behrendt et al. 1991). On the other hand, East Antarctica and the central part of Greenland are considered to be stable consisting of Precambrian cratons (e.g. Moores 1991; Henriksen et al. 2009). Though these general geologic features are well accepted, details of the geological composition and crustal properties beneath polar ice sheets remain unknown. More recent studies (Dahl-Jensen et al. 2003; Winberry & Anandakrishnan 2004; Lawrence et al. 2006; Hansen et al. 2009, 2010; Chaput et al. 2014; Ramirez et al. 2016; Luthra et al. 2016) have provided initial information of subglacial crustal structure, but significant differences between these previous results indicate that more work is needed before our picture of the subsurface is robust. Polar ice sheets pose a challenge to crustal structure imaging. The kilometres-thick ice sheets that mask most of Antarctica and Greenland (Fig. 1) prevent direct observation of geologic features. In addition, seismic waves reverberate strongly within the ice column due to the large acoustic impedance contrast at the base. Those reverberated signals mask the arrivals from crustal structure (e.g. Anandakrishnan & Winberry 2004; Winberry & Anandakrishnan 2004; Wilson & Aster 2005; Heeszel et al. 2016). The teleseismic P-wave receiver function (PRF) techniques used commonly for data from stations installed on bedrock (Vinnik 1977; Langston 1979) yield poor results for data from stations installed on thick ice cover. The strong interference makes construction and assessment of the subice structure using PRFs difficult, leading to large uncertainties (Julià et al. 2004; Hansen et al. 2009). S-wave receiver functions (SRFs) have been used in Antarctica for subsurface structure estimation because the Moho converted phase arrives before the direct S-wave arrival and before the ice-layer reverberations (Hansen et al. 2009, 2010, 2016; Ramirez et al. 2016). However, teleseismic S-wave receiver functions are often noisier than their P-wave counterparts. Methods to improve P-wave PRFs in ice-covered regions can help us better constrain subice crustal structure. Figure 1. Open in new tabDownload slide Station locations (black triangles) and ice thickness maps in Antarctica (left) and Greenland (right). The colours on land show ice thickness variations. In oceans, the colours indicate the bathymetry. Station codes are listed adjacent to the corresponding location triangles. Figure 1. Open in new tabDownload slide Station locations (black triangles) and ice thickness maps in Antarctica (left) and Greenland (right). The colours on land show ice thickness variations. In oceans, the colours indicate the bathymetry. Station codes are listed adjacent to the corresponding location triangles. Ice thickness is well-known in much of Antarctica and Greenland from ice-penetrating radar and reflection seismic measurements (Lythe & Vaughan 2001; Fretwell et al. 2013). The elastic properties of ice have been measured in the laboratory (Thiel & Ostenso 1961), and in situ with active source surveys (Robinson 1968) and passive seismic surveys (Wittlinger & Farra 2015). The result of these efforts is a reasonably accurate and well-established model of the ice sheet’s elastic properties. As shown below, this knowledge we have about ice sheets is sufficient to remove the effects of the ice from the teleseismic signals and to isolate the deeper responses using wavefield continuation and decomposition (WCD) techniques. WCD techniques were developed to remove large amplitude signals generated by near-surface sedimentary structure in order to extract seismic responses from deeper structure (Langston 2011). Using deconvolution, the downward continued signals can be transformed into transfer functions (Langston 2011) and subsurface PRFs (Tao et al. 2014) that are less-contaminated by the reverberations in the shallow structure. The transfer functions or subsurface PRFs can be analysed using modest modifications to the standard receiver function (RF) analysis procedures. Stacking (e.g. Zhu 2000; Zhu & Kanamori 2000), grid search, linearized and stochastic inversion (e.g. Owens et al. 1984; Ammon et al. 1990; Mosegaard & Tarantola 1995; Levin & Park 1997; Sambridge 1999a,b; Nicholson et al. 2005; Shen et al. 2013) have been applied to infer structural information from RF waveforms. Stacking-based techniques use a portion of the waveform and assume simple reference velocity model parametrizations. The exponential increase in the number of simulations needed for each new parameter limits the usage of grid search methods for RF data for complicated models. Linearized inversions that use only the RF data suffer from the possibility that their solution could be trapped in a local minimum (Ammon et al. 1990). Stochastic inversions are gaining in popularity because of improved and less-expensive computational resources. All methods that depend solely on using the receiver functions for the structure inversions suffer from a trade-off between depth and slowness that is difficult to resolve with the range of horizontal slownesses available for teleseismic analyses (Langston 1979; Ammon et al. 1990). With prior information on wave speed, the crustal thickness can be estimated; even with the limitations imposed by the trade-off, the first-order information provided by RFs indicate differences in the subsurface structure among the stations studied. We use data from stations of the Polar Earth Observing Network (POLENET), the Transantarctic Mountains Seismic Experiment (TAMSEIS), the Gamburtsev Mountains Seismic Experiment (GAMSEIS), the Global Seismographic Network (GSN) in Antarctica and the international Greenland Icesheet Seismic Network (GLISN) network in Greenland. These large-scale field experiments have collected data that provide an opportunity to image structure in regions that have heretofore had poor crustal structure estimates. However, in order to fully benefit from these data, new approaches that account for ice-reverberation interference are needed. In the following we apply a subsurface PRF technique (Reading et al. 2003; Langston 2011; Tao et al. 2014) to those data and develop a simplified waveform inversion technique to provide a first-order interpretation of subsurface PRFs. By coupling known ice-layer properties and wave propagation theory, the effects of large amplitude shallow reflections and refractions can be minimized in the teleseismic P-wave signals. We illustrate the procedure using nine stations in Antarctica and one from Greenland. For each, we first compute subsurface PRFs and measure effective subglacial shear velocities. Then we invert the PRFs for first-order crustal structure using a stochastic Markov chain Monte Carlo (McMC) approach. Our goal is to illustrate the procedure and to define a simple inversion tool to guide subsurface PRF interpretation. A sensible longer-term strategy is to combine the subsurface PRFs with complementary information contained in surface-wave dispersion (e.g. Özalaybey et al. 1997; Julià et al. 2000; Chai et al. 2015), S-wave RFs (e.g. Farra & Vinnik 2000; Yuan et al. 2006; Hansen et al. 2009), and Rayleigh wave ellipticity (e.g. Tanimoto & Rivera 2008; Lin et al. 2012). 2 SUBSURFACE RECEIVER FUNCTIONS As described above, near-surface layers with low seismic velocities such as glacial ice and sediments cause large amplitude reverberations. The PRFs calculated from those data have poor resolution of deeper structure. We begin with a numerical illustration of the problem and our solution, comparing the responses of an ice–crust–mantle model and a crust–mantle model to explore this issue. The model parameters are listed in Table 1. The crust–mantle model is the same as the ice–crust–mantle model but with the ice-cover removed. The ice–crust–mantle model has a large ice-bedrock property contrast, in order to illustrate the methodology. The surface and the subsurface PRFs are shown in Fig. 2. As expected, the 2 km thick ice layer produces substantial ice-layer reverberations (see grey curves in Fig. 2a for example ray paths) that complicate and dominate the surface PRFs (waveforms 1 and 2 in Fig. 2b). The dashed lines show the arrival time of the crust–mantle boundary P-to-S conversion (Ps). The converted phase from the crust–mantle boundary (see the dashed line in Fig. 2a for ray path) and multiples from the Moho discontinuity are difficult to identify in both narrowband (Gaussian filter parameter width of 1.0, corresponding to a pulse width of 1.67 s at half the maximum) and broad-band (Gaussian filter parameter width of 5.0, corresponding to a pulse width of 0.75 s at half the maximum) surface PRFs (waveforms 1 and 2 in Fig. 2b). To compute subsurface PRFs from surface seismograms, we downward continue and decompose the vertical and radial motions to obtain upgoing P, downgoing P, upgoing S, and downgoing S wavevectors in the subsurface. Subsurface PRFs are computed by deconvolving the upgoing P wavevector from the upgoing S wavevector using time domain iterative deconvolution (Ligorría & Ammon 1999). The synthetic subsurface PRFs computed using WCD (waveforms 3 and 4 in Fig. 2b) show a clear converted phase and accurate low-frequency multiples from the crust and mantle transition. The Ps arrival time in the subsurface PRFs is the same as that in a surface PRF that is simulated using the model without the ice layer. The test shows that the approach works well for converted phases but that the responses later in the signal (the multiples) include the effects of the ice (multiples are produced at the ice base and at Earth’s surface). This complication is not a problem if the observations and predictions are computed with the same approach and the model includes the relatively well-constrained ice layer. An intuitive view of the procedure is that a virtual station is constructed at the base of the ice layer (the black triangle in Fig. 2a) by downward continuation of the waveforms recorded at the surface (the grey triangle in Fig. 2a). Figure 2. Open in new tabDownload slide Ray paths and synthetic waveforms using the ice–crust–mantle model and the crust–mantle model. (a) Ray paths of the direct P-wave (the thin black line), the converted phase from Moho (the dashed line), and ice-layer reverberations (grey lines) using the ice–crust–mantle model. Thick black lines identify layer boundaries. The grey and black triangles indicate the depth of corresponding waveforms in (b). (b) Synthetic waveforms computed using the ice–crust–mantle model (1)–(4) and the crust–mantle model (5). The waveforms represent (1) narrowband (Gaussian 1.0) surface PRFs; (2) broad-band (Gaussian 5.0) surface PRFs; (3) narrowband subsurface PRFs at the base of the ice layer; (4) broad-band subsurface RFs at the base of the glacier; (5) surface PRFs. Narrowband waveforms (1) and (3) are amplified by 2.8 times to adjust for amplitude differences caused by the Gaussian filters. The two dashed lines indicate the expected arrival time of Moho converted phase at the surface (right line) and at the base of the ice layer (left line). Figure 2. Open in new tabDownload slide Ray paths and synthetic waveforms using the ice–crust–mantle model and the crust–mantle model. (a) Ray paths of the direct P-wave (the thin black line), the converted phase from Moho (the dashed line), and ice-layer reverberations (grey lines) using the ice–crust–mantle model. Thick black lines identify layer boundaries. The grey and black triangles indicate the depth of corresponding waveforms in (b). (b) Synthetic waveforms computed using the ice–crust–mantle model (1)–(4) and the crust–mantle model (5). The waveforms represent (1) narrowband (Gaussian 1.0) surface PRFs; (2) broad-band (Gaussian 5.0) surface PRFs; (3) narrowband subsurface PRFs at the base of the ice layer; (4) broad-band subsurface RFs at the base of the glacier; (5) surface PRFs. Narrowband waveforms (1) and (3) are amplified by 2.8 times to adjust for amplitude differences caused by the Gaussian filters. The two dashed lines indicate the expected arrival time of Moho converted phase at the surface (right line) and at the base of the ice layer (left line). To illustrate the robustness of subsurface PRFs, synthetic surface and subsurface PRFs computed using a series of crust-mantle models and ice–crust–mantle models with varying thickness are shown in Fig. 3. The simple linear arrival-time moveout (dashed lines) expected for changing crustal thickness (Fig. 3b) cannot be seen in surface PRFs for the ice–crust–mantle models (Fig. 3a). Such observable data patterns are valuable for assessing seismic models. Converted phases on the subsurface PRFs recover the expected arrival-time moveout (Figs 3c and d). The ice layer complicates Moho multiples in subsurface PRFs (Fig. 3d, also noted by Tao et al. 2014), but the subsurface PRFs are much simpler than surface PRFs. Considering the synthetic PRFs computed for two different crustal thicknesses (35 km and 40 km) shown in Figs 4(a) and (b), reducing the ice contribution in a PRF increases the importance of signals from the deeper structure and thus enhances sensitivity to the deeper structure. In Fig. 4(c), we illustrate the improved sensitivity using a measure of changes in the PRFs induced by a change in crustal thickness. The vertical axis shows the integral of the squared difference of the PRFs corresponding to a change in crustal thickness. The subsurface PRFs are more sensitive to a crustal thickness perturbation. Wittlinger & Farra (2012, 2015) suggested that first-order variations can exist within the ice, which we have assumed uniform. We investigated two-layer isotopic and anisotropic models based on the work of Wittlinger & Farra (2012, 2015). Details (using isotropic and anisotropic models) are available in the supplement. The calculations suggest that strong variations within the ice will affect the WCD, but we see no evidence (artefacts in the subsurface PRFs) in our subsurface PRFs similar to those that were observed in synthetic PRFs, when an incorrect ice model was used. Figure 3. Open in new tabDownload slide Comparison of surface and subsurface PRF synthetic waveforms for a range of crustal thickness. (a) Narrowband surface PRFs computed using ice–crust–mantle models. The amplitude of all narrowband waveforms is multiplied by 2.8 to account for the Gaussian filter differences. The dashed line shows the expected arrival time for the converted phase from the crust–mantle boundary (we account for the extra traveltime in the ice layer). (b) Broad-band surface PRFs simulated using crust–mantle models. The dashed lines indicate the expected arrival time at the surface of crust–mantle models for the Moho converted phase and multiples. (c) Narrowband subsurface PRFs calculated from ice–crust–mantle models. The amplitude for these narrowband subsurface PRFs is amplified by 2.8 times. (d) Broad-band subsurface PRFs computed from ice–crust–mantle models. The dashed lines in (c) and (d) indicate the expected arrival time at the base of the ice layer of ice–crust–mantle models for the Moho converted phase and multiples. Figure 3. Open in new tabDownload slide Comparison of surface and subsurface PRF synthetic waveforms for a range of crustal thickness. (a) Narrowband surface PRFs computed using ice–crust–mantle models. The amplitude of all narrowband waveforms is multiplied by 2.8 to account for the Gaussian filter differences. The dashed line shows the expected arrival time for the converted phase from the crust–mantle boundary (we account for the extra traveltime in the ice layer). (b) Broad-band surface PRFs simulated using crust–mantle models. The dashed lines indicate the expected arrival time at the surface of crust–mantle models for the Moho converted phase and multiples. (c) Narrowband subsurface PRFs calculated from ice–crust–mantle models. The amplitude for these narrowband subsurface PRFs is amplified by 2.8 times. (d) Broad-band subsurface PRFs computed from ice–crust–mantle models. The dashed lines in (c) and (d) indicate the expected arrival time at the base of the ice layer of ice–crust–mantle models for the Moho converted phase and multiples. Figure 4. Open in new tabDownload slide Synthetic waveform comparison of a 35 km (black lines) and a 45 km (grey lines) crust for surface PRFs (a) and subsurface PRFs (b). (c) Signals differences (per cent difference between a 35 km crust PRF and a PRF computed for a variation in crust thickness normalized by the 35 km crust PRF signal power) as a function of crustal thickness difference respect to 35 km for both surface PRFs and subsurface PRFs. The narrower quadratic for subsurface PRFs indicates the subsurface signals are more sensitivity to crustal thickness changes than surface PRFs. Figure 4. Open in new tabDownload slide Synthetic waveform comparison of a 35 km (black lines) and a 45 km (grey lines) crust for surface PRFs (a) and subsurface PRFs (b). (c) Signals differences (per cent difference between a 35 km crust PRF and a PRF computed for a variation in crust thickness normalized by the 35 km crust PRF signal power) as a function of crustal thickness difference respect to 35 km for both surface PRFs and subsurface PRFs. The narrower quadratic for subsurface PRFs indicates the subsurface signals are more sensitivity to crustal thickness changes than surface PRFs. We performed numerous additional tests using the approach with more complicated (but still known) models. The experiments show that responses caused by the sediments may dominate the subsurface PRFs when shallow interfaces such as thick, low-velocity sediments exist beneath the ice layer. This is no different to computing surface receiver functions using data recorded on a thick, low-speed sedimentary basin. Like the surface receiver functions, the signals can be downward continued through the sedimentary material if that structure is known. But resolving the sedimentary structure using only receiver functions is at least as difficult as it is on a non-ice covered sedimentary basin. 2.1 Application to ice-sheet observations The synthetic tests show that for one-dimensional layered models, WCD allows us to remove shallow P-wave reverberations caused by structure above a reference depth if the seismic structure from the surface to the reference depth is known. To do so, the elastic properties (P-wave and S-wave speeds and density) and layer thickness of all the layers above the reference depth as well as elastic properties just beneath the reference depth are required. For this study, we set the reference depth at the base of the ice layer. Though multiple layers can be included above the reference depth (Supporting Information Fig. S1), a constant-velocity layer is used to represent the ice for simplicity. According to laboratory measurements (Thiel & Ostenso 1961) and active source surveys (Robinson 1968), representative P-wave and S-wave speeds of ice are 3.8 and 1.9 km s−1, respectively. The ice thickness beneath a station is extracted from the BEDMAP2 model (Fretwell et al. 2013) for Antarctic and the ETOPO1 global relief model for Greenland. We also need the average speed at the reference depth for the analysis, which in our case is in the crystalline basement (below the ice and beneath any substantial sedimentary material beneath it). This average basement speed itself has great value for investigating ice-stream or ice sheet basal friction (e.g. Anandakrishnan et al. 1998; Anandakrishnan & Winberry 2004; Leeman et al. 2016). We estimate the average upper-crustal subglacial speed beneath the ice sheet using a grid search for each station (details will be described in following sections). The ice density is fixed at 900 kg m−3. The basement (including sediments if exist) and mantle densities are estimated with empirical relationships from Brocher (2005). Then using standard plane-wave propagator-matrix theory (e.g. Thomson 1950; Haskell 1953) and the expressions in Appendix A, the teleseismic P-waveforms observed on top of the ice are back propagated to the reference depth and the signals are split into upgoing and downgoing P- and S-wave contributions. Basic theory of wavefield downward continuation is reviewed in Appendix A. Our interest is in the wavefield travelling upwards from the deeper structure. A standard P-wave receiver function is the result of the deconvolution of the vertical from the radial component. We compute a subsurface receiver function by deconvolving the upgoing P wavevector from the upgoing S wavevector using time domain iterative deconvolution (Ligorría & Ammon 1999), which is similar to the surface approach (e.g. Langston 1979) that isolates the local response from near source effects. Synthetic subsurface PRFs are computed in the same way but using predicted displacements in place of observed ground motion records. Tao et al. (2014) performed a similar continuation and decomposition and called the results subsurface PRFs and used the PRFs in a multi-stage two-layer H-β searching procedure. We invert the subsurface PRF waveforms instead of searching for minimized upgoing S energy. Waveform inversion works for multiple layers, which can be laborious for the multi-stage H-β stacking. 2.2 Estimating effective subglacial shear wave speed Computing subsurface PRFs require the elastic properties of both the ice layer and the material beneath. Using our estimates of the subsurface PRFs, we are able to measure the shear velocities beneath ice sheets (also referred to as the effective subglacial shear wave velocity). Experiments with one-dimensional numerical models showed that the early part of the subsurface PRF has zero amplitude when the near-true subglacial shear wave speed is used; a positive initial amplitude corresponds to an overestimate of the shear wave speed; and a negative initial amplitude corresponds to an underestimate of the subglacial shear wave speed (Fig. 5a). The same pattern is observed with real data and is illustrated in Fig. 5(b) using observations recorded at seismic station BYRD. We found that we could measure the effect using the energy in the subsurface PRF waveforms prior to 0 s (the Gaussian used in the deconvolution is acausal). This early arrival energy is minimum at the true subglacial shear velocity (Fig. 5c). In Fig. 5(d), the minimum early arrival energy suggests an effective subglacial shear speed of 3.4 km s−1 at station BYRD by measuring early arrival energy of waveforms in Fig. 5(b). Figure 5. Open in new tabDownload slide Subsurface PRFs computed using different subglacial shear velocities and the corresponding early arrival energy measurements for a synthetic model and observations from station BYRD. The model true subglacial shear velocity is 3.4 km s−1. (a) Stacked subsurface PRFs processed using a group of subglacial shear velocities. Each subsurface PRF waveform was averaged with eight subsurface synthetic PRFs using different ray parameters. (b) Stacked subsurface PRFs computed using a set of subglacial shear velocities. Each subsurface PRF was averaged with all the acceptable signals from the analysis of station BYRD. (c) Normalized early arrival (before 0 s) energy measured from waveforms in (a) for different subglacial shear velocities. The early arrival energy is normalized by the maximum energy of the curve; you can see the relative amount of energy in the raw signals. (d) Similar to (c) but using real data from station BYRD. Each sample in (d) was computed from 74 waveforms. Figure 5. Open in new tabDownload slide Subsurface PRFs computed using different subglacial shear velocities and the corresponding early arrival energy measurements for a synthetic model and observations from station BYRD. The model true subglacial shear velocity is 3.4 km s−1. (a) Stacked subsurface PRFs processed using a group of subglacial shear velocities. Each subsurface PRF waveform was averaged with eight subsurface synthetic PRFs using different ray parameters. (b) Stacked subsurface PRFs computed using a set of subglacial shear velocities. Each subsurface PRF was averaged with all the acceptable signals from the analysis of station BYRD. (c) Normalized early arrival (before 0 s) energy measured from waveforms in (a) for different subglacial shear velocities. The early arrival energy is normalized by the maximum energy of the curve; you can see the relative amount of energy in the raw signals. (d) Similar to (c) but using real data from station BYRD. Each sample in (d) was computed from 74 waveforms. The use of the energy in the early part of the PRF helps attack the pragmatic problem of downward continuation. In addition, because the material properties of the layers beneath the ice are also of great interest to ice sheet modellers and those interested in glaciation, we explored our ability to extract information on the structure from the representative values. Unfortunately, heterogeneity (vertical and lateral) in the near surface and the limited signal bandwidth complicates a direct interpretation of our estimated subglacial shear wave speed. Numerical experiments (see the supplement) suggest that the measured subglacial shear wave speed using a signal bandwidth shaped by a Gaussian filter with a width parameter of 1.0 relate to the structure of a varying thickness (up to 10 km beneath the ice) and may represent an average of rock, ice, and sediments (if exist). The exact nature of the average depends on the upper-crustal structure. Using broader-band PRFs can increase resolution, but adds sensitivity to ice structure and lateral variations in the subsurface. Not surprisingly, interpretation of the shallowest structure has to be evaluated on a case-by-case basis. Despite this limitation, the estimated effective subglacial shear velocities are useful for extracting signals originated below the upper crust. 2.3 Results We estimated the subglacial average shear wave speed for the ten stations with thick ice cover. Numerical experiments show that early arrival energy measurements are not sensitive to P-wave speed beneath the ice. We used an empirical relationship (Brocher 2005) between S-wave and P-wave speed to estimate subglacial P-wave velocity. The assumed ice thicknesses and estimated effect on the shear wave speeds are listed in Table 2 and early arrival energy curves are shown in (Fig. 6). Each point in the curves is computed from a few tens to over one hundred P-waveforms. For example, every energy value for station BYRD in Fig. 6 is measured from 74 subsurface PRFs. All early arrival energy curves show a roughly quadratic shape as a function of shear wave speed. Most of the results produce reasonable upper crustal shear wave speeds between 3.0–3.4 km s−1. The subglacial shear speeds are in the range of 3.2–3.6 km s−1 for station BRYD, SUMG, GM02, GM05, N140, N215, and P061. Since these values are comparable to typical upper-crustal speeds (Christensen & Mooney 1995), we infer that these stations are underlain by layers of low-speed material of thickness that is less than a few kilometres. Given our bandwidth and because of the trade-off between ice thickness and near-surface speed we cannot provide better resolution or precision. In contrast, low effective subglacial shear velocities were measured at the South Pole (QSPA) and station E028, (2.9 km s−1 and 3.0 km s−1 respectively), which suggests the existence of relatively thick sediments beneath these sites. Our results compare well with those of Winberry & Anandakrishnan (2004) who imaged a sedimentary layer beneath the ice sheet at the South Pole (station SPA) but no sediments at station BYRD. Station E028 is located above the Wilkes Subglacial Basin (Pyle et al. 2010). Airborne gravity measurements suggest the presence of over 1 km thick sediments across most of the Wilkes Basin (Frederick et al. 2016), which agrees with our subglacial shear speed estimate for the material beneath station E028. The very low effective subglacial shear speed at station SIPL is consistent with the previous investigation by Chaput et al. (2014) and implies significant low-speed material beneath the site. Our estimates of effective subglacial shear speed are consistent with averaged upper crustal speeds from Chaput et al. (2014) for stations BYRD and SIPL with uncertainties on the order of 0.2 km s−1. We did not find reliable upper crustal speeds beneath other stations from previous studies. Figure 6. Open in new tabDownload slide Normalized early arrival energy as a function of subglacial shear velocity for all the stations examined in this study. To facilitate comparison, the arrival energy values were normalized using the maximum energy of each curve (the value differs from curve to curve in the range of 0.007–0.016 s−2 in which station GM02 is the smallest). Each data point in a curve was computed from a minimum of dozens to a maximum of more than one hundred signals. Figure 6. Open in new tabDownload slide Normalized early arrival energy as a function of subglacial shear velocity for all the stations examined in this study. To facilitate comparison, the arrival energy values were normalized using the maximum energy of each curve (the value differs from curve to curve in the range of 0.007–0.016 s−2 in which station GM02 is the smallest). Each data point in a curve was computed from a minimum of dozens to a maximum of more than one hundred signals. Figure 7. Open in new tabDownload slide Surface and subsurface PRFs at station BYRD as a function of ray parameter (steeper incidence angles at the bottom). Grey solid lines are PRFs computed from individual teleseismic events; the black solid lines are averaged PRFs within ray-parameter ranges: 0.06–0.08 s km−1 for the upper black line, 0.05–0.06 s km−1 for the central black line, and 0.04–0.05 s km−1 for the bottom black line. Dashed lines show expected arrival time of different phases using the one-layer velocity model from Winberry & Anandakrishnan (2004). (a) Surface PRFs at station BYRD. Black dashed lines show the arrival times of Ps, PpPs and PpSs+PsPs phases at the surface. Grey dashed lines are for reference. (b) Subsurface PRFs at station BYRD. Black dashed lines show the arrival times of Ps, PpPs and PpSs+PsPs phases at the base of the ice layer. Figure 7. Open in new tabDownload slide Surface and subsurface PRFs at station BYRD as a function of ray parameter (steeper incidence angles at the bottom). Grey solid lines are PRFs computed from individual teleseismic events; the black solid lines are averaged PRFs within ray-parameter ranges: 0.06–0.08 s km−1 for the upper black line, 0.05–0.06 s km−1 for the central black line, and 0.04–0.05 s km−1 for the bottom black line. Dashed lines show expected arrival time of different phases using the one-layer velocity model from Winberry & Anandakrishnan (2004). (a) Surface PRFs at station BYRD. Black dashed lines show the arrival times of Ps, PpPs and PpSs+PsPs phases at the surface. Grey dashed lines are for reference. (b) Subsurface PRFs at station BYRD. Black dashed lines show the arrival times of Ps, PpPs and PpSs+PsPs phases at the base of the ice layer. Table 1. Ice–crust–mantle test model. Description . Vs (km s−1) . Vp (km s−1) . Vp/Vs . Thickness (km) . Ice 1.9 3.8 2.00 2.0 Crust 3.5 6.0 1.71 35.0 Mantle 4.6 8.0 1.74 – Description . Vs (km s−1) . Vp (km s−1) . Vp/Vs . Thickness (km) . Ice 1.9 3.8 2.00 2.0 Crust 3.5 6.0 1.71 35.0 Mantle 4.6 8.0 1.74 – Open in new tab Table 1. Ice–crust–mantle test model. Description . Vs (km s−1) . Vp (km s−1) . Vp/Vs . Thickness (km) . Ice 1.9 3.8 2.00 2.0 Crust 3.5 6.0 1.71 35.0 Mantle 4.6 8.0 1.74 – Description . Vs (km s−1) . Vp (km s−1) . Vp/Vs . Thickness (km) . Ice 1.9 3.8 2.00 2.0 Crust 3.5 6.0 1.71 35.0 Mantle 4.6 8.0 1.74 – Open in new tab Table 2. Summary of results obtained and comparison with previous studies. The Moho Depths listed are 95 per cent confidence intervals and referenced to the base of ice sheets. The quality labels for each station are based on waveform complexity and waveform fits. *Ice thicknesses are from BEDMAP2 (Fretwell et al. 2013). **a. Ramirez et al. (2016); b. Chaput et al. (2014); c. Hansen et al. (2009, 2010); d. Kumar et al. (2007); e. Lawrence et al. (2006); f. Winberry & Anandakrishnan (2004); g. Dahl-Jensen et al. (2003). Station . Latitude . Longitude . Ice* . Subglacial . Moho . Previous studies** . No. of PRFs . Quality . . . . (km) . Vs; Vp . depth . Moho depth . surface; . signal; . . . . . (km s−1) . (km) . (km) . subsurface . result . BYRD −80.02 −119.47 2.2 3.4; 5.8 28 ± 2 24 ± 2b; 27f 85; 74 A; A QSPA −89.93 144.44 2.9 2.9; 4.8 36 ± 4 34f 195; 192 A; A SUMG 72.57 −38.46 3.1 3.3; 5.6 48 ± 4 39d; 49g 132; 113 B; B GM02 −79.43 97.58 2.8 3.3; 5.5 52 ± 3 42 ± 5a; 39 ± 1c 50; 46 B; C GM05 −81.18 51.16 3.5 3.3; 5.5 57 ± 5 49 ± 6a; 47 ± 2c 33; 24 C; C N140 −82.01 96.77 2.8 3.3; 5.6 49 ± 4 46 ± 4a; 46 ± 1c 135; 97 B; C N215 −78.90 59.99 3.5 3.5; 5.9 41 ± 3 49 ± 6a; 45 ± 2c 71; 68 C; C P061 −84.50 77.22 3.2 3.3; 5.6 44 ± 4 43 ± 5a; 43 ± 1c 84; 71 B; B SIPL −81.64 −148.96 1.0 1.8; 3.3 22 ± 3 27 ± 10b 77; 77 B; C E028 −76.31 154.04 1.6 3.1; 5.2 44 ± 20 45 ± 4a; 44 ± 2c; 37e 16; 13 B; D Station . Latitude . Longitude . Ice* . Subglacial . Moho . Previous studies** . No. of PRFs . Quality . . . . (km) . Vs; Vp . depth . Moho depth . surface; . signal; . . . . . (km s−1) . (km) . (km) . subsurface . result . BYRD −80.02 −119.47 2.2 3.4; 5.8 28 ± 2 24 ± 2b; 27f 85; 74 A; A QSPA −89.93 144.44 2.9 2.9; 4.8 36 ± 4 34f 195; 192 A; A SUMG 72.57 −38.46 3.1 3.3; 5.6 48 ± 4 39d; 49g 132; 113 B; B GM02 −79.43 97.58 2.8 3.3; 5.5 52 ± 3 42 ± 5a; 39 ± 1c 50; 46 B; C GM05 −81.18 51.16 3.5 3.3; 5.5 57 ± 5 49 ± 6a; 47 ± 2c 33; 24 C; C N140 −82.01 96.77 2.8 3.3; 5.6 49 ± 4 46 ± 4a; 46 ± 1c 135; 97 B; C N215 −78.90 59.99 3.5 3.5; 5.9 41 ± 3 49 ± 6a; 45 ± 2c 71; 68 C; C P061 −84.50 77.22 3.2 3.3; 5.6 44 ± 4 43 ± 5a; 43 ± 1c 84; 71 B; B SIPL −81.64 −148.96 1.0 1.8; 3.3 22 ± 3 27 ± 10b 77; 77 B; C E028 −76.31 154.04 1.6 3.1; 5.2 44 ± 20 45 ± 4a; 44 ± 2c; 37e 16; 13 B; D Open in new tab Table 2. Summary of results obtained and comparison with previous studies. The Moho Depths listed are 95 per cent confidence intervals and referenced to the base of ice sheets. The quality labels for each station are based on waveform complexity and waveform fits. *Ice thicknesses are from BEDMAP2 (Fretwell et al. 2013). **a. Ramirez et al. (2016); b. Chaput et al. (2014); c. Hansen et al. (2009, 2010); d. Kumar et al. (2007); e. Lawrence et al. (2006); f. Winberry & Anandakrishnan (2004); g. Dahl-Jensen et al. (2003). Station . Latitude . Longitude . Ice* . Subglacial . Moho . Previous studies** . No. of PRFs . Quality . . . . (km) . Vs; Vp . depth . Moho depth . surface; . signal; . . . . . (km s−1) . (km) . (km) . subsurface . result . BYRD −80.02 −119.47 2.2 3.4; 5.8 28 ± 2 24 ± 2b; 27f 85; 74 A; A QSPA −89.93 144.44 2.9 2.9; 4.8 36 ± 4 34f 195; 192 A; A SUMG 72.57 −38.46 3.1 3.3; 5.6 48 ± 4 39d; 49g 132; 113 B; B GM02 −79.43 97.58 2.8 3.3; 5.5 52 ± 3 42 ± 5a; 39 ± 1c 50; 46 B; C GM05 −81.18 51.16 3.5 3.3; 5.5 57 ± 5 49 ± 6a; 47 ± 2c 33; 24 C; C N140 −82.01 96.77 2.8 3.3; 5.6 49 ± 4 46 ± 4a; 46 ± 1c 135; 97 B; C N215 −78.90 59.99 3.5 3.5; 5.9 41 ± 3 49 ± 6a; 45 ± 2c 71; 68 C; C P061 −84.50 77.22 3.2 3.3; 5.6 44 ± 4 43 ± 5a; 43 ± 1c 84; 71 B; B SIPL −81.64 −148.96 1.0 1.8; 3.3 22 ± 3 27 ± 10b 77; 77 B; C E028 −76.31 154.04 1.6 3.1; 5.2 44 ± 20 45 ± 4a; 44 ± 2c; 37e 16; 13 B; D Station . Latitude . Longitude . Ice* . Subglacial . Moho . Previous studies** . No. of PRFs . Quality . . . . (km) . Vs; Vp . depth . Moho depth . surface; . signal; . . . . . (km s−1) . (km) . (km) . subsurface . result . BYRD −80.02 −119.47 2.2 3.4; 5.8 28 ± 2 24 ± 2b; 27f 85; 74 A; A QSPA −89.93 144.44 2.9 2.9; 4.8 36 ± 4 34f 195; 192 A; A SUMG 72.57 −38.46 3.1 3.3; 5.6 48 ± 4 39d; 49g 132; 113 B; B GM02 −79.43 97.58 2.8 3.3; 5.5 52 ± 3 42 ± 5a; 39 ± 1c 50; 46 B; C GM05 −81.18 51.16 3.5 3.3; 5.5 57 ± 5 49 ± 6a; 47 ± 2c 33; 24 C; C N140 −82.01 96.77 2.8 3.3; 5.6 49 ± 4 46 ± 4a; 46 ± 1c 135; 97 B; C N215 −78.90 59.99 3.5 3.5; 5.9 41 ± 3 49 ± 6a; 45 ± 2c 71; 68 C; C P061 −84.50 77.22 3.2 3.3; 5.6 44 ± 4 43 ± 5a; 43 ± 1c 84; 71 B; B SIPL −81.64 −148.96 1.0 1.8; 3.3 22 ± 3 27 ± 10b 77; 77 B; C E028 −76.31 154.04 1.6 3.1; 5.2 44 ± 20 45 ± 4a; 44 ± 2c; 37e 16; 13 B; D Open in new tab 3 SUBSURFACE RECEIVER FUNCTION ANALYSIS Based on experience on processing surface PRFs from the EarthScope Transportable Array (e.g. Chai et al. 2015), a large number of receiver functions obtained even from high-quality stations are quite noisy and some that fit the convolutional model may actually be mostly noise. Waveform selection criteria are necessary to separate signals from noise. We screen the data using surface PRFs. We excluded P-waveforms with a signal-to-noise ratio below ten and surface PRFs with convolution misfits less than 85 per cent of the radial component signal power. Nevertheless, some noisy or contaminated signals may pass through these selection criteria. With a dense network, spatial smoothing can be used to minimize the effects of noisy RFs (Chai et al. 2015). Since the stations in this study were widely separated, we used a clustering analysis to identify and to exclude abnormal waveforms. The clustering analysis arranges surface PRFs into signal-defined groups based on the Euclidean distance between the observed signals. For each station, we selected the largest group for further analysis. Smaller PRF groups usually corresponded to azimuths with few observations. Similar to surface PRF analyses (e.g. Chai et al. 2015), we observe azimuthal dependence in the subsurface PRFs and note that a single ice thickness may not work for groups from different azimuths. However, most of the signals could be used to estimate subsurface PRFs when reasonable estimates of the shallow structure are available or can be calculated. Teleseismic seismograms surviving the surface-PRF screening process were downward continued, decomposed, and deconvolved to estimate subsurface PRFs. Only subsurface PRFs with convolution misfits better than 85 per cent of the upgoing S component signal power were used in further analyses. With one exception, surface and subsurface PRFs were stacked into three ray-parameter bins (0.04–0.05, 0.05–0.06, and 0.06–0.08 s km−1) to further improve the signal-to-noise ratio. Fewer events were available for station E028 due to its relatively short deployment time, so for that station the surface and subsurface PRFs were averaged into single surface and subsurface PRFs. Surface PRFs and subsurface PRFs processed for station BYRD are shown in Figs 7(a) and (b) as a function of ray parameter. This case is particularly interesting since at first look, the surface receiver functions do not look that complicated. However, what one might construe as a Ps conversion from the crust–mantle boundary delayed 5 s, is in fact an interference of ice reverberations. A more likely converted phase is apparent near 2.5–3 s in the subsurface receiver functions. 3.1 Monte Carlo Markov chain inversion Interpretations of RFs are non-unique (e.g. Langston 1979; Ammon et al. 1990; Ammon 1991), but first order information suitable for regional and global comparisons can be extracted using a simple crustal model parametrizations and constraints on the crustal seismic wave speeds (e.g. Zandt & Ammon 1995; Zhu & Kanamori 2000). A range of velocity models could explain PRF observations almost equally well. Formalizing uncertainties for simple models leads naturally to stochastic approaches. We use McMC searches (e.g. Gregory 2010) to quantify the range of seismic structural parameters consistent with the subsurface PRFs and simple one- two-, or three-layer crustal models. As always, the uncertainty estimates are linked to the model parametrization and using only a few layers in the crustal models can lead to subtle differences compared with models including realistic gradients and smooth layer transitions. Uncertainty estimates are derived from the posterior probability distribution, which is computed from the prior information on model parameters and the likelihood function according to Bayes’ theorem. Model parameter priors are truncated to exclude unrealistic models based on global results of active-seismic surveys (Christensen & Mooney 1995) and laboratory measurements (Brocher 2005). Specifically, we remain conservative and restrict crustal shear velocities to the range 2.0–4.5 km s−1 and we also exclude models with crustal P-wave speeds larger than 9.0 km s−1 or smaller than 3.3 km s−1 (these are average speeds in relatively thick crustal layers). Mantle P-wave and S-wave speeds are limited to the range of 7.2–9.6 km s−1 and 4.3–4.8 km s−1, respectively. Based on Brocher’s (2005) compilation, Vp/Vs ratios are limited between 1.53 and 2.00. The crust’s thickness may range between 10 and 75 km. Since we include only a few layers (up to three) in the crust, seismic speeds are forced to increase with depth. We start each search with a simple starting model and then randomly perturb the seismic wave speeds and thickness of up to three crustal layers by sampling from truncated Gaussian prior distributions. Once a perturbed model is generated, subsurface PRFs are computed for both the initial model and the perturbed model. The McMC chain is constructed using a standard Metropolis–Hasting algorithm using a squared-misfit-based likelihood function including a data-derived covariance matrix accounting for correlated noise (e.g. Sambridge 1999a; Agostinetti & Malinverno 2010; Bodin et al. 2012). Details on the likelihood function and the covariance matrix can be found in Appendix B. The choice of the covariance matrix affects the estimated model uncertainty. We compared three types of widely used covariance matrix forms (uniform diagonal, nonuniform diagonal and full matrix) in Appendix B. Not surprisingly, numerical experiments (see Appendix B) suggest that a full covariance matrix can simulate errors in stacked subsurface (and surface) PRFs better than diagonal covariance matrices. Using uniform variance or uncorrelated variances may result in underestimated model parameter uncertainties. Our experiments suggest the uncertainties estimated using a diagonal matrix can be as small as 50 per cent of the uncertainties estimated using more realistic noise statistics (see Appendix B). Full covariance matrix was used in our analysis. For each McMC search, we compute the probability of accepting the newly generated model as the ratio of the likelihood between the new model and the previous model. If the ratio is higher than 1 (i.e. the likelihood of the new model is higher than the likelihood of the previous model), the new model is always accepted. If the ratio is between 0 and 1, the new model will be accepted with such probability (i.e. Lnew/Lpervious). The next perturbation starts with the accepted model. The processes is repeated and the accepted models form a chain. Constructing multiple (5–10) chains with hundreds of thousands of sampled models (about two million models in total for each station) provides a suite of reasonably well-fitting models and an importance sampling of the model-parameter posterior (under the assumption of the model parametrization—no more than three layers and no seismic wave speed inversions). We can marginalize the posterior to estimate the posterior distributions for crustal thickness and average crustal shear speed using the data derived covariance matrix to account for data noise. We summarize the marginal distributions using confidence ellipses calculated using algorithms from Hoover & Rockville (1984). Fig. 8 is a comparison of the McMC search results using synthetic surface and subsurface PRFs. Covariance matrices used in the test surface PRFs and subsurface PRFs were derived from observations at station BYRD. We fixed ice-layer properties in both searches to assure a fair comparison. We can directly compare McMC results of surface PRFs with those from subsurface PRFs. Both posterior distributions are centred on the synthetic target model average speed and crustal thickness. The centre of error ellipses for both surface PRFs and subsurface PRFs are very close to the true values. However, error ellipses for subsurface PRFs are smaller and the marginals for crustal thickness (Fig. 8c) and average crust shear wave speed (Fig. 8d) are narrower for subsurface PRFs. Formally, since uncertainties caused by incorrect ice-layer properties are not included in the McMC calculation, the uncertainties computed from McMC do not include the effects of ice property uncertainty. We used synthetic tests to explore the impact of using incorrect ice-layer properties to compute subsurface PRFs (see Supporting Information Fig. S7). For modest variations in the ice properties (a 5 per cent, 0.1 km error in ice thickness or a 5 per cent, 0.1 km s−1 error in ice shear velocity) the McMC results converge to the correct range of model parameters. Incorrect ice-layer properties shift the error ellipses a distance smaller than the original estimated uncertainties. Large errors in the assumed ice structure produce notable artefacts in the PRFs and thus can be identified as part of the data processing. Figure 8. Open in new tabDownload slide McMC search results using synthetic surface PRFs and subsurface PRFs using data-derived covariances compute for station BYRD. The velocity model can be found in Table S4 in the electronic supplement to this article. Model parameters of the ice layer were fixed. (a) Crustal thickness of all the velocity models constrained by surface PRFs as a function of averaged crust shear wave speed, Vs. Dashed lines indicate true values. The ellipses are error ellipses represent 95 per cent confidence interval (outer) and 68 per cent confidence interval (inner). (b) same as (a), but the velocity models were constrained by subsurface PRFs. (c) Prior and marginal distributions of crustal thickness. (d) Prior and marginal distributions of average crust shear wave speed. The prior distributions in (c) and (d) extend beyond the plot limits. Dashed lines indicate true values. Figure 8. Open in new tabDownload slide McMC search results using synthetic surface PRFs and subsurface PRFs using data-derived covariances compute for station BYRD. The velocity model can be found in Table S4 in the electronic supplement to this article. Model parameters of the ice layer were fixed. (a) Crustal thickness of all the velocity models constrained by surface PRFs as a function of averaged crust shear wave speed, Vs. Dashed lines indicate true values. The ellipses are error ellipses represent 95 per cent confidence interval (outer) and 68 per cent confidence interval (inner). (b) same as (a), but the velocity models were constrained by subsurface PRFs. (c) Prior and marginal distributions of crustal thickness. (d) Prior and marginal distributions of average crust shear wave speed. The prior distributions in (c) and (d) extend beyond the plot limits. Dashed lines indicate true values. 3.2 Results We show waveforms and velocity models estimated using the observations from station BYRD in Fig. 9 as an example (complete results for the other nine stations are documented in the electronic supplement). For BYRD, better fitting models (higher probability, darker dots in Fig. 9a) show that average crustal shear wave speed correlates negatively with crustal thickness, which is a pattern observed for all other stations. This observation can be attributed to the velocity-thickness trade-off intrinsic to receiver functions with a limited ray-parameter range. Sample predicted waveforms associated with the average of all models in the McMC chain are shown in Fig. 9(b). Since the model-data relationship is nonlinear, averaging models is not guaranteed to produce a good-fitting model, but in this case, the model fits comparably to the better fitting models in the chain. To account for layer thickness differences, the models were sampled at 1 km intervals during averaging. The result (Fig. 9c) is a higher-dimensional model than those included in the search, but the consistency of the better fitting models in the chain reduces complexity associate with this mapping. The average model includes a mid-crust transition and a relative sharp crust-mantle transition. Although surface PRFs were not included in the inversion, we computed synthetic surface PRFs using the averaged model (Fig. 9d) as an additional consistency check. The predicted surface PRF waveforms agree well with the stacked surface PRFs. The same is not true at all stations. Figure 9. Open in new tabDownload slide McMC search results for signals from station BYRD. (a) Crustal thickness of all the velocity models as a function of averaged crust shear wave speed, Vs. The ellipses are error ellipses represent 95 per cent confidence interval (outer) and 68 per cent confidence interval (inner). Each dot represents one velocity model. Darker dots have larger likelihood than lighter ones. Top and right panels show the marginal distributions. (b) Predicted subsurface PRFs (black lines) of the McMC-averaged model in (c) and stacked subsurface PRF observations (grey lines, same as in Fig. 7a). (c) The Shear velocities profile of the McMC-averaged model is shown with the black line. Grey lines indicate two standard derivations (95 per cent confidence interval). (d) Predicted surface PRFs (black lines) of the McMC-averaged model and stacked surface PRF observations (grey lines, same as in Fig. 7b) (surface PRFs were not used in the McMC analysis). Figure 9. Open in new tabDownload slide McMC search results for signals from station BYRD. (a) Crustal thickness of all the velocity models as a function of averaged crust shear wave speed, Vs. The ellipses are error ellipses represent 95 per cent confidence interval (outer) and 68 per cent confidence interval (inner). Each dot represents one velocity model. Darker dots have larger likelihood than lighter ones. Top and right panels show the marginal distributions. (b) Predicted subsurface PRFs (black lines) of the McMC-averaged model in (c) and stacked subsurface PRF observations (grey lines, same as in Fig. 7a). (c) The Shear velocities profile of the McMC-averaged model is shown with the black line. Grey lines indicate two standard derivations (95 per cent confidence interval). (d) Predicted surface PRFs (black lines) of the McMC-averaged model and stacked surface PRF observations (grey lines, same as in Fig. 7b) (surface PRFs were not used in the McMC analysis). For the station E028, synthetic surface PRFs do not agree with the surface PRF observations well and we suspect that the differences are caused by sedimentary layers that are not well modelled using thick layers and subsurface PRFs (see Supporting Information Fig. S16). The waveform fits for both surface and subsurface PRFs are also not optimal at station SIPL suggesting significant heterogeneity in the crust and possibly in the mantle as well (see Supporting Information Fig. S15). At station SIPL, we observed secondary peaks in marginal distribution that are possibly caused by cycle skipping that correspond to unlikely models (extreme average speeds and thicknesses) and are easy to eliminate on a geologic basis. Crustal thickness and related uncertainties are listed in Table 2. In general, our results compare well with those results from other studies, including independent S-wave receiver function analyses. At station SUMG, we estimate a crust–mantle transition depth between 44 and 52 km, which is typical for a stable Precambrian shield. The result agrees well with the estimate from Dahl-Jensen et al. (2003), but is thicker than the value of Kumar et al. (2007). For the stations GM02 and GM05, our estimates of Moho depths are deeper than previous studies. The crustal thickness estimates for most of the stations have relatively large uncertainties (∼±4 km) due to the correlation (trade-off) between the layer’s seismic velocity and layer thickness. A exception is station E028, the uncertainty in Moho depth is 20 km due to limited data. 4 DISCUSSION AND CONCLUSIONS Our primary goal in this work is to explore and to demonstrate the potential for subsurface receiver function analysis in regions of thick ice cover. Using ten examples we have shown how subsurface PRF analysis can isolate the response from the lower crust and upper mantle while removing the large-amplitude reverberations from shallow low-velocity ice layers. Subsurface PRF arrival times of converted phase and multiples are the same as those of a surface PRF recorded at a virtual station located at the reference depth. Although the amplitude of the subsurface PRF is different from the surface PRF, the subsurface PRF is affected by the structural parameters in the same manner as that of the surface PRF. In general, a subsurface PRF can be treated as a scaled surface PRF recorded at a virtual station at a reference depth with initial P-wave removed and some added complexity associated with the ice cover in the multiples. The effective subglacial shear velocity measured by minimizing early arrival energy on subsurface PRFs is an indicator of subglacial upper-crustal structure and can identify stations with thick near-surface sedimentary structures. Our effective subglacial shear velocity measurements agree with previous estimates at the same locales. Unfortunately, but not surprisingly, due to limited bandwidth, subsurface PRFs alone cannot uniquely determine detailed structure of subglacial sediments (the same is true for surface observations without ice). However, our numerical experiments and the examples from ten seismic stations located on ice sheets suggest the effective subglacial shear velocity measurement contribute some useful information on the subglacial near-surface structure. We used a simple McMC search using subsurface PRFs to constrain simple models of the crust and to quantify the associated uncertainties. Although using a few constant-velocity layers is a rough approximation of the often-complex subsurface geology, the structural estimations derived from subsurface PRFs beneath ice-covered stations provide first-order information on subglacial structure. Most of our estimates from the ten example stations corroborate previous findings in this regard. Integration of the subsurface PRFs with complementary information such as surface PRFs, shear wave receiver functions, surface-wave dispersion observation and ellipticity measurements can reduce uncertainties in structural investigations. With no new information added to the process, careful investigations of the subsurface structure using a fixed ice layer have been able to resolve some of the subsurface features (e.g. Anandakrishnan & Winberry 2004). Isolating the response of the deeper structure with a subsurface PRF allows the analyst to see the signals that constrain the deeper structure and better assess the appropriateness of a particular model of the deeper structure. Therefore, the approach has the potential to greatly ease the analysis of many well-recorded P-wave observations that have been difficult to interpret using standard surface-receiver functions approaches. Acknowledgments We thank editor Lapo Boschi and two anonymous reviewers for their constructive criticism to greatly improve the manuscript. The U.S. National Science Foundation supported this work through the grant EAR-1053484. The NSF Office of Polar Programs grants 0632230, 0632239, 0652322, 0632335, 0632136, 0632209 and 0632185 support the POLENET-NAET program that collected most of the data used in this work. We thank the Pennsylvania State University for providing computing resources though Lion-X clusters that were used for our data processing. The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to waveforms, related metadata, and/or derived products used in this study. IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope (SAGE) Proposal of the National Science Foundation under Cooperative Agreement EAR-1261681. We used seismic data from the Global Seismographic Network (doi:10.7914/SN/IU), IPY POLENET-Antarctica (doi:10.7914/SN/YT_2007), a broad-band seismic experiment to image the lithosphere beneath the Ganburtsev Mountains (doi:10.7914/SN/ZM_2007), a broad-band seismic investigation of deep continental structure across the East-West Antarctic Boundary (TAMSEIS, doi:10.7914/SN/XP_2000), and the GEOFON network (doi:10.14470/TR560404). The ETOPO1 Global Relief Model (https://www.ngdc.noaa.gov, doi:10.7289/V5C8276M, last accessed October 2012) was used to extract ice thicknesses in Greenland. Ice thicknesses in Antarctica were extracted using the Bedmap2 Toolbox for Matlab version 4.2 (http://www.mathworks.com/matlabcentral/fileexchange/42353-bedmap2-toolbox-for-matlab, last accessed September 2015) and the Bedmap2 ice thickness model (https://www.bas.ac.uk/project/bedmap-2, (Fretwell et al. 2013)). Some plots were made using the Generic Mapping Tools version 5.2.1 (http://gmt.soest.hawaii.edu, Wessel et al. 2013) and Matplotlib version 1.5.1 (Hunter 2007). Obspy version 1.0 (Beyreuther et al. 2010) and Numpy version 1.10.4 (van der Walt et al. 2011) were used to process data. The h5py version 2.5.0 was used to store models and data (http://www.h5py.org, last accessed May 2016). RFSYN package was used to compute synthetic waveforms for anisotropic models (Levin & Park 1998). REFERENCES Agostinetti N.P. , Malinverno A., 2010 . Receiver function inversion by trans-dimensional Monte Carlo sampling , Geophys. J. Int. , 181 ( 2 ), 858 – 872 . OpenURL Placeholder Text WorldCat Ammon C.J. , 1991 . The isolation of receiver effects from teleseismic P waveforms , Bull. seism. Soc. Am. , 81 ( 6 ), 2504 – 2510 . OpenURL Placeholder Text WorldCat Ammon C.J. , Randall G.E., Zandt G., 1990 . On the nonuniqueness of receiver function inversions , J. geophys. Res. , 95 ( B10 ), 15 303 – 15 318 . Google Scholar Crossref Search ADS WorldCat Anandakrishnan S. , Winberry J.P., 2004 . Antarctic subglacial sedimentary layer thickness from receiver function analysis , Glob. Planet. Change , 42 ( 1–4 ), 167 – 176 . Google Scholar Crossref Search ADS WorldCat Anandakrishnan S. , Blankenship D.D., Alley R.B., Stoffa P.L., 1998 . Influence of subglacial geology on the position of a West Antarctic ice stream from seismic observations , Nature , 394 ( 6688 ), 62 – 65 . Google Scholar Crossref Search ADS WorldCat Behrendt J. , 1999 . Crustal and lithospheric structure of the West Antarctic Rift System from geophysical investigations—a review , Glob. Planet. Change , 23 ( 1–4 ), 25 – 44 . Google Scholar Crossref Search ADS WorldCat Behrendt J.C. , LeMasurier W.E., Cooper A.K., Tessensohn F., Tréhu A., Damaske D., 1991 . Geophysical studies of the West Antarctic Rift System , Tectonics , 10 ( 6 ), 1257 – 1273 . Google Scholar Crossref Search ADS WorldCat Beyreuther M. , Barsch R., Krischer L., Megies T., Behr Y., Wassermann J., 2010 . ObsPy: A Python toolbox for seismology , Seismol. Res. Lett. , 81 ( 3 ), 530 – 533 . Google Scholar Crossref Search ADS WorldCat Bodin T. , Sambridge M., Tkalčić H., Arroucau P., Gallagher K., Rawlinson N., 2012 . Transdimensional inversion of receiver functions and surface wave dispersion , J. geophys. Res. , 117, B02301, doi:10.1029/2011JB008560 . OpenURL Placeholder Text WorldCat Brocher T.M. , 2005 . Empirical relations between elastic wavespeeds and density in the Earth’s crust , Bull. seism. Soc. Am. , 95 ( 6 ), 2081 – 2092 . Google Scholar Crossref Search ADS WorldCat Chai C. , Ammon C.J., Maceira M., Herrmann R.B., 2015 . Inverting interpolated receiver functions with surface wave dispersion and gravity: application to the western U.S. and adjacent Canada and Mexico , Geophys. Res. Lett. , 42 ( 11 ), 4359 – 4366 . Google Scholar Crossref Search ADS WorldCat Chaput J. et al. , 2014 . The crustal thickness of West Antarctica , J. geophys. Res. , 119 ( 1 ), 378 – 395 . Google Scholar Crossref Search ADS WorldCat Christensen N.I. , Mooney W.D., 1995 . Seismic velocity structure and composition of the continental crust: a global view , J. geophys. Res. , 100 ( B6 ), 9761 – 9788 . Google Scholar Crossref Search ADS WorldCat Dahl-Jensen T. et al. , 2003 . Depth to Moho in Greenland: receiver-function analysis suggests two Proterozoic blocks in Greenland , Earth planet. Sci. Lett. , 205 ( 3–4 ), 379 – 393 . Google Scholar Crossref Search ADS WorldCat Farra V. , Vinnik L., 2000 . Upper mantle stratification by P and S receiver functions , Geophys. J. Int. , 141 ( 3 ), 699 – 712 . Google Scholar Crossref Search ADS WorldCat Frederick B.C. , Young D.A., Blankenship D.D., Richter T.G., Kempf S.D., Ferraccioli F., Siegert M.J., 2016 . Distribution of subglacial sediments across the Wilkes Subglacial Basin, East Antarctica , J. geophys. Res. , 121 ( 4 ), 790 – 813 . Google Scholar Crossref Search ADS WorldCat Fretwell P. et al. , 2013 . Bedmap2: improved ice bed, surface and thickness datasets for Antarctica , Cryosphere , 7 ( 1 ), 375 – 393 . Google Scholar Crossref Search ADS WorldCat Gregory P. , 2010 . Bayesian Logical Data Analysis for the Physical Sciences: A comparative approach with Mathematica® Support, Cambridge Univ. Press . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Gregory J.M. , Huybrechts P., 2006 . Ice-sheet contributions to future sea-level change , Phil. Trans. R. Soc. A , 364 ( 1844 ), 1709 – 1732 . Google Scholar Crossref Search ADS WorldCat Hansen S.E. , Julià J., Nyblade A.A., Pyle M.L., Wiens D.A., Anandakrishnan S., 2009 . Using S wave receiver functions to estimate crustal structure beneath ice sheets: an application to the Transantarctic Mountains and East Antarctic craton , Geochem. Geophys. Geosyst. , 10 , Q08014 , doi:10.1029/2009GC002576 . Google Scholar Crossref Search ADS WorldCat Hansen S.E. , Nyblade A.A., Heeszel D.S., Wiens D.A., Shore P., Kanao M., 2010 . Crustal structure of the Gamburtsev Mountains, East Antarctica, from S-wave receiver functions and Rayleigh wave phase velocities , Earth planet. Sci. Lett. , 300 ( 3–4 ), 395 – 401 . Google Scholar Crossref Search ADS WorldCat Hansen S.E. , Kenyon L.M., Graw J.H., Park Y., Nyblade A.A., 2016 . Crustal structure beneath the Northern Transantarctic Mountains and Wilkes Subglacial Basin: implications for tectonic origins , J. geophys. Res. , 121 ( 2 ), 812 – 825 . Google Scholar Crossref Search ADS WorldCat Haskell N.A. , 1953 . The dispersion of surface waves on multilayered media , Bull. seism. Soc. Am. , 43 ( 1 ), 17 – 34 . OpenURL Placeholder Text WorldCat Heeszel D.S. et al. , 2016 . Upper mantle structure of central and West Antarctica from array analysis of Rayleigh wave phase velocities , J. geophys. Res. , 121 ( 3 ), 1758 – 1775 . Google Scholar Crossref Search ADS WorldCat Henriksen N. , Higgins A.K., Kalsbeek F., Pulvertaft T. C.R., 2009 . Greenland from Archaean to Quaternary: descriptive text to the 1995 geological map of Greenland, 1:2500000, 2nd edition , Geol. Surv. Denmark and Greenland Bull. , 18 , 1 – 126 . OpenURL Placeholder Text WorldCat Hoover W.E. , Rockville M.D., 1984 . Algorithms for confidence circles and ellipses, Tech. rep., U. S. Department of Commerce , National Oceanic and Atmospheric Administration , NOAA Technical Report NOS 107 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hunter J.D. , 2007 . Matplotlib: A 2D graphics environment , Comput, Sci. Eng. , 9 ( 3 ), 90 – 95 . Google Scholar Crossref Search ADS WorldCat Huybrechts P. , de Wolde J., 1999 . The dynamic response of the Greenland and Antarctic ice sheets to multiple-century climatic warming , J. Climate , 12 ( 8 ), 2169 – 2188 . Google Scholar Crossref Search ADS WorldCat Julià J. , Ammon C.J., Herrmann R.B., Correig A.M., 2000 . Joint inversion of receiver function and surface wave dispersion observations , Geophys. J. Int. , 143 ( 1 ), 99 – 112 . Google Scholar Crossref Search ADS WorldCat Julià J. , Herrmann R.B., Ammon C.J., Akinci A., 2004 . Evaluation of deep sediment velocity structure in the New Madrid Seismic Zone , Bull. seism. Soc. Am. , 94 ( 1 ), 334 – 340 . Google Scholar Crossref Search ADS WorldCat Kumar P. , Kind R., Priestley K., Dahl-Jensen T., 2007 . Crustal structure of Iceland and Greenland from receiver function studies , J. geophys. Res. , 112 , B03301 , doi:10.1029/2005JB003991. OpenURL Placeholder Text WorldCat Langston C.A. , 1979 . Structure under Mount Rainier, Washington, inferred from teleseismic body waves , J. geophys. Res. , 84 ( B9 ), 4749 – 4762 . Google Scholar Crossref Search ADS WorldCat Langston C.A. , 2011 . Wave-Field continuation and decomposition for passive seismic imaging under deep unconsolidated sediments , Bull. seism. Soc. Am. , 101 ( 5 ), 2176 – 2190 . Google Scholar Crossref Search ADS WorldCat Lawrence J.F. , Wiens D.A., Nyblade A.A., Anandakrishnan S., Shore P.J., Voigt D., 2006 . Crust and upper mantle structure of the Transantarctic Mountains and surrounding regions from receiver functions, surface waves, and gravity: implications for uplift models , Geochem. Geophys. Geosyst. , 7 ( 10 ), 1 – 23 . OpenURL Placeholder Text WorldCat Leeman J.R. , Valdez R.D., Alley R.B., Anandakrishnan S., Saffer D.M., 2016 . Mechanical and hydrologic properties of Whillans Ice Stream till: implications for basal strength and stick-slip failure , J. geophys. Res. , 121 ( 7 ), 1295 – 1309 . Google Scholar Crossref Search ADS WorldCat Levin V. , Park J., 1997 . P–SH conversions in a flat-layered medium with anisotropy of arbitrary orientation , Geophys. J. Int. , 131 ( 2 ), 253 – 266 . Google Scholar Crossref Search ADS WorldCat Levin V. , Park J., 1998 . P–SH Conversions in layered media with hexagonally symmetric anisotropy: a CookBook , Pure appl. Geophys. , 151 ( 2-4 ), 669 – 697 . Google Scholar Crossref Search ADS WorldCat Ligorría J.P. , Ammon C.J., 1999 . Iterative deconvolution and receiver-function estimation , Bull. seism. Soc. Am. , 89 ( 5 ), 1395 – 1400 . OpenURL Placeholder Text WorldCat Lin F.-C. , Schmandt B., Tsai V.C., 2012 . Joint inversion of Rayleigh wave phase velocity and ellipticity using USArray: constraining velocity and density structure in the upper crust , Geophys. Res. Lett. , 39 ( 12 ), 1 – 7 . Google Scholar Crossref Search ADS WorldCat Luthra T. , Anandakrishnan S., Winberry J.P., Alley R.B., Holschuh N., 2016 . Basal characteristics of the main sticky spot on the ice plain of Whillans Ice Stream, Antarctica , Earth planet. Sci. Lett. , 440 , 12 – 19 . Google Scholar Crossref Search ADS WorldCat Lythe M.B. , Vaughan D.G., 2001 . BEDMAP: A new ice thickness and subglacial topographic model of Antarctica , J. geophys. Res. , 106 ( B6 ), 11 335 – 11 351 . Google Scholar Crossref Search ADS WorldCat Moores E.M. , 1991 . Southwest U.S.-East Antarctic (SWEAT) connection: a hypothesis , Geology , 19 ( 5 ), 425 – 428 . Google Scholar Crossref Search ADS WorldCat Mosegaard K. , Tarantola A., 1995 . Monte Carlo sampling of solutions to inverse problems , J. geophys. Res. , 100 ( B7 ), 12 431 – 12 447 . Google Scholar Crossref Search ADS WorldCat Nicholson T. , Bostock M., Cassidy J.F., 2005 . New constraints on subduction zone structure in northern Cascadia , Geophys. J. Int. , 161 ( 3 ), 849 – 859 . Google Scholar Crossref Search ADS WorldCat Owens T.J. , Zandt G., Taylor S.R., 1984 . Seismic evidence for an ancient rift beneath the Cumberland Plateau, Tennessee: a detailed analysis of broadband teleseismic P waveforms , J. geophys. Res. , 89 ( B9 ), 7783 – 7795 . Google Scholar Crossref Search ADS WorldCat Özalaybey S. , Savage M.K., Sheehan A.F., Louie J.N., Brune J.N., 1997 . Shear wave velocity structure in the northern Basin and Range province from the combined analysis of receiver functions and surface waves , Bull. seism. Soc. Am. , 87 ( 1 ), 183 – 199 . OpenURL Placeholder Text WorldCat Pyle M.L. , Wiens D.A., Nyblade A.A., Anandakrishnan S., 2010 . Crustal structure of the Transantarctic Mountains near the Ross Sea from ambient seismic noise tomography , J. geophys. Res. , 115 , B11310 , doi:10.1029/2009JB007081 . Google Scholar Crossref Search ADS WorldCat Ramirez C. et al. , 2016 . Crustal and upper-mantle structure beneath ice-covered regions in Antarctica from S-wave receiver functions and implications for heat flow , Geophys. J. Int. , 204 ( 3 ), 1636 – 1648 . Google Scholar Crossref Search ADS WorldCat Reading A. , Kennett B., Sambridge M., 2003 . Improved inversion for seismic structure using transformed, S-wavevector receiver functions: removing the effect of the free surface , Geophys. Res. Lett. , 30 ( 19 ), 1981 , doi:10.1029/2003GL018090 . Google Scholar Crossref Search ADS WorldCat Robinson E.S. , 1968 . Seismic Wave Propagation on a Heterogeneous Polar Ice Sheet , J. geophys. Res. , 73 ( 2 ), 739 – 753 . Google Scholar Crossref Search ADS WorldCat Sambridge M. , 1999a . Geophysical inversion with a neighbourhood algorithm—II. Appraising the ensemble , Geophys. J. Int. , 138 ( 3 ), 727 – 746 . Google Scholar Crossref Search ADS WorldCat Sambridge M. , 1999b . Geophysical inversion with a neighbourhood algorithm—I. Searching a parameter space , Geophys. J. Int. , 138 ( 2 ), 479 – 494 . Google Scholar Crossref Search ADS WorldCat Schoof C. , Hewitt I., 2013 . Ice-Sheet dynamics , Annual Review of Fluid Mechanics , 45 ( 1 ), 217 – 239 . Google Scholar Crossref Search ADS WorldCat Shen W. , Ritzwoller M.H., Schulte Pelkum V., Lin F.-C., 2013 . Joint inversion of surface wave dispersion and receiver functions: a Bayesian Monte-Carlo approach , Geophys. J. Int. , 192 ( 2 ), 807 – 836 . Google Scholar Crossref Search ADS WorldCat Tanimoto T. , Rivera L., 2008 . The ZH ratio method for long-period seismic data: sensitivity kernels and observational techniques , Geophys. J. Int. , 172 ( 1 ), 187 – 198 . Google Scholar Crossref Search ADS WorldCat Tao K. , Liu T., Ning J., Niu F., 2014 . Estimating sedimentary and crustal structure using wavefield continuation: theory, techniques and applications , Geophys. J. Int. , 197 ( 1 ), 443 – 457 . Google Scholar Crossref Search ADS WorldCat Thiel E. , Ostenso N.A., 1961 . Seismic studies on Antarctic ice shelves , Geophysics , 26 ( 6 ), 706 – 715 . Google Scholar Crossref Search ADS WorldCat Thomas R. et al. , 2004 . Accelerated sea-level rise from West Antarctica , Science , 306 ( 5694 ), 255 – 258 . Google Scholar Crossref Search ADS PubMed WorldCat Thomson W.T. , 1950 . Transmission of elastic waves through a stratified solid medium , J. appl. Phys. , 21 ( 2 ), 89 – 93 . Google Scholar Crossref Search ADS WorldCat van der Walt S. , Colbert S.C., Varoquaux G., 2011 . The NumPy Array: A structure for efficient numerical computation , Comput. Sci. Eng. , 13 ( 2 ), 22 – 30 . Google Scholar Crossref Search ADS WorldCat Vinnik L.P. , 1977 . Detection of waves converted from P to SV in the mantle , Phys. Earth planet. Inter. , 15 ( 1 ), 39 – 45 . Google Scholar Crossref Search ADS WorldCat Wessel P. , Smith W. H.F., Scharroo R., Luis J., Wobbe F., 2013 . Generic Mapping Tools: improved version released , EOS, Trans. Am. geophys. Un. , 94 ( 45 ), 409 – 410 . Google Scholar Crossref Search ADS WorldCat Wilson D. , Aster R., 2005 . Seismic imaging of the crust and upper mantle using regularized joint receiver functions, frequency–wave number filtering, and multimode Kirchhoff migration , J. geophys. Res. , 110 , B05305 , doi:10.1029/2004JB003430 . OpenURL Placeholder Text WorldCat Winberry J.P. , Anandakrishnan S., 2004 . Crustal structure of the West Antarctic rift system and Marie Byrd Land hotspot , Geology , 32 ( 11 ), 977 – 980 . Google Scholar Crossref Search ADS WorldCat Wittlinger G. , Farra V., 2012 . Observation of low shear wave velocity at the base of the polar ice sheets: evidence for enhanced anisotropy , Geophys. J. Int. , 190 ( 1 ), 391 – 405 . Google Scholar Crossref Search ADS WorldCat Wittlinger G. , Farra V., 2015 . Evidence of unfrozen liquids and seismic anisotropy at the base of the polar ice sheets , Polar Sci. , 9 ( 1 ), 66 – 79 . Google Scholar Crossref Search ADS WorldCat Yuan X.H. , Kind R., Li X.Q., Wang R.J., 2006 . The S receiver functions: synthetics and data example , Geophys. J. Int. , 165 ( 2 ), 555 – 564 . Google Scholar Crossref Search ADS WorldCat Zandt G. , Ammon C.J., 1995 . Continental crust composition constrained by measurements of crustal Poisson’s ratio , Nature , 374 ( 6518 ), 152 – 154 . Google Scholar Crossref Search ADS WorldCat Zhu L. , 2000 . Crustal structure across the San Andreas Fault, southern California from teleseismic converted waves , Earth planet. Sci. Lett. , 179 ( 1 ), 183 – 190 . Google Scholar Crossref Search ADS WorldCat Zhu L.P. , Kanamori H., 2000 . Moho depth variation in southern California from teleseismic receiver functions , J. geophys. Res. , 105 ( B2 ), 2969 – 2980 . Google Scholar Crossref Search ADS WorldCat SUPPORTING INFORMATION Supplementary data are available at GJIRAS online. Figure S1. Synthetic waveforms of subsurface PRF processed using incorrect ice thicknesses (a) and incorrect ice shear velocities (b) in the reference model. The Vp/Vs ratio is fixed for (b). The velocity models are perturbed based on the Ice-Crust-Mantle model (Table 1) in the paper. Results for the correct ice model are shown as the middle trace in each panel. Figure S2. Comparison of surface and subsurface PRF synthetic waveforms for a range of crustal thickness using two layers in the ice. (a) Surface PRFs computed using the model S1 in Table S1. (b) Subsurface PRFs downward continued to the base of the first layer of the ice using the model S1. (c) Subsurface PRFs downward continued to the base of the second layer of the ice using the model S1. (d) Same as (c) but using the model S2 in Table S2 for downward continuation. Figure S3. Comparison of surface (a, b) and subsurface (c, d) PRF synthetic waveforms for isotropic (model S1) and anisotropic ice models (model S1 with 4 per cent Vp anisotropy and −6 per cent Vs anisotropy in the second ice layer with the c-axis in the vertical direction) using two bandwidths Gaussian 1.0 (a, c) and Gaussian 5.0 (b, d). The anisotropic subsurface PRFs used synthetic waveforms from the anisotropic ice model but downward continued with the isotropic model S1. Figure S4. Effective subglacial shear velocities measured from synthetic waveforms using model S3 in Table S3 with a series of different upper crust thickness (a) and shear velocity (b). A Gaussian filter with a width parameter of 1.0 was used. Numbers adjacent to curves in (a) represent upper crust layer thickness of the true model. Numbers adjacent to curves in (b) show upper crust shear velocity of the true model. Gray lines indicate the effective subglacial shear velocities that have minimum early arrival energy. The upper crust shear velocity is fixed at 2.8 km s−1 in (a). The upper crust layer thickness in (b) is fixed at 10 km. The estimated speed depends on the shallow structure, but represents a vertical average of the material properties. Figure S5. Similar to S3, but use Gaussian filter with a width parameter of 5.0. The upper crust shear velocity is fixed at 2.8 km s−1 in (a). The upper crust layer thickness in (b) is fixed at 3 km. Figure S6. Prior distributions of crustal thickness and average crust shear-wave velocity. Gray dots show crustal thickness of all the sampled velocity models as a function of averaged crust shear-wave speed. Top and right panels show the prior distributions. Figure S7. McMC search results processed using incorrect ice-layer properties. The velocity models are changed based on Model S4 (Table S4). The format of the figure is similar to Fig. 8(b). (a) The reference model has an ice thickness of 0.1 km smaller than the correct value. (b) The reference model has an ice thickness of 0.1 km larger than the correct value. (c) The reference model has an ice shear-wave speed of 0.1 km s−1 smaller than the correct value. (d) The reference model has an ice shear-wave speed of 0.1 km s−1 larger than the correct value. Ice thickness uncertainty estimates are often less than 0.1 km. Figure S8. McMC search results for signals from station QSPA. See supplement introduction for layout description. We consider these results quality A. The predictions are slightly higher frequency than the observations because of the limited number of model parameters, but generally fit the main features in the observations well. The models are relatively consistent and the surface receiver function fits are acceptable. Figure S9. McMC search results for signals from station SUMG. We consider these results quality B primarily because the subsurface receiver functions are relatively featureless. The signals suggest minimal sharp seismic boundaries beneath the site including a gradual crust-mantle transition. The surface receiver function fits are quite good, but the signals come from almost all ice. Figure S10. McMC search results for signals from station GM02. See supplement introduction for layout description. We consider these results quality C since our fit to the subsurface receiver functions is poor for the initial part of the signal. The data may require more detailed modeling of the ice-basement transition. The fit to the surface receiver functions is not bad, but degrades with time this suggests lateral variation in the near-surface structure (ice or shallow subglacial basement). Figure S11. McMC search results for signals from station GM05. See supplement introduction for layout description. We consider these results quality C since our fit to the subsurface receiver functions is marginal. The data may require more detailed modeling of the ice-basement transition. The fit to the surface receiver functions is not bad the difference in frequency content suggests our shallow, sharp boundaries overestimate the abruptness of the near-surface material property changes. Figure S12. McMC search results for signals from station N140. See supplement introduction for layout description. We consider these results quality C. We do not fit the subsurface PRFs very well. The fit to the surface receiver functions is reasonable. Figure S13. McMC search results for signals from station N215. See supplement introduction for layout description. We consider these results quality C since our fit to the subsurface receiver functions fails to match the large amplitude signals. The fit to the surface receiver functions is not bad the difference in frequency content suggests our shallow, sharp boundaries overestimate the abruptness of the near-surface material property changes. Figure S14. McMC search results for signals from station P061. See supplement introduction for layout description. We consider these results quality B. The subsurface receiver functions are small, but the major features are matched and the McMC constraints reasonably tight. The fit to the surface receiver functions is good, but the signals come almost exclusively from the ice. Figure S15. McMC search results for signals from station SIPL. See supplement introduction for layout description. We consider these results quality C because the subsurface receiver functions are complex and the fits to the surface receiver functions are not very good. Still, the long-period features are matched and the large amplitude signals require a strong contrast between the glacial and subglacial domains. The deeper parts of the model are less well constrained. Figure S16. McMC search results for signals from station E028. See supplement introduction for layout description. We consider these results quality D because of the limited data available and large uncertainties in model parameters. The subsurface receiver functions are well matched, but the fit to the surface receiver functions is not very good. Table S1. Model S1, the model parameters are from Wittlinger & Farra (2015). Table S2. Model S2, the model parameters are from Wittlinger & Farra (2015). Table S3. Model S3. Table S4. Model S4. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. APPENDIX A: WAVEFIELD DOWNWARD CONTINUATION The formulas for wavefield downward continuation have been documented elsewhere such as Langston (2011) and Tao et al. (2014). We briefly review the equations for wavefield downward continuation for completeness. Assuming a plain-layer isotropic velocity model with m layers (layer 1, 2, 3, … m starting from the surface), wavevectors at the base of layer m can be represented as a function (eq. A1) of surface displacements, horizontal phase velocity and model parameters of the velocity model in frequency domain (Haskell 1953). \begin{equation} {\left[\begin{array}{cccc}\Delta ^{\prime }_m + \Delta ^{\prime \prime }_m \Delta ^{\prime }_m - \Delta ^{\prime \prime }_m \Omega ^{\prime }_m - \Omega ^{\prime \prime }_m \Omega ^{\prime }_m + \Omega ^{\prime \prime }_m \end{array}\right]} = \boldsymbol{J} {\left[\begin{array}{cccc}\dot{u}_0 / c \dot{w}_0 / c 0 0 \end{array}\right]} \end{equation} (A1) In eq. (A1), |$\Delta ^{\prime }_m$| is downgoing P wavevector, |$\Delta ^{\prime \prime }_m$| is upgoing P wavevector, |$\Omega ^{\prime }_m$| is downgoing S wavevector, |$\Omega ^{\prime \prime }_m$| is upgoing S wavevector, u0 is surface radial displacement, w0 is surface vertical displacement, c is the horizontal phase velocity and |$\boldsymbol{J}$| is a 4 by 4 matrix. The matrix |$\boldsymbol{J}$| can be computed using \begin{equation} \boldsymbol{J} = \boldsymbol{E^{-1}_m} \boldsymbol{a_{m-1}} \boldsymbol{a_{m-2}} \cdots \boldsymbol{a_1} \end{equation} (A2) with the explicit form of matrix |$\boldsymbol{E^{-1}_m}$| and |$\boldsymbol{a_m}$| given in Haskell (1953). Assuming a reference velocity model, all the terms on the right side of eq. (A1) are known. P and S wavevectors at a reference depth (base of layer m) can be calculated with eq. (A1). APPENDIX B: LIKELIHOOD FUNCTION AND RECEIVER FUNCTION ERROR STATISTICS The likelihood function |$L(\boldsymbol{m})$| is a function of misfit and the data covariance matrix. \begin{equation} L(\boldsymbol{m}) = A \exp \left(-\frac{1}{2} (G(\boldsymbol{m}) - \boldsymbol{d}_{{\rm obs}})^T \boldsymbol{C}_e^{-1} \left(G(\boldsymbol{m}) - \boldsymbol{d}_{{\rm obs}}\right)\right) \end{equation} (B1) in which |$\boldsymbol{C}_e$| is the data covariance matrix, A is a factor related to |$\boldsymbol{C}_e$|⁠, |$G(\boldsymbol{m})$| is the predicted waveform for model |$\boldsymbol{m}$|⁠, |$\boldsymbol{d}_{{\rm obs}}$| is the observed waveform, and T represents vector transpose. Since only the ratio of likelihood is used in the McMC search, we did not compute the factor A explicitly. Earlier studies (e.g. Sambridge 1999a; Agostinetti & Malinverno 2010; Bodin et al. 2012) have suggested that errors in surface PRFs are temporally correlated. Three types of data covariance matrices are commonly used including uniform-variance diagonal, nonuniform-variance diagonal and full matrix forms. Using the uniform variance covariance matrix assumes that data noise in PRFs is uniformly distributed with time and independent. A uniform-variance covariance matrix equals the identity matrix multiplied by a constant estimate of the data variance (average variance of the data). A non-uniform-diagonal covariance matrix accounts for the temporal variation of data noise but still assumes the noise is uncorrelated—diagonal elements equal the time-dependent data variance. The full covariance matrix takes the temporal variation of data variance and the temporal correlation of noise into consideration. The full covariance matrix has off-diagonal elements that are computed directly from observations. We used the Numpy function numpy.cov to estimate the covariance matrix for the suite of PRFs—a matrix is formed with columns consisting of the PRFs and the covariance matrix is estimated using the product of the matrix with its transpose. To visualize differences of these three types of covariance matrix, we compute synthetic error using all three types of covariance matrix derived from observations at Station BYRD. As shown in Fig. B1, synthetic error estimates computed using this covariance matrix exhibit similar variations as the data, much better than the diagonal approximations. Figure B1. Open in new tabDownload slide Comparison of synthetic data using different types of covariance matrices and real data. Thick black lines are stacked subsurface PRF observations at station BYRD for ray parameters between 0.05 and 0.06 s km−1. Thin black lines show one standard deviation error bounds. Red lines represent synthetic noise computed with uniform variance, nonuniform diagonal, and full covariance matrix on top of the stacked waveform for (a), (b) and (d), respectively. Red lines in (c) show individual subsurface PRFs prior to stacking. Figure B1. Open in new tabDownload slide Comparison of synthetic data using different types of covariance matrices and real data. Thick black lines are stacked subsurface PRF observations at station BYRD for ray parameters between 0.05 and 0.06 s km−1. Thin black lines show one standard deviation error bounds. Red lines represent synthetic noise computed with uniform variance, nonuniform diagonal, and full covariance matrix on top of the stacked waveform for (a), (b) and (d), respectively. Red lines in (c) show individual subsurface PRFs prior to stacking. The full covariance matrix can be used to account for the noise correlations in the likelihood calculations during the McMC search. To compute likelihood using the full covariance matrix, the inverse of the covariance matrix is required. Since the covariance matrix is usually singular, we compute the inverse of the covariance matrix using singular value decomposition (Sambridge 1999a). Eigenvalues that are more than 10−3 times (Sambridge 1999a) smaller than the highest eigenvalue are truncated. We compared McMC search results that used different types of covariance matrix (Fig. B2). The results indicate that using uniform variance or diagonal covariance matrix in the likelihood function underestimates the uncertainties of model parameters. For the example shown in Fig. B2, uncertainties obtained with uniform variance covariance matrix are only about 50 per cent of those using full covariance matrix. Figure B2. Open in new tabDownload slide Comparison of McMC results using different covariance matrices. Synthetic subsurface PRFs from the same model as in Fig. 8 were used. The covariance matrices were obtained from observations at station BYRD. (a), (d) and (g) correspond to a uniform variance assumption; (b), (e) and (h) correspond to a nonuniform, diagonal covariance matrix assumption; and (c), (f) and (i) correspond to a full covariance estimate. (a)–(c) show accepted velocity models from McMC searches. The error ellipses represent 95 per cent confidence interval (outer) and 68 per cent confidence interval (inner). The colour scales for (a), (b) and (c) are the same as that in Fig. 8. Marginal distributions of crustal thickness are shown in (d)–(f). Marginal distributions of average crustal shear wave velocity are shown in (g)–(i). Dashed lines indicate true values. Figure B2. Open in new tabDownload slide Comparison of McMC results using different covariance matrices. Synthetic subsurface PRFs from the same model as in Fig. 8 were used. The covariance matrices were obtained from observations at station BYRD. (a), (d) and (g) correspond to a uniform variance assumption; (b), (e) and (h) correspond to a nonuniform, diagonal covariance matrix assumption; and (c), (f) and (i) correspond to a full covariance estimate. (a)–(c) show accepted velocity models from McMC searches. The error ellipses represent 95 per cent confidence interval (outer) and 68 per cent confidence interval (inner). The colour scales for (a), (b) and (c) are the same as that in Fig. 8. Marginal distributions of crustal thickness are shown in (d)–(f). Marginal distributions of average crustal shear wave velocity are shown in (g)–(i). Dashed lines indicate true values. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Estimating subglacial structure using P-wave receiver functions
JF - Geophysical Journal International
DO - 10.1093/gji/ggx075
DA - 2017-05-01
UR - https://www.deepdyve.com/lp/oxford-university-press/estimating-subglacial-structure-using-p-wave-receiver-functions-iiBh4ywVbi
SP - 1064
VL - 209
IS - 2
DP - DeepDyve
ER -