# Wavefield iterative deconvolution to remove multiples and produce phase specific Ps receiver functions

Wavefield iterative deconvolution to remove multiples and produce phase specific Ps receiver... Summary Common conversion point stacking or migration of receiver functions (RFs) and H–k (H is depth and k is Vp/Vs) stacking of RFs has become a common method to study the crust and upper mantle beneath broad-band three-component seismic stations. However, it can be difficult to interpret Pds RFs due to interference between the Pds, PPds and PSds phases, especially in the mantle portion of the lithosphere. We propose a phase separation method to isolate the prominent phases of the RFs and produce separate Pds, PPds and PSds ‘phase specific’ receiver functions (referred to as PdsRFs, PPdsRFs and PSdsRFs, respectively) by deconvolution of the wavefield rather than single seismograms. One of the most important products of this deconvolution method is to produce Ps receiver functions (PdsRFs) that are free of crustal multiples. This is accomplished by using H–k analysis to identify specific phases in the wavefield from all seismograms recorded at a station which enables development of an iterative deconvolution procedure to produce the above-mentioned phase specific RFs. We refer to this method as wavefield iterative deconvolution (WID). The WID method differentiates and isolates different RF phases by exploiting their differences in moveout curves across the entire wave front. We tested the WID by applying it to synthetic seismograms produced using a modified version of the PREM velocity model. The WID effectively separates phases from each stacked RF in synthetic data. We also applied this technique to produce RFs from seismograms recorded at ARU (a broad-band station in Arti, Russia). The phase specific RFs produced using WID are easier to interpret than traditional RFs. The PdsRFs computed using WID are the most improved, owing to the distinct shape of its moveout curves as compared to the moveout curves for the PPds and PSds phases. The importance of this WID method is most significant in reducing interference between phases for depths of less than 300 km. Phases from deeper layers (i.e. P660s as compared to PP220s) are less likely to be misinterpreted because the large amount of moveout causes the appropriate phases to stack coherently if there is sufficient distribution in ray parameter. WID is most effective in producing clean PdsRFs that are relatively free of reverberations whereas PPdsRFs and PSdsRFs retain contamination from reverberations. Composition and structure of the continental crust, Composition and structure of the mantle, Body waves, Crustal imaging 1 INTRODUCTION Receiver function (RF) techniques are widely used to make 2-D and 3-D stacked images of crust and upper mantle structure beneath networks of three component broad-band seismic stations (Dueker & Sheehan 1998; Simmons & Gurrola 2000; Gilbert et al. 2003; Liu et al. 2003; Gao et al. 2002; Andrews & Deuss 2008). Single station analysis of crust and upper mantle structure is commonly performed by inversion of RF waveforms (Ammon et al. 1990; Julia et al. 2000; Bodin et al. 2012), H–k analysis to estimate Moho depths and Vp/Vs (Zhu & Kanamori 2000; Dugda et al. 2005; Nair et al. 2006) and studies of the upper mantle transition zone (Gurrola et al. 1994; Gurrola & Minster 1998; Gilbert et al. 2003; Lawrance & Shearer 2006; Tauzin et al. 2008). A problem with all of these methods is that there can be interference (Fig. 1) that can complicate the interpretation of these data. Figure 1. View largeDownload slide Illustration of ray paths for P, Pds, PPds and PSds phases. Figure 1. View largeDownload slide Illustration of ray paths for P, Pds, PPds and PSds phases. In RF analysis, multiple traces can be grouped and stacked after correcting for moveout (time delays at different ray parameters). The arrival times of the different Ps phases are delayed with increasing ray parameter, whereas the moveout for reverberations (i.e. PPds and PSds phases) is advanced with increasing ray parameter. For 1-D structure, these moveout times are easily computed (Vinnik 1977). By stacking moveout corrected RFs from a range of ray parameters using a range of reasonable velocity models, one can estimate the depth and average velocity structure to the boundary via velocity spectrum stacking (VSS, Gurrola et al. 1994; Gurrola & Minster 1998) or estimate depth and Vp/Vs ratio (H–k method, Zhu & Kanamori 2000) to these boundaries. However, in all these traditional Ps RF methods, reverberations from shallow depths (e.g. mid-crust and Moho) can be mis-stacked as Ps phases or interfere with Ps phases from greater depths. To avoid misinterpretation of phases and to better stack RFs, we developed a method to separate Pds, PPds and PSds phases from RFs and produce separate, single phase type receiver functions (PdsRFs, PPdsRFs and PSdsRFs) from each seismogram. The wavefield iterative deconvolution (WID) method is introduced to differentiate and isolate different RF phases by exploiting differences in moveout curves for each of the common phase types. This is a phase separation method to produce separate, ‘phase specific’ Pds, PPds and PSds receiver functions. We proceed by first illustrating the potential effectiveness of the WID method using synthetic RFs computed using raytracing (Aki & Richards 2002). Then, we apply the WID technique to data from the broad-band seismic station at Arti, Russia (ARU). We chose ARU as a test station because this station has high quality data from a wide range of ray parameters and the crust and mantle structure there is relatively uncomplicated as compared to other stations around the world. 2 METHOD Ps RF analysis is a method to isolate Pds converted phases and reverberations of the P-phase that end in an S-phase that originate from the structure beneath a three-component seismic station. Since Pds conversions are SV phases, we separate the SV and SH phases by rotating these data from the N--S and E–W components to compute radial and transverse components (Plesinger et al. 1986). Deconvolving the radial and transverse components by the vertical component produces RFs that isolate the Pds phases and reverberations ending in an S phase (Lang 1979; Ammon 1991). Frequency domain deconvolution by spectral division of the horizontal components by the vertical component after pre-whitening is one of the oldest methods of deconvolution used to produce RFs (Langston 1977; Owens et al. 1988). The time-domain iterative deconvolution method, introduced to RF analysis by Ligorria & Ammon (1999), lends itself well to analysis that incorporates the shape of the moveout curves to distinguish Pds, PPds and PSds phases from one another. We accomplish this by including all the seismograms recorded at a seismic station into a simultaneous deconvolution that uses a modified H–k analysis to identify the wavefield of the strongest remaining phase in each iteration and construct three, phase specific RFs from each seismogram in the data set. We adopt the regularized iterative deconvolution method (Wang & Pavlis 2016) that employs a regularized deconvolution of the seismogram, rather than cross correlation, to identify phases in each iteration of the deconvolution. Wang & Pavlis (2016) give a detailed discussion of the benefits of the regularized iterative deconvolution over the use of cross-correlation. In a head-to-head comparison of five common methods of computing RFs, Huang et al. (2016) also found the regularized iterative deconvolution to be the most robust of the methods tested. Fig. 1 illustrates the ray paths of converted Pds phases and major reverberations ending in converted S-phases (PPds and PSds). The relative arrival times of any phase in a RF for a horizontal interface is the arrival time of the phase of interest minus the arrival time of the direct P phase (Langston 1977; Vinnik 1977; Gurrola et al. 1994) and, for 1-D velocity models, can be represented by the following equations:   $$\Delta {T_{Ps}} = Z*\left( {\sqrt {{V_s}^{ - 2} - {p^2}} - \sqrt {{V_p}^{ - 2} - {p^2}} } \right)$$ (1)  $$\Delta {T_{PPds}} = Z*\left( {\sqrt {{V_s}^{ - 2} - {p^2}} + \sqrt {{V_p}^{ - 2} - {p^2}} } \right)$$ (2)  $$\Delta {T_{PSds}} = 2*Z\left( {\sqrt {{V_s}^{ - 2} - {p^2}} } \right),$$ (3)where p is the ray parameter of the incident wave, Vs is the S-wave velocity, and Z is the depth to the interface. In RF work, Vp is typically replaced by k*Vs (where k = Vp/Vs) in eqs (1) and (2) so that we can compute traveltimes from a Vs model and an assumed (or estimated) value for k. By writing PPds and PSds traveltime equations (eqs 2 and 3) as a function of k and Vs rather than as a function of Vp and Vs, we can simply use one equation for both types of reverberations and use k as the discriminating criteria between PPds and PSds phases:   $$\Delta {T_{PXds}} = Z*\left( {\sqrt {{V_s}^{ - 2} - {p^2}} + \sqrt {{{(k*{V_s})}^{ - 2}} - {p^2}} } \right).$$ (4)In this case, k = 1 will correspond to PSds phases and eq. (4) would reduce to be the same as eq. (3). The WID algorithm described below identifies Pds, PPds and PSds phases from single-phase H–k analysis to produce phase specific RFs. To perform WID we: Calculate two ‘single phase’ H–k matrices, one for Pds phases and the other for PXds phases using traveltime eqs (1) and (4), respectively. Any reasonable 1-D S-wave velocity model can be used as a reference model in H–k stacking. To get optimal results from the H–k analysis, it is necessary to have data from several ray parameter bins. The algorithm searches for the largest amplitude across the two H–k matrices and identifies the phase by:   \begin{equation*} Phase\ type = \left( {\begin{array}{@{}*{1}{c}@{}} {\rm Pds \ if \ largest \ phase \ is \ on \ the \ Pds} \ H - k \ {\rm matrix}\\ {\rm PPds \ if \ largest \ phase \ is \ on \ the \ PXds} \ H - k \ {\rm matrix \ and} \ 1.5 < k < 1.95\\ {\rm PSds \ if \ largest \ phase \ is \ on \ the \ PXds} \ H - k \ {\rm matrix \ and} \ 0.9 < k < 1.1 \end{array}} \right)\end{equation*} We then use cross-correlation over a narrow time window centred around the expected arrival time of this phase to estimate the ‘real’ arrival time on each seismogram. To assure that the arrival times and amplitudes for the identified phase on each of the phase specific RFs reflects the original data, we only use steps one and two to estimate the depth and k value for each phase to identify its phase type. We then construct 3 phase specific RFs (spike functions) for each seismogram by placing the amplitude of the identified phase at the time delays (determined in step 3) on the phase specific RFs. The phase specific RFs are convolved by the vertical component to produce the three different, phase dependent seismograms. All three of the phase dependent seismograms are subtracted from the observed seismograms to produce modified ‘observed’ seismograms with all previously identified phases removed (Ligorria & Ammon 1999). Steps 1 through 6 are repeated for many iterations until all significant phases are removed from the original seismograms. The number of iterations needed to complete this process is determined by calculating the sum of the standard deviations of the modified seismograms during each iteration. The algorithm stops when sum of the standard deviations between iterations becomes small or nearly the same. 2.1 Test of WID with synthetic data To test the WID method described above, we applied WID to synthetic seismograms computed for the teleseismic range of ray parameters (0.04–0.08 s km–1, Fig. 2). These synthetics were computed by ray tracing through the modified PREM velocity model (Dziewonski & Anderson 1981). We modified PREM by adding a 2 km thick sedimentary surface layer to simulate complications that may occur at many seismic stations. We also added a 75 km deep low velocity layer to simulate arrivals that are of interest in recent RF investigations of the lithosphere. The amplitude of each phase was determined by reflection and transmission coefficient equations from Aki & Richards (2002) and a 0.5 Hz Gaussian low-pass filter was applied. Figure 2. View largeDownload slide Synthetic receiver function as a function of ray parameter computed from a modified version of the PREM velocity model. We labelled the prominent Pds phases from the Moho, the 75 km low velocity layer (that we added to the velocity model), the 220, 410 and 660 km discontinuities. The phases from the surface layer and Mid-crustal discontinuities were not labelled so as not to clutter the figure. Reverberations from Moho and higher order reverberations are also labelled. Figure 2. View largeDownload slide Synthetic receiver function as a function of ray parameter computed from a modified version of the PREM velocity model. We labelled the prominent Pds phases from the Moho, the 75 km low velocity layer (that we added to the velocity model), the 220, 410 and 660 km discontinuities. The phases from the surface layer and Mid-crustal discontinuities were not labelled so as not to clutter the figure. Reverberations from Moho and higher order reverberations are also labelled. The Pds phases from the surface layer, mid-crust, 75 km low velocity layer (LVL), Moho, 220, 410 and 660 km discontinuities are clearly visible in these synthetics (Fig. 2). Pds traveltimes increase as ray parameter increases, which is opposite to the moveout of the PXds phases that arrive earlier as ray parameter increases. The curvature of Pds traveltimes (moveout) becomes larger as conversion depth increases. For example, it is harder to see the moveout of the mid-crust and Moho Pds phases than for the Pds phase from deeper layers such as P220s, P410s and P660s. The difference in the shape of the moveout curve between Pds and PXds phases can be used to distinguish them from one another. For example, notice that the P660s and PP220s phases curve in opposite directions. We produced ‘single phase’ H–k plots (Pds and PXds H–k plots) using the synthetic data shown in Fig. 2. Fig. 3(a) is the single phase H–k plot computed by stacking these RFs using the Pds traveltime equation (eq.1). We refer to this as a single phase H–k plot because it is a plot of stacked Pds phase amplitudes as a function of H and k values and should not be confused with the H–k method by Zhu & Kanamori (2000) that calculates the summation of the weighted stack of Pds, PPds and PSds phases. The largest stacked amplitudes are for Pds phases from the near surface layer, Moho, 220, 410 and 660 km discontinuities and are centred near k = 1.81 (marked by stars). However, the peaks corresponding to the Pds phases from mid-crust and 75 km LVL appear at a k value of ∼1.77. Figure 3. View largeDownload slide H–k plot using Pds moveout curve (a). Pds phases from surface layer, Mid-crust, Moho, 75 km low velocity layer, 220, 410 and 660 discontinuities are labelled; H–k plot using PXds moveout curve (b). Mid-crust, Moho and 220 PPds, Mid-crust and Moho PSds are labelled. Stars represent location of highest amplitude for each phase. Figure 3. View largeDownload slide H–k plot using Pds moveout curve (a). Pds phases from surface layer, Mid-crust, Moho, 75 km low velocity layer, 220, 410 and 660 discontinuities are labelled; H–k plot using PXds moveout curve (b). Mid-crust, Moho and 220 PPds, Mid-crust and Moho PSds are labelled. Stars represent location of highest amplitude for each phase. We see in Fig. 3 that, even for noise free synthetics, the peak for the largest amplitude is elongated and not sharp. There can be large errors in picking this peak, but in WID the estimated value for k only needs to be accurate enough to distinguish between the three phase types. If the phase is a PPds or PSds phase, it will have its largest peak along the left edge of the Pds H–k plot (Fig. 3a); otherwise it will be near the correct value. Likewise, a Pds phase will plot on the left edge of the PXds H–k plot. Thus, it is unlikely that errors in the estimated k value will cause a Pds phase to be incorrectly identified as a PXds phase. Highly accurate values for k are also not essential to distinguish between the PSds and PPds phases on the PXs H–k plot (Fig. 3b), since the PSds phase is identified by a k value near 1, whereas the PPds is identified by having a value greater than 1.5. As a result, an error analysis of the k value is not essential because our method is effective even in the presence of relatively large errors in the estimated value for k. Fig. 3(b) is the H–k plot computed by stacking RFs using the PXds traveltime equation (eq. 4). To better understand how each phase is identified and placed on the appropriate phase specific RF, we produced plots of the amplitudes, H and k values for each pulse added to each of the phase specific RFs (Fig. 4). Each circle in Fig. 4(a) represents the H and k values that correspond to the largest stacked amplitude of a Pds phase from each iteration of the WID. Fig. 4 shows that many of the peaks on these amplitude versus H–k plots do not match the value of 1.8 used to produce the synthetics in Fig. 2. These errors in k values are due to interference between different phases (i.e. interference between the PPds phase and higher order reverberations from the surface layer). A k value of 1.81 was found for the Moho and 220 PPds phases, but the largest amplitude of the mid crust PPds phase occurred at k = 1.6. The peak amplitude of the PSds phase for the mid crust and Moho is between k = 0.9 and k = 1 in Figs 4(b) and (c). Errors in the estimated k values for the PPds and PSds phases from the surface layer (k 1.5 and k 1.06, respectively) do not prevent them from being identified as the appropriate phase. In general, we expect the errors for k to be larger for shallow layers since they have the flattest moveout curves. Figure 4. View largeDownload slide H & k value and amplitudes for phases identified in each iteration of the phase separation portion of the WID. The relative amplitude for the phase is illustrated by the size of the circle. Separate plots (a, b and c) correspond to the Pds, PPds and PSds, respectively. Figure 4. View largeDownload slide H & k value and amplitudes for phases identified in each iteration of the phase separation portion of the WID. The relative amplitude for the phase is illustrated by the size of the circle. Separate plots (a, b and c) correspond to the Pds, PPds and PSds, respectively. We applied the WID method discussed above to the synthetics to produce three separate phase dependent RFs (PdsRFs, PPdsRFs and PSdsRFs are shown in Figs 5a–c respectively). Fig. 5(d) is the residual of the original seismograms after subtracting all the phase specific RFs from the original synthetics. The Pds phases from the surface layer, mid crust, Moho, 75 km LVL, 220, 410 and 660 km are clearly visible on the PdsRFs (Fig. 5a). The PPds phases from the mid crust, Moho and 220 are the strongest phases on PPdsRFs but many higher order reverberations of these phases from shallower layers are present as well (Fig. 5b). The PPds phase from 75 km LVL is present and strong but is hidden among the higher order reverberations from shallower layers. The PPdsRFs are contaminated with more high order reverberations than PSdsRFs. The PdsRFs are the least contaminated by other phases. The residual RFs (Fig. 4d) show that there are some artefacts in the data after removing the strong RF phases, but these artefacts would likely be below the noise level of real data. Figure 5. View largeDownload slide Receiver functions produced from the synthetics in Fig. 2 after phase separation: (a) Pds receiver function (b) PPds receiver function (c) PSds receiver function and (d) Residual receiver function Figure 5. View largeDownload slide Receiver functions produced from the synthetics in Fig. 2 after phase separation: (a) Pds receiver function (b) PPds receiver function (c) PSds receiver function and (d) Residual receiver function The utility of WID is best illustrated by a comparison of the stacks of the phase specific RFs after applying moveout corrections. The stack of the PdsRFs (Fig. 6 trace b) only has Pds phases from the surface layer, mid-crust, Moho, 75 km LVL, 410 and 660 km discontinuities, whereas the stacked RF produced from traditional RFs (Fig. 6 trace a) is contaminated by mis-stacked PPds phases from the mid-crust and Moho that that appear at 90 and 140 km depths. Likewise, PSds phases from the shallow layer, mid-crust and Moho appear as mis-stacked negative phases at about 5, 110 and 190 km depths. The higher order reverberations are mis-stacked as phases between the 250 and 350 km depths. The PPds phase from the 220 has sufficient moveout so as not to ‘mis-stack’ on any of these stacked RFs. So, WID is most effective for differentiating phases from shallower depths. Figure 6. View largeDownload slide Pds stacks of the synthetic RF produced from the noise free synthetics without WID phase separation (a) and with WID phase separation (b). Also shown are stacked Pds RFs computed from noisy synthetic without WID phase separation (c) and with WID phase separation (d). Figure 6. View largeDownload slide Pds stacks of the synthetic RF produced from the noise free synthetics without WID phase separation (a) and with WID phase separation (b). Also shown are stacked Pds RFs computed from noisy synthetic without WID phase separation (c) and with WID phase separation (d). To test the robustness of WID in the deconvolution of noisy data, we added 10 per cent noise to the synthetic data and repeated the above tests of the WID method (Fig. 6). The noise added to the synthetics is the residual after WID of the data from the seismic station ARU (see the discussion below) after scaling the residual noise to have a standard deviation that is 10 per cent of the peak signal in the synthetics. The stacked PdsRFs computed with added noise (Fig. 6d) and without added noise (Fig. 6d) are nearly identical. Fig. 7 shows the stacked RFs after applying PPds moveout corrections and converting to depth. Fig. 7(a) is the stack of traditional RFs, and traces b, c and d are stacks of the PPdsRFs produced using WID with different criteria. Once again, the stacked traditional RFs (Fig. 7a) are contaminated by mis-stacked phases, whereas the phase separated PPdsRFs (Figs 7b–d) have much less contamination from other mis-stacked phases. The shallowest phases on the traditional PPds RF (Fig. 7a) are mis-stacked Pds phases. These shallow phases have very little moveout so they can appear strongly on the wrong stacked RFs. There are several phases in the 50–200 km depth ranges that have very small amplitudes. These are either higher order reverberations of the PXds phases or mis-stacked Pds phases from greater depths. By changing the phase identification criteria in the WID, we may be able to clean up the PPdsRFs. By narrowing the range of allowable k values to 1.75 through 1.85, the stacked PPdsRF (Fig. 7c) contains fewer higher order reverberations, but it is also missing some of the expected PPds phases (i.e. the PPds phases from the mid-crustal and the 75 km LVL). The stack of PPdsRFs with only the 10 largest PPds phases (Fig. 7d) identified by WID includes more of the expected phases but also retains some of the large amplitude higher order reverberations. From this analysis, it is obvious that WID is not as effective in producing ‘clean’ PPdsRFs as it is in producing a relatively uncontaminated stacked PdsRFs. An argument can be made that WID is only effective in producing clean PdsRFs, but this analysis shows that WID can be used to extract useful information from the PPds phase and produce a better stacked PPdsRFs than traditional processing. The PPdsRFs produced from the same noisy synthetics used in the analysis of the PdsRFs above have a few small phases not present in the respective noise free phase specific stacked PPdsRFs, but these are very small peaks and may be the result of reverberations that were not completely removed from the ARU residuals. Figure 7. View largeDownload slide Stacked PPds receiver functions produced from the synthetics in Fig. 2. Traces a and e are the stacked RF produced from traditional RFs without noise and with real seismic noise added, respectively. Traces b, c, d and f are stacked, phase specific PPdsRFs produced using the WID method and noise free PPdsRFs including phases with k values between 1.5 and 1.9 (b), noise free PPdsRFs including phases with k values between 1.75 and 1.85 (c), noise free PPdsRFs with only the phases with the 10 largest amplitudes (d), and PPdsRFs produced from synthetics with real seismic noise added and including phases with k values between 1.5 and 1.9 (f). Figure 7. View largeDownload slide Stacked PPds receiver functions produced from the synthetics in Fig. 2. Traces a and e are the stacked RF produced from traditional RFs without noise and with real seismic noise added, respectively. Traces b, c, d and f are stacked, phase specific PPdsRFs produced using the WID method and noise free PPdsRFs including phases with k values between 1.5 and 1.9 (b), noise free PPdsRFs including phases with k values between 1.75 and 1.85 (c), noise free PPdsRFs with only the phases with the 10 largest amplitudes (d), and PPdsRFs produced from synthetics with real seismic noise added and including phases with k values between 1.5 and 1.9 (f). Fig. 8 shows the stack of PSdsRFs from Fig. 5(c) and the stack of traditional RFs. Once again, the PSdsRFs (Fig. 8b) only show legitimate, strong PSds phases from the mid-crust, Moho and 75 km LVL. The traditional stacked PSdsRFs (Fig. 8a) have strong phases at shallow depths that are mis-stacked Pds and PPds phases and include the poorly stacked Pds phases from greater depths in the 130–180 km depth interval. The PSdsRFs (Fig. 8d) produced using the noisy synthetics are relatively unaffected by the noise. Figure 8. View largeDownload slide Stacked PSds receiver functions produced from the synthetics in Fig. 2. Traces a and c are the stacked RF produced from traditional RFs without noise and with real seismic noise added, respectively. Traces b and d are stacked, phase specific PSdsRFs produced using the WID method with noise free PSdsRFs (b) and PSdsRFs produced from synthetics with real seismic noise added. Figure 8. View largeDownload slide Stacked PSds receiver functions produced from the synthetics in Fig. 2. Traces a and c are the stacked RF produced from traditional RFs without noise and with real seismic noise added, respectively. Traces b and d are stacked, phase specific PSdsRFs produced using the WID method with noise free PSdsRFs (b) and PSdsRFs produced from synthetics with real seismic noise added. WID is effective in selecting the expected phase type for each of the phase specific stacked RFs computed from synthetic data. In the case of the PPdsRFs and PSdsRFs, higher order reverberations of the expected phases are still present, but the PdsRFs are free of such reverberations. The importance of this WID method is most significant in cleaning up the RFs at depths of less than 300 km. The deeper converted phases are less likely to be misinterpreted in stacked RFs because there is less cross stacking between Pds and PXds phases since the amount of moveout is much greater than the pulse width of these deeper phases. However, PPds and PSds phases can cross contaminate each other to greater depths because the moveout curves are similar in shape. 2.2 Test of WID using real data We used data from the GDSN station ARU (Arti, Russia) to test the utility of the WID method applied to real data. We analysed all the data recorded from 1990 through 2016 at ARU and found 649 events that produced RFs with high signal-to-noise ratios (Fig. 9). The epicentral distance for this data set ranges from 30° to 90° (corresponding to a range of ray parameters of 0.04 to 0.08 s km–1). The backazimuths with data coverage across a wide range of ray parameters are mainly from 80° to 140°. Data were windowed from 30 s prior to and 150 s after the direct-P wave arrival. A data stream of this length is sufficient to include all primary reverberations from lithospheric depths and Pds phases through the upper mantle transition zone. We rotated these seismograms into the L-Q-T coordinate system after estimating the incident angle of incoming P wave using singular value decomposition to find the angle that maximizes the amplitude of P wave energy on the Q component. By doing so, we maximize the isolation of SV and P-wave energy on the L and Q component, respectively (Plesinger et al. 1986). An advantage of using the L-Q-T system is that it makes it possible to detect phases from shallow depths by avoiding side lobes from the direct P. By minimizing the contamination of the P phase (Q component) with SV energy, we also have a cleaner source function for deconvolution. Figure 9. View largeDownload slide Station and earthquakes epicentres recorded from 1990 to 2016 used in our study. Black triangle shows ARU station location; black circles show the 649 earthquakes used in this study. Figure 9. View largeDownload slide Station and earthquakes epicentres recorded from 1990 to 2016 used in our study. Black triangle shows ARU station location; black circles show the 649 earthquakes used in this study. We used simultaneous frequency domain water level deconvolution to produce the unfiltered SV RFs by deconvolving the Q component from the L component of the data. Prior to deconvolution, we applied a bandpass filter to the data using corner frequencies of 0.02 and 0.5 Hz. We clustered the events into ray parameter bins 0.001 s km–1 wide to enable simultaneous deconvolution of all the events in a bin (Gurrola et al. 1995). Fig. 10 shows the traditional RFs computed for the ARU station. These RFs were shifted so the first arrival would be at 10 s. The dominant phases are Pds phases from the surface layer and the Moho, which appear at 11 and 14 s. Other phases are present but they are masked by the noise. In the next section, we apply the WID method to these RFs to produce PdsRFs, PPdsRFs and PSdsRFs. Figure 10. View largeDownload slide Receiver functions computed for ARU station with 0.02–0.5 Hz bandpass filter applied. Surface layer and Moho Pds phases are the only clearly observable phases. Figure 10. View largeDownload slide Receiver functions computed for ARU station with 0.02–0.5 Hz bandpass filter applied. Surface layer and Moho Pds phases are the only clearly observable phases. The number of iterations included in the WID was determined by computing the sum of the standard deviation of the residual seismograms after each iteration of the WID (Fig. 11), and the WID was stopped when the sum of the standard deviation changed by very little between iterations. In our case, 120 iterations appeared to be a sufficient number of iterations. Figure 11. View largeDownload slide Sum of the standard deviations of all receiver function traces versus number of iteration through iterative deconvolution. Figure 11. View largeDownload slide Sum of the standard deviations of all receiver function traces versus number of iteration through iterative deconvolution. Fig. 12(a) is a ‘single phase’ H–k plot computed by stacking RFs computed from ARU data using Pds traveltime equations for a range of k values using PREM S-velocities as the reference model. The largest stacking amplitude occurs for the Moho at ∼38 km in depth with k = 1.87. The mid lithosphere at 95 km depth, 410 and 660 Pds phases are also identified at k = 1.87, 1.88 and 1.89, respectively. The k values found for ARU data are different than those found for the PREM synthetics. PREM velocities are not an exact match to the velocity model appropriate for ARU so we would not expect accurate k values. However, k is only used here as a discriminant to identify phase types; any reasonable 1-D velocity model works for WID. Figure 12. View largeDownload slide Single phase H--k plot using Pds moveout curve for station ARU receiver function (a) and the H–k plot using PXds moveout curve (b). Pds phases from Moho, mid-lithosphere, 410 and 660 discontinuities are labelled; Moho PPds and PSds phases are also labelled. Figure 12. View largeDownload slide Single phase H--k plot using Pds moveout curve for station ARU receiver function (a) and the H–k plot using PXds moveout curve (b). Pds phases from Moho, mid-lithosphere, 410 and 660 discontinuities are labelled; Moho PPds and PSds phases are also labelled. To verify that the WID is not dependent on variations in a reasonable velocity model, we preformed WID on the synthetics using a modified version of PREM (changing PREM by ±10 per cent) and found that stacked phase specific RFs were virtually identical regardless of which velocity model was used. Fig. 12(b) is an H–k plot computed by stacking RFs using the PXds traveltime equation (eq. 4). We find the maximum amplitude for the Moho PPds phase is at k = 1.5, and the peak amplitude for Moho PSds occurs at k = 0.95 (Figs 12b, 14b and c). These H–k plots appear to have more phases than expected for the number of layers in our velocity model. These extra phases are stacked higher order reverberations whose amplitude are smaller than expected phases. However, the expected PPds and PSds phases are clearly observable on the phase specific RFs. Fig. 13 shows the PdsRFs (a), PPdsRFs (b) and PSdsRFs (c) produced by WID for the ARU data. We are able to identify Pds phases from a shallow surface layer, the Moho, a mid-lithosphere discontinuity (at 95 km depth), the 410 discontinuity and the 660 km discontinuity in the PdsRFs (Fig. 13a). All these phases are clearly delayed as ray parameter increases, which is expected for the Pds phase. The unlabelled phases at around 20 s are higher order Pds reverberations from shallower layers. The amplitude versus offset (AVO) relations show the expected pattern of increasing amplitude with larger ray parameter. Figure 13. View largeDownload slide Individual phase receiver function produced from ARU data after phase separation (a), (b) and (c) shows the Pds, PPds and PSds receiver function, respectively. Figure 13. View largeDownload slide Individual phase receiver function produced from ARU data after phase separation (a), (b) and (c) shows the Pds, PPds and PSds receiver function, respectively. Fig. 14 shows depths, amplitudes and k values of all the phases in the RFs produced using WID on the ARU data. Fig. 14(a) shows the H and k values used to construct the PdsRFs in Fig. 13(a). All the dots for the Psd phase in Fig. 14(a) are at depths consistent with expected discontinuities, but the k values vary, reflecting errors due to the widths of the peaks on the H–k plots, errors in the assumed velocity model, contamination by noise and interference from other phases such as higher order reverberations. The multiple k values for the same depth (H), such as the Moho depth and the surface layer, reflect that, in WID, the same depth can be revisited between iterations. It is noteworthy that largest amplitudes found for the Moho Pds (Fig. 14a), PPds (Fig. 14b) and PSds (Fig. 14c) phases are at a depth of ∼38 km. The k value of 1.88 found for the P410s and P660s phases is higher than computed from PREM (which is 1.85 for 410 and 1.84 for 660), but in WID we only use the k for phase identification and not to interpret the local geology. Errors in the assumed velocity model will leak into the errors in the estimated k. However, our goal here is to clean up RFs to produce individual phase specific PdsRFs, PPdsRFs and PSdsRFs in the time domain, not to estimate k. The PPdsRFs (Fig. 13b) appear to have strong phases from the Moho (at 15 s) and mid lithosphere (at 35 s), but, as with the synthetics, the PPdsRFs are more heavily contaminated by higher order reverberations than the PdsRFs or the PSdsRFs. The range of k values for the PPds phases (Fig. 14b) is larger than found for the Pds phases, because the PPds plot includes many higher order reverberations, which should not stack at the correct k value. We can suppress reverberations by using a narrower range of k values, but we may not want to do so, because these reverberations may be useful in waveform modelling. Figure 14. View largeDownload slide H & k value for each iteration during phase separation (a), Pds receiver function amplitudes, (b) PPds receiver function amplitudes and (c) PSds receiver function amplitudes. Figure 14. View largeDownload slide H & k value for each iteration during phase separation (a), Pds receiver function amplitudes, (b) PPds receiver function amplitudes and (c) PSds receiver function amplitudes. The Moho and the mid lithosphere phases in the PSdsRFs (Fig. 13c) are easily identified since they have reversed amplitudes compared to PdsRFs and PPdsRFs. The k values for PSds are centred around a k 1 (Fig. 14c), which is consistent with what we expected for PSds phases in setting up the PXds phase moveout eq. (4). Fig. 15 shows the Pds, PPds and PSds produced by WID and traditional RF methods. The Moho and mid lithosphere Pds phases are separated with little interference from other phases in the PdsRFs (Fig. 15a). There appears to be two phases from the Moho depth on the traditional Pds RF (Fig. 15b), but one of the two is likely a mis-stacked PPds phase from a shallower layer. In the PdsRFs (Fig. 15a), the shallower of the two phases were removed. Moho PPds and PSds are mis-stacked at depths of about 150 and 180 km, respectively, in the traditional Pds RFs (Fig. 15b), but these phases were removed from the PdsRFs (Fig. 15b) computed using WID. Figure 15. View largeDownload slide Stacks of receiver functions using all 649 events recorded at ARU for depths of up to 200 km: (a) stack of PdsRFs produced using WID, (b) stack of Pds RFs using traditional RFs, (c) stack of WID PPdsRF only including the phases with the eight largest amplitudes, (d) stack of WID PPdsRFs only including those phases identified with Vp/Vs ratios of between 1.85 to 1.89, (e) stack of WID PPdsRFs that includes phases identified with Vp/Vs ratios between 1.4 and 1.9, (f) stack of traditional PPds RFs, (g) stack of PSdsRFs produced using WID, and (h) stack of traditional PSds RFs. Note that Moho and mid-lithosphere phases are labelled since they are shown in all three stacks. Moho phases stack coherently at a depth of ∼38 km and a mid-lithosphere discontinuity is at 95 km depth. Figure 15. View largeDownload slide Stacks of receiver functions using all 649 events recorded at ARU for depths of up to 200 km: (a) stack of PdsRFs produced using WID, (b) stack of Pds RFs using traditional RFs, (c) stack of WID PPdsRF only including the phases with the eight largest amplitudes, (d) stack of WID PPdsRFs only including those phases identified with Vp/Vs ratios of between 1.85 to 1.89, (e) stack of WID PPdsRFs that includes phases identified with Vp/Vs ratios between 1.4 and 1.9, (f) stack of traditional PPds RFs, (g) stack of PSdsRFs produced using WID, and (h) stack of traditional PSds RFs. Note that Moho and mid-lithosphere phases are labelled since they are shown in all three stacks. Moho phases stack coherently at a depth of ∼38 km and a mid-lithosphere discontinuity is at 95 km depth. As discussed for the synthetic tests above, the phase separated PPdsRF stacks have more peaks than PdsRFs (Fig. 16a) and PSdsRFs (Fig. 16g) when phases from a wide range of k are included (Fig. 15c). When we further restrict the range of allowable k values (Fig. 15d) and amplitudes (Fig. 15c) in phase selection during WID, we can produce cleaner PPdsRFs. Fig. 15(c) shows the PPdsRFs with only the eight largest amplitude arrivals included, resulting in PPdsRFs with phases from the Moho, but without the mid lithosphere phase which has a small amplitude. By using a tighter range of k values (Fig. 15d), we are able retain the phases from the mid-lithosphere and an apparent double phase from the Moho depth. However, this PPdsRF has less contamination from higher order reverberations than are visible in the PPdsRF computed using a wider range of k values. It is also possible that some low amplitude phases are side lobes or artefacts leaked in during the WID. Nevertheless, the Moho PPds and PSds phases are stacked at a depth of about 38 km, which is consistent with the Moho Pds. The phase from a depth of about 95 km is apparent in all the stacks, so there is good reason to conclude that these are real phases from a mid-lithosphere discontinuity or a shallow LAB. Figure 16. View largeDownload slide (a) Stacked PdsRF produced from WID processing of all 649 events recorded at ARU. (b) Stacked PdsRF produced using WID applied to ARU data from 2014 to 2015 (55 events), (c) Stacked PdsRF produced using WID applied to ARU data from 2012 (43 events). (d) Stack of traditional RF produced all 649 RFs computed from data recorded at ARU. Figure 16. View largeDownload slide (a) Stacked PdsRF produced from WID processing of all 649 events recorded at ARU. (b) Stacked PdsRF produced using WID applied to ARU data from 2014 to 2015 (55 events), (c) Stacked PdsRF produced using WID applied to ARU data from 2012 (43 events). (d) Stack of traditional RF produced all 649 RFs computed from data recorded at ARU. To determine how well WID may work for temporary broad-band deployments, which typically last for 1–2 yr, we preformed additional tests of WID for 1 to 2 yr subsets of the ARU data. The most useful product of WID appears to be the PdsRFs, so we only focused on Pds phases in our discussion of tests using 1 and 2 yrs’ worth of data. Fig. 16(a) shows the PdsRFs using all 649 events recorded from ARU data. Figs 16(b) and (c) were produced using 2 yr of data (55 events) and 1 yr of data (43 events), respectively. Fig. 16(d) is the Pds stack of all 649 RFs from ARU using traditional RF processing. All the phase specific PdsRFs recovered Pds phases from the expected depths (Moho, Mid-lithosphere, 410 and 660). Phases from deeper layers, such as 410 and 660 km, have sufficient curvature so that they do not mis-stack as PXds phases on the raw RFs (Fig. 16b). There are more mis-stacked phases (or noise) throughout the transition zone depths on the traditional RFs than are present on all three of the PdsRFs computed using WID (Figs 16a–c). The positive Pds phase between 160 and 180 km on the PdsRFs computed from all the ARU data is not present on the PdsRFs computed from 2 yr of data and have a negative polarity on the PdsRFs produced from 1 yr worth of data. This phase is likely real since it is strong on the PdsRFs produced from all the data, but the differences on the short term RFs may reflect azimuthal variations due to differences in data distribution for the three different PdsRFs. There are phases that are unique to the PdsRFs from the 1 or 2 yr data sets that are not common to all three WID PdsRFs. The difference between these PdsRFs may reflect limitations in the distribution for the 1 and 2 yr data set that result in the inclusion of different miss-stacked reverberations or azimuthal variations in the local geology. However, the point of this paper is not to analyse ARU geology but to demonstrate that WID can produce RFs that are relatively free of higher order reverberations from shallow layers. It is clear that all the PdsRFs computed for different amounts of data have less noise and fewer mis-stacked reverberations than the traditional RFs, which demonstrates that WID can be effective for short-term deployments. To determine if there are substantial interpretable phases in the residual waveform from WID (Fig. 17), we stacked RFs made from the residual of the seismograms left after WID RFs (Fig. 18). We see there are no distinctive phases in the stacked residual data, except possibly a small residual of the Moho Pds in Fig. 18(a). This possible residual Moho amplitude reflects that the extracted amplitude from iterative process may not completely remove all amplitude from an interface. There are no obvious phases remaining in the stacks of all three residual that stand out significantly above the noise level. Figure 17. View largeDownload slide Residual receiver function for station ARU after WID. Figure 17. View largeDownload slide Residual receiver function for station ARU after WID. Figure 18. View largeDownload slide Stacks of residual receiver function (a) using Pds moveout corrections, (b) PPds moveout corrections and (c) PSds moveout corrections. Figure 18. View largeDownload slide Stacks of residual receiver function (a) using Pds moveout corrections, (b) PPds moveout corrections and (c) PSds moveout corrections. 3 CONCLUSION WID produces phase specific RFs that are easier to interpret than the stacks of traditional RFs. The WID method is the most useful for shallow phases where the pulse width is close to or larger than the amount of moveout. To get an optimal result from WID, RFs from a wide range of distances are required. PdsRFs are the most valuable of the phase specific RFs because the PdsRFs tend to have very little contamination from higher order reverberation compared to the PSdsRFs and especially the PPdsRFs. WID also works best using data with high signal-to-noise ratios. We have tested WID at many stations and found that it did not work well for many of the South American stations. However, it produced clear 3-D stacked images using California data from permanent and temporary deployments with Moho and LAB phases (Ainiwaer & Gurrola, in preparation) that are similar to recent published Sp work (Lekic et al. 2011). We also found that WID works well using transportable array data from Texas (Gurrola et al., in preparation). The Pds stacked RFs from WID (PdsRFs) are the most improved, owing to the distinct shape of the Pds moveout curve as compared to the moveout curves for the PPds and PSds phases and reverberations. The H–k analysis in WID is only used to identify a phase, but arrival times and amplitudes for specific phases on the RFs are estimated using cross-correlation of each individual seismogram, independently. Therefore, the resulting RFs can still be used to model 2-D or 3-D Earth structure. Acknowledgements We appreciate the comments from fellow students to improve the manuscript. We are grateful the IRIS data centre for providing seismological data. We thank the reviewers (G. Pavlis and A. Frederiksen) for comments and suggestions that greatly improved this manuscript. REFERENCES Aki K., Richards P.G., 2002. Quantitative Seismology . 2nd edn, pp. 144– 145, University Science Books. Andrews J., Deuss A. 2008. Detailed nature of the 660 km region of the mantle from global receiver function data, J. geophys. Res.  113( B6), B06304, doi:10.1029/2007JB005111. https://doi.org/10.1029/2007JB005111 Google Scholar CrossRef Search ADS   Ammon C.J., 1991. The isolation of receiver effects from teleseismic P waveforms, Bull. seism. Soc. Am.  81 2504– 2510. 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. https://doi.org/10.1029/JB095iB10p15303 Google Scholar CrossRef Search ADS   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( B2), doi:10.1029/2011JB008560. https://doi.org/10.1029/2011JB008560 Dueker K.G., Sheehan A.F. 1998. Mantle discontinuity structure beneath the colorado rocky mountains and high plains, J. geophys. Res.  103( B4), 7153– 7169. https://doi.org/10.1029/97JB03509 Google Scholar CrossRef Search ADS   Dugda M.T. 2005. Crustal structure in Ethiopia and Kenya from receiver function analysis: implications for rift development in eastern Africa, J. geophys. Res.  110( B1), B01303, doi:10.1029/2004JB003065. https://doi.org/10.1029/2004JB003065 Google Scholar CrossRef Search ADS   Dziewonski A.M., Anderson D.L., 1981. Preliminary reference Earth model, Phys. Earth planet. Inter.  25( 4), 297– 356. https://doi.org/10.1016/0031-9201(81)90046-7 Google Scholar CrossRef Search ADS   Gao S.S., Silver P.G., Liu K.H. & the Kaapvaal Seismic Group, 2002. Mantle discontinuities beneath Southern Africa, Geophys. Res. Lett.  29( 10), doi:10.1029/2001GL013834. https://doi.org/10.1029/2001GL013834 Gilbert H.J., Sheehan A.F., Dueker K.G., Molnar P. 2003. Receiver functions in the western United States, with implications for upper mantle structure and dynamics, J. geophys. Res.  108( B5), 2229, doi:10.1029/2001JB001194. https://doi.org/10.1029/2001JB001194 Google Scholar CrossRef Search ADS   Gurrola H., Simmons N.A. 2000. Multiple seismic discontinuities near the base of the transition zone in the Earth's mantle, Nature  405( 6786), 559– 562. https://doi.org/10.1038/35014589 Google Scholar CrossRef Search ADS PubMed  Gurrola H., Minster J.B., Owens T. 1994. The use of velocity spectrum for stacking receiver functions and imaging upper mantle discontinuities, Geophys. J. Int.  117( 2), 427– 440. https://doi.org/10.1111/j.1365-246X.1994.tb03942.x Google Scholar CrossRef Search ADS   Gurrola H., Baker G.E., Minster J.B. 1995. Simultaneous time-domain deconvolution with application to the computation of receiver functions, Geophys. J. Int.  120( 3), 537– 543. https://doi.org/10.1111/j.1365-246X.1995.tb01837.x Google Scholar CrossRef Search ADS   Gurrola H., Bernard Minster J. 1998, Thickness estimates of the upper-mantle transition zone from bootstrapped velocity spectrum stacks of receiver functions, Geophys. J. Int.  133( 1), 31– 43. https://doi.org/10.1046/j.1365-246X.1998.1331470.x Google Scholar CrossRef Search ADS   Huang X., Gurrola H., Parker J., 2016. A comparative study of popular methods to compute receiver functions applied to synthetic data, Seism. Res. Lett. , 129–121–129–4. 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. https://doi.org/10.1046/j.1365-246x.2000.00217.x Google Scholar CrossRef Search ADS   Langston C.A., 1977. The effect of planar dipping structure on source and receiver responses for constant ray parameter, Bull. seism. Soc. Am.  67 1029– 1050. Langston C.A. 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves, J. geophys. Res.  84( B9), 4749– 4762. https://doi.org/10.1029/JB084iB09p04749 Google Scholar CrossRef Search ADS   Lawrance J.F., Shearer P.M. 2006. A global study of transition zone thickness using receiver functions, J. geophys. Res.  111( B6), n/a. https://doi.org/10.1029/2005JB003973 Lekic V., French S.W., Fischer K.M. 2011. Science , 334 783– 787. https://doi.org/10.1126/science.1208898 CrossRef Search ADS PubMed  Ligorria J.P., Ammon C.J., 1999. Iterative deconvolution and receiver-function estimation, Bull. seism. Soc. Am.  89 1395– 1400. Liu K.H., Gao S.S., Silver P.G., Zhang Y. 2003. Mantle layering across central South America, J. geophys. Res.  108( B11), 2510, doi:10.1029/2002JB002208, B11. https://doi.org/10.1029/2002JB002208 Google Scholar CrossRef Search ADS   Nair S.K., Gao S.S., Liu K.H., Silver P.G., 2006. Southern African crustal evolution and composition: constraints from receiver function studies, J. geophys. Res.  111( B2), doi:10.1029/2005JB003802. https://doi.org/10.1029/2005JB003802 Owens T.J., Crosson R.S., Hendrickson M.A., 1988. Constraints on the subduction geometry beneath western Washington from broadband teleseismic waveform modeling, Bull. seism. Soc. Am.  78 1319– 1334. Plešinger A., Hellweg M., Seidl D., 1986. Interactive high-resolution polarization analysis of broadband seismograms, J. Geophys.  59 129– 139. Tauzin B., Debayle E., Wittlinger G. 2008. The mantle transition zone as seen by global Pds phases: no clear evidence for a thin transition zone beneath hotspots, J. geophys. Res.  113( B8), B08309, doi:10.1029/2007JB005364. https://doi.org/10.1029/2007JB005364 Google Scholar CrossRef Search ADS   Vinnik L.P. 1977. Detection of waves converted from P to SV in the mantle, Phys. Earth planet. Inter.  15( 1), 39– 45. https://doi.org/10.1016/0031-9201(77)90008-5 Google Scholar CrossRef Search ADS   Wang Y., Pavlis G.L. 2016. Generalized iterative deconvolution for receiver function estimation. Geophys. J. Int.  204( 2), 1086– 1099. https://doi.org/10.1093/gji/ggv503 Google Scholar CrossRef Search ADS   Zhu L., Kanamori H. 2000. Moho depth variation in southern California from teleseismic receiver functions, J. geophys. Res.  105( B2), 2969– 2980. https://doi.org/10.1029/1999JB900322 Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Wavefield iterative deconvolution to remove multiples and produce phase specific Ps receiver functions

, Volume 212 (3) – Mar 1, 2018
15 pages

/lp/ou_press/wavefield-iterative-deconvolution-to-remove-multiples-and-produce-O5cBArogZE
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx454
Publisher site
See Article on Publisher Site

### Abstract

Summary Common conversion point stacking or migration of receiver functions (RFs) and H–k (H is depth and k is Vp/Vs) stacking of RFs has become a common method to study the crust and upper mantle beneath broad-band three-component seismic stations. However, it can be difficult to interpret Pds RFs due to interference between the Pds, PPds and PSds phases, especially in the mantle portion of the lithosphere. We propose a phase separation method to isolate the prominent phases of the RFs and produce separate Pds, PPds and PSds ‘phase specific’ receiver functions (referred to as PdsRFs, PPdsRFs and PSdsRFs, respectively) by deconvolution of the wavefield rather than single seismograms. One of the most important products of this deconvolution method is to produce Ps receiver functions (PdsRFs) that are free of crustal multiples. This is accomplished by using H–k analysis to identify specific phases in the wavefield from all seismograms recorded at a station which enables development of an iterative deconvolution procedure to produce the above-mentioned phase specific RFs. We refer to this method as wavefield iterative deconvolution (WID). The WID method differentiates and isolates different RF phases by exploiting their differences in moveout curves across the entire wave front. We tested the WID by applying it to synthetic seismograms produced using a modified version of the PREM velocity model. The WID effectively separates phases from each stacked RF in synthetic data. We also applied this technique to produce RFs from seismograms recorded at ARU (a broad-band station in Arti, Russia). The phase specific RFs produced using WID are easier to interpret than traditional RFs. The PdsRFs computed using WID are the most improved, owing to the distinct shape of its moveout curves as compared to the moveout curves for the PPds and PSds phases. The importance of this WID method is most significant in reducing interference between phases for depths of less than 300 km. Phases from deeper layers (i.e. P660s as compared to PP220s) are less likely to be misinterpreted because the large amount of moveout causes the appropriate phases to stack coherently if there is sufficient distribution in ray parameter. WID is most effective in producing clean PdsRFs that are relatively free of reverberations whereas PPdsRFs and PSdsRFs retain contamination from reverberations. Composition and structure of the continental crust, Composition and structure of the mantle, Body waves, Crustal imaging 1 INTRODUCTION Receiver function (RF) techniques are widely used to make 2-D and 3-D stacked images of crust and upper mantle structure beneath networks of three component broad-band seismic stations (Dueker & Sheehan 1998; Simmons & Gurrola 2000; Gilbert et al. 2003; Liu et al. 2003; Gao et al. 2002; Andrews & Deuss 2008). Single station analysis of crust and upper mantle structure is commonly performed by inversion of RF waveforms (Ammon et al. 1990; Julia et al. 2000; Bodin et al. 2012), H–k analysis to estimate Moho depths and Vp/Vs (Zhu & Kanamori 2000; Dugda et al. 2005; Nair et al. 2006) and studies of the upper mantle transition zone (Gurrola et al. 1994; Gurrola & Minster 1998; Gilbert et al. 2003; Lawrance & Shearer 2006; Tauzin et al. 2008). A problem with all of these methods is that there can be interference (Fig. 1) that can complicate the interpretation of these data. Figure 1. View largeDownload slide Illustration of ray paths for P, Pds, PPds and PSds phases. Figure 1. View largeDownload slide Illustration of ray paths for P, Pds, PPds and PSds phases. In RF analysis, multiple traces can be grouped and stacked after correcting for moveout (time delays at different ray parameters). The arrival times of the different Ps phases are delayed with increasing ray parameter, whereas the moveout for reverberations (i.e. PPds and PSds phases) is advanced with increasing ray parameter. For 1-D structure, these moveout times are easily computed (Vinnik 1977). By stacking moveout corrected RFs from a range of ray parameters using a range of reasonable velocity models, one can estimate the depth and average velocity structure to the boundary via velocity spectrum stacking (VSS, Gurrola et al. 1994; Gurrola & Minster 1998) or estimate depth and Vp/Vs ratio (H–k method, Zhu & Kanamori 2000) to these boundaries. However, in all these traditional Ps RF methods, reverberations from shallow depths (e.g. mid-crust and Moho) can be mis-stacked as Ps phases or interfere with Ps phases from greater depths. To avoid misinterpretation of phases and to better stack RFs, we developed a method to separate Pds, PPds and PSds phases from RFs and produce separate, single phase type receiver functions (PdsRFs, PPdsRFs and PSdsRFs) from each seismogram. The wavefield iterative deconvolution (WID) method is introduced to differentiate and isolate different RF phases by exploiting differences in moveout curves for each of the common phase types. This is a phase separation method to produce separate, ‘phase specific’ Pds, PPds and PSds receiver functions. We proceed by first illustrating the potential effectiveness of the WID method using synthetic RFs computed using raytracing (Aki & Richards 2002). Then, we apply the WID technique to data from the broad-band seismic station at Arti, Russia (ARU). We chose ARU as a test station because this station has high quality data from a wide range of ray parameters and the crust and mantle structure there is relatively uncomplicated as compared to other stations around the world. 2 METHOD Ps RF analysis is a method to isolate Pds converted phases and reverberations of the P-phase that end in an S-phase that originate from the structure beneath a three-component seismic station. Since Pds conversions are SV phases, we separate the SV and SH phases by rotating these data from the N--S and E–W components to compute radial and transverse components (Plesinger et al. 1986). Deconvolving the radial and transverse components by the vertical component produces RFs that isolate the Pds phases and reverberations ending in an S phase (Lang 1979; Ammon 1991). Frequency domain deconvolution by spectral division of the horizontal components by the vertical component after pre-whitening is one of the oldest methods of deconvolution used to produce RFs (Langston 1977; Owens et al. 1988). The time-domain iterative deconvolution method, introduced to RF analysis by Ligorria & Ammon (1999), lends itself well to analysis that incorporates the shape of the moveout curves to distinguish Pds, PPds and PSds phases from one another. We accomplish this by including all the seismograms recorded at a seismic station into a simultaneous deconvolution that uses a modified H–k analysis to identify the wavefield of the strongest remaining phase in each iteration and construct three, phase specific RFs from each seismogram in the data set. We adopt the regularized iterative deconvolution method (Wang & Pavlis 2016) that employs a regularized deconvolution of the seismogram, rather than cross correlation, to identify phases in each iteration of the deconvolution. Wang & Pavlis (2016) give a detailed discussion of the benefits of the regularized iterative deconvolution over the use of cross-correlation. In a head-to-head comparison of five common methods of computing RFs, Huang et al. (2016) also found the regularized iterative deconvolution to be the most robust of the methods tested. Fig. 1 illustrates the ray paths of converted Pds phases and major reverberations ending in converted S-phases (PPds and PSds). The relative arrival times of any phase in a RF for a horizontal interface is the arrival time of the phase of interest minus the arrival time of the direct P phase (Langston 1977; Vinnik 1977; Gurrola et al. 1994) and, for 1-D velocity models, can be represented by the following equations:   $$\Delta {T_{Ps}} = Z*\left( {\sqrt {{V_s}^{ - 2} - {p^2}} - \sqrt {{V_p}^{ - 2} - {p^2}} } \right)$$ (1)  $$\Delta {T_{PPds}} = Z*\left( {\sqrt {{V_s}^{ - 2} - {p^2}} + \sqrt {{V_p}^{ - 2} - {p^2}} } \right)$$ (2)  $$\Delta {T_{PSds}} = 2*Z\left( {\sqrt {{V_s}^{ - 2} - {p^2}} } \right),$$ (3)where p is the ray parameter of the incident wave, Vs is the S-wave velocity, and Z is the depth to the interface. In RF work, Vp is typically replaced by k*Vs (where k = Vp/Vs) in eqs (1) and (2) so that we can compute traveltimes from a Vs model and an assumed (or estimated) value for k. By writing PPds and PSds traveltime equations (eqs 2 and 3) as a function of k and Vs rather than as a function of Vp and Vs, we can simply use one equation for both types of reverberations and use k as the discriminating criteria between PPds and PSds phases:   $$\Delta {T_{PXds}} = Z*\left( {\sqrt {{V_s}^{ - 2} - {p^2}} + \sqrt {{{(k*{V_s})}^{ - 2}} - {p^2}} } \right).$$ (4)In this case, k = 1 will correspond to PSds phases and eq. (4) would reduce to be the same as eq. (3). The WID algorithm described below identifies Pds, PPds and PSds phases from single-phase H–k analysis to produce phase specific RFs. To perform WID we: Calculate two ‘single phase’ H–k matrices, one for Pds phases and the other for PXds phases using traveltime eqs (1) and (4), respectively. Any reasonable 1-D S-wave velocity model can be used as a reference model in H–k stacking. To get optimal results from the H–k analysis, it is necessary to have data from several ray parameter bins. The algorithm searches for the largest amplitude across the two H–k matrices and identifies the phase by:   \begin{equation*} Phase\ type = \left( {\begin{array}{@{}*{1}{c}@{}} {\rm Pds \ if \ largest \ phase \ is \ on \ the \ Pds} \ H - k \ {\rm matrix}\\ {\rm PPds \ if \ largest \ phase \ is \ on \ the \ PXds} \ H - k \ {\rm matrix \ and} \ 1.5 < k < 1.95\\ {\rm PSds \ if \ largest \ phase \ is \ on \ the \ PXds} \ H - k \ {\rm matrix \ and} \ 0.9 < k < 1.1 \end{array}} \right)\end{equation*} We then use cross-correlation over a narrow time window centred around the expected arrival time of this phase to estimate the ‘real’ arrival time on each seismogram. To assure that the arrival times and amplitudes for the identified phase on each of the phase specific RFs reflects the original data, we only use steps one and two to estimate the depth and k value for each phase to identify its phase type. We then construct 3 phase specific RFs (spike functions) for each seismogram by placing the amplitude of the identified phase at the time delays (determined in step 3) on the phase specific RFs. The phase specific RFs are convolved by the vertical component to produce the three different, phase dependent seismograms. All three of the phase dependent seismograms are subtracted from the observed seismograms to produce modified ‘observed’ seismograms with all previously identified phases removed (Ligorria & Ammon 1999). Steps 1 through 6 are repeated for many iterations until all significant phases are removed from the original seismograms. The number of iterations needed to complete this process is determined by calculating the sum of the standard deviations of the modified seismograms during each iteration. The algorithm stops when sum of the standard deviations between iterations becomes small or nearly the same. 2.1 Test of WID with synthetic data To test the WID method described above, we applied WID to synthetic seismograms computed for the teleseismic range of ray parameters (0.04–0.08 s km–1, Fig. 2). These synthetics were computed by ray tracing through the modified PREM velocity model (Dziewonski & Anderson 1981). We modified PREM by adding a 2 km thick sedimentary surface layer to simulate complications that may occur at many seismic stations. We also added a 75 km deep low velocity layer to simulate arrivals that are of interest in recent RF investigations of the lithosphere. The amplitude of each phase was determined by reflection and transmission coefficient equations from Aki & Richards (2002) and a 0.5 Hz Gaussian low-pass filter was applied. Figure 2. View largeDownload slide Synthetic receiver function as a function of ray parameter computed from a modified version of the PREM velocity model. We labelled the prominent Pds phases from the Moho, the 75 km low velocity layer (that we added to the velocity model), the 220, 410 and 660 km discontinuities. The phases from the surface layer and Mid-crustal discontinuities were not labelled so as not to clutter the figure. Reverberations from Moho and higher order reverberations are also labelled. Figure 2. View largeDownload slide Synthetic receiver function as a function of ray parameter computed from a modified version of the PREM velocity model. We labelled the prominent Pds phases from the Moho, the 75 km low velocity layer (that we added to the velocity model), the 220, 410 and 660 km discontinuities. The phases from the surface layer and Mid-crustal discontinuities were not labelled so as not to clutter the figure. Reverberations from Moho and higher order reverberations are also labelled. The Pds phases from the surface layer, mid-crust, 75 km low velocity layer (LVL), Moho, 220, 410 and 660 km discontinuities are clearly visible in these synthetics (Fig. 2). Pds traveltimes increase as ray parameter increases, which is opposite to the moveout of the PXds phases that arrive earlier as ray parameter increases. The curvature of Pds traveltimes (moveout) becomes larger as conversion depth increases. For example, it is harder to see the moveout of the mid-crust and Moho Pds phases than for the Pds phase from deeper layers such as P220s, P410s and P660s. The difference in the shape of the moveout curve between Pds and PXds phases can be used to distinguish them from one another. For example, notice that the P660s and PP220s phases curve in opposite directions. We produced ‘single phase’ H–k plots (Pds and PXds H–k plots) using the synthetic data shown in Fig. 2. Fig. 3(a) is the single phase H–k plot computed by stacking these RFs using the Pds traveltime equation (eq.1). We refer to this as a single phase H–k plot because it is a plot of stacked Pds phase amplitudes as a function of H and k values and should not be confused with the H–k method by Zhu & Kanamori (2000) that calculates the summation of the weighted stack of Pds, PPds and PSds phases. The largest stacked amplitudes are for Pds phases from the near surface layer, Moho, 220, 410 and 660 km discontinuities and are centred near k = 1.81 (marked by stars). However, the peaks corresponding to the Pds phases from mid-crust and 75 km LVL appear at a k value of ∼1.77. Figure 3. View largeDownload slide H–k plot using Pds moveout curve (a). Pds phases from surface layer, Mid-crust, Moho, 75 km low velocity layer, 220, 410 and 660 discontinuities are labelled; H–k plot using PXds moveout curve (b). Mid-crust, Moho and 220 PPds, Mid-crust and Moho PSds are labelled. Stars represent location of highest amplitude for each phase. Figure 3. View largeDownload slide H–k plot using Pds moveout curve (a). Pds phases from surface layer, Mid-crust, Moho, 75 km low velocity layer, 220, 410 and 660 discontinuities are labelled; H–k plot using PXds moveout curve (b). Mid-crust, Moho and 220 PPds, Mid-crust and Moho PSds are labelled. Stars represent location of highest amplitude for each phase. We see in Fig. 3 that, even for noise free synthetics, the peak for the largest amplitude is elongated and not sharp. There can be large errors in picking this peak, but in WID the estimated value for k only needs to be accurate enough to distinguish between the three phase types. If the phase is a PPds or PSds phase, it will have its largest peak along the left edge of the Pds H–k plot (Fig. 3a); otherwise it will be near the correct value. Likewise, a Pds phase will plot on the left edge of the PXds H–k plot. Thus, it is unlikely that errors in the estimated k value will cause a Pds phase to be incorrectly identified as a PXds phase. Highly accurate values for k are also not essential to distinguish between the PSds and PPds phases on the PXs H–k plot (Fig. 3b), since the PSds phase is identified by a k value near 1, whereas the PPds is identified by having a value greater than 1.5. As a result, an error analysis of the k value is not essential because our method is effective even in the presence of relatively large errors in the estimated value for k. Fig. 3(b) is the H–k plot computed by stacking RFs using the PXds traveltime equation (eq. 4). To better understand how each phase is identified and placed on the appropriate phase specific RF, we produced plots of the amplitudes, H and k values for each pulse added to each of the phase specific RFs (Fig. 4). Each circle in Fig. 4(a) represents the H and k values that correspond to the largest stacked amplitude of a Pds phase from each iteration of the WID. Fig. 4 shows that many of the peaks on these amplitude versus H–k plots do not match the value of 1.8 used to produce the synthetics in Fig. 2. These errors in k values are due to interference between different phases (i.e. interference between the PPds phase and higher order reverberations from the surface layer). A k value of 1.81 was found for the Moho and 220 PPds phases, but the largest amplitude of the mid crust PPds phase occurred at k = 1.6. The peak amplitude of the PSds phase for the mid crust and Moho is between k = 0.9 and k = 1 in Figs 4(b) and (c). Errors in the estimated k values for the PPds and PSds phases from the surface layer (k 1.5 and k 1.06, respectively) do not prevent them from being identified as the appropriate phase. In general, we expect the errors for k to be larger for shallow layers since they have the flattest moveout curves. Figure 4. View largeDownload slide H & k value and amplitudes for phases identified in each iteration of the phase separation portion of the WID. The relative amplitude for the phase is illustrated by the size of the circle. Separate plots (a, b and c) correspond to the Pds, PPds and PSds, respectively. Figure 4. View largeDownload slide H & k value and amplitudes for phases identified in each iteration of the phase separation portion of the WID. The relative amplitude for the phase is illustrated by the size of the circle. Separate plots (a, b and c) correspond to the Pds, PPds and PSds, respectively. We applied the WID method discussed above to the synthetics to produce three separate phase dependent RFs (PdsRFs, PPdsRFs and PSdsRFs are shown in Figs 5a–c respectively). Fig. 5(d) is the residual of the original seismograms after subtracting all the phase specific RFs from the original synthetics. The Pds phases from the surface layer, mid crust, Moho, 75 km LVL, 220, 410 and 660 km are clearly visible on the PdsRFs (Fig. 5a). The PPds phases from the mid crust, Moho and 220 are the strongest phases on PPdsRFs but many higher order reverberations of these phases from shallower layers are present as well (Fig. 5b). The PPds phase from 75 km LVL is present and strong but is hidden among the higher order reverberations from shallower layers. The PPdsRFs are contaminated with more high order reverberations than PSdsRFs. The PdsRFs are the least contaminated by other phases. The residual RFs (Fig. 4d) show that there are some artefacts in the data after removing the strong RF phases, but these artefacts would likely be below the noise level of real data. Figure 5. View largeDownload slide Receiver functions produced from the synthetics in Fig. 2 after phase separation: (a) Pds receiver function (b) PPds receiver function (c) PSds receiver function and (d) Residual receiver function Figure 5. View largeDownload slide Receiver functions produced from the synthetics in Fig. 2 after phase separation: (a) Pds receiver function (b) PPds receiver function (c) PSds receiver function and (d) Residual receiver function The utility of WID is best illustrated by a comparison of the stacks of the phase specific RFs after applying moveout corrections. The stack of the PdsRFs (Fig. 6 trace b) only has Pds phases from the surface layer, mid-crust, Moho, 75 km LVL, 410 and 660 km discontinuities, whereas the stacked RF produced from traditional RFs (Fig. 6 trace a) is contaminated by mis-stacked PPds phases from the mid-crust and Moho that that appear at 90 and 140 km depths. Likewise, PSds phases from the shallow layer, mid-crust and Moho appear as mis-stacked negative phases at about 5, 110 and 190 km depths. The higher order reverberations are mis-stacked as phases between the 250 and 350 km depths. The PPds phase from the 220 has sufficient moveout so as not to ‘mis-stack’ on any of these stacked RFs. So, WID is most effective for differentiating phases from shallower depths. Figure 6. View largeDownload slide Pds stacks of the synthetic RF produced from the noise free synthetics without WID phase separation (a) and with WID phase separation (b). Also shown are stacked Pds RFs computed from noisy synthetic without WID phase separation (c) and with WID phase separation (d). Figure 6. View largeDownload slide Pds stacks of the synthetic RF produced from the noise free synthetics without WID phase separation (a) and with WID phase separation (b). Also shown are stacked Pds RFs computed from noisy synthetic without WID phase separation (c) and with WID phase separation (d). To test the robustness of WID in the deconvolution of noisy data, we added 10 per cent noise to the synthetic data and repeated the above tests of the WID method (Fig. 6). The noise added to the synthetics is the residual after WID of the data from the seismic station ARU (see the discussion below) after scaling the residual noise to have a standard deviation that is 10 per cent of the peak signal in the synthetics. The stacked PdsRFs computed with added noise (Fig. 6d) and without added noise (Fig. 6d) are nearly identical. Fig. 7 shows the stacked RFs after applying PPds moveout corrections and converting to depth. Fig. 7(a) is the stack of traditional RFs, and traces b, c and d are stacks of the PPdsRFs produced using WID with different criteria. Once again, the stacked traditional RFs (Fig. 7a) are contaminated by mis-stacked phases, whereas the phase separated PPdsRFs (Figs 7b–d) have much less contamination from other mis-stacked phases. The shallowest phases on the traditional PPds RF (Fig. 7a) are mis-stacked Pds phases. These shallow phases have very little moveout so they can appear strongly on the wrong stacked RFs. There are several phases in the 50–200 km depth ranges that have very small amplitudes. These are either higher order reverberations of the PXds phases or mis-stacked Pds phases from greater depths. By changing the phase identification criteria in the WID, we may be able to clean up the PPdsRFs. By narrowing the range of allowable k values to 1.75 through 1.85, the stacked PPdsRF (Fig. 7c) contains fewer higher order reverberations, but it is also missing some of the expected PPds phases (i.e. the PPds phases from the mid-crustal and the 75 km LVL). The stack of PPdsRFs with only the 10 largest PPds phases (Fig. 7d) identified by WID includes more of the expected phases but also retains some of the large amplitude higher order reverberations. From this analysis, it is obvious that WID is not as effective in producing ‘clean’ PPdsRFs as it is in producing a relatively uncontaminated stacked PdsRFs. An argument can be made that WID is only effective in producing clean PdsRFs, but this analysis shows that WID can be used to extract useful information from the PPds phase and produce a better stacked PPdsRFs than traditional processing. The PPdsRFs produced from the same noisy synthetics used in the analysis of the PdsRFs above have a few small phases not present in the respective noise free phase specific stacked PPdsRFs, but these are very small peaks and may be the result of reverberations that were not completely removed from the ARU residuals. Figure 7. View largeDownload slide Stacked PPds receiver functions produced from the synthetics in Fig. 2. Traces a and e are the stacked RF produced from traditional RFs without noise and with real seismic noise added, respectively. Traces b, c, d and f are stacked, phase specific PPdsRFs produced using the WID method and noise free PPdsRFs including phases with k values between 1.5 and 1.9 (b), noise free PPdsRFs including phases with k values between 1.75 and 1.85 (c), noise free PPdsRFs with only the phases with the 10 largest amplitudes (d), and PPdsRFs produced from synthetics with real seismic noise added and including phases with k values between 1.5 and 1.9 (f). Figure 7. View largeDownload slide Stacked PPds receiver functions produced from the synthetics in Fig. 2. Traces a and e are the stacked RF produced from traditional RFs without noise and with real seismic noise added, respectively. Traces b, c, d and f are stacked, phase specific PPdsRFs produced using the WID method and noise free PPdsRFs including phases with k values between 1.5 and 1.9 (b), noise free PPdsRFs including phases with k values between 1.75 and 1.85 (c), noise free PPdsRFs with only the phases with the 10 largest amplitudes (d), and PPdsRFs produced from synthetics with real seismic noise added and including phases with k values between 1.5 and 1.9 (f). Fig. 8 shows the stack of PSdsRFs from Fig. 5(c) and the stack of traditional RFs. Once again, the PSdsRFs (Fig. 8b) only show legitimate, strong PSds phases from the mid-crust, Moho and 75 km LVL. The traditional stacked PSdsRFs (Fig. 8a) have strong phases at shallow depths that are mis-stacked Pds and PPds phases and include the poorly stacked Pds phases from greater depths in the 130–180 km depth interval. The PSdsRFs (Fig. 8d) produced using the noisy synthetics are relatively unaffected by the noise. Figure 8. View largeDownload slide Stacked PSds receiver functions produced from the synthetics in Fig. 2. Traces a and c are the stacked RF produced from traditional RFs without noise and with real seismic noise added, respectively. Traces b and d are stacked, phase specific PSdsRFs produced using the WID method with noise free PSdsRFs (b) and PSdsRFs produced from synthetics with real seismic noise added. Figure 8. View largeDownload slide Stacked PSds receiver functions produced from the synthetics in Fig. 2. Traces a and c are the stacked RF produced from traditional RFs without noise and with real seismic noise added, respectively. Traces b and d are stacked, phase specific PSdsRFs produced using the WID method with noise free PSdsRFs (b) and PSdsRFs produced from synthetics with real seismic noise added. WID is effective in selecting the expected phase type for each of the phase specific stacked RFs computed from synthetic data. In the case of the PPdsRFs and PSdsRFs, higher order reverberations of the expected phases are still present, but the PdsRFs are free of such reverberations. The importance of this WID method is most significant in cleaning up the RFs at depths of less than 300 km. The deeper converted phases are less likely to be misinterpreted in stacked RFs because there is less cross stacking between Pds and PXds phases since the amount of moveout is much greater than the pulse width of these deeper phases. However, PPds and PSds phases can cross contaminate each other to greater depths because the moveout curves are similar in shape. 2.2 Test of WID using real data We used data from the GDSN station ARU (Arti, Russia) to test the utility of the WID method applied to real data. We analysed all the data recorded from 1990 through 2016 at ARU and found 649 events that produced RFs with high signal-to-noise ratios (Fig. 9). The epicentral distance for this data set ranges from 30° to 90° (corresponding to a range of ray parameters of 0.04 to 0.08 s km–1). The backazimuths with data coverage across a wide range of ray parameters are mainly from 80° to 140°. Data were windowed from 30 s prior to and 150 s after the direct-P wave arrival. A data stream of this length is sufficient to include all primary reverberations from lithospheric depths and Pds phases through the upper mantle transition zone. We rotated these seismograms into the L-Q-T coordinate system after estimating the incident angle of incoming P wave using singular value decomposition to find the angle that maximizes the amplitude of P wave energy on the Q component. By doing so, we maximize the isolation of SV and P-wave energy on the L and Q component, respectively (Plesinger et al. 1986). An advantage of using the L-Q-T system is that it makes it possible to detect phases from shallow depths by avoiding side lobes from the direct P. By minimizing the contamination of the P phase (Q component) with SV energy, we also have a cleaner source function for deconvolution. Figure 9. View largeDownload slide Station and earthquakes epicentres recorded from 1990 to 2016 used in our study. Black triangle shows ARU station location; black circles show the 649 earthquakes used in this study. Figure 9. View largeDownload slide Station and earthquakes epicentres recorded from 1990 to 2016 used in our study. Black triangle shows ARU station location; black circles show the 649 earthquakes used in this study. We used simultaneous frequency domain water level deconvolution to produce the unfiltered SV RFs by deconvolving the Q component from the L component of the data. Prior to deconvolution, we applied a bandpass filter to the data using corner frequencies of 0.02 and 0.5 Hz. We clustered the events into ray parameter bins 0.001 s km–1 wide to enable simultaneous deconvolution of all the events in a bin (Gurrola et al. 1995). Fig. 10 shows the traditional RFs computed for the ARU station. These RFs were shifted so the first arrival would be at 10 s. The dominant phases are Pds phases from the surface layer and the Moho, which appear at 11 and 14 s. Other phases are present but they are masked by the noise. In the next section, we apply the WID method to these RFs to produce PdsRFs, PPdsRFs and PSdsRFs. Figure 10. View largeDownload slide Receiver functions computed for ARU station with 0.02–0.5 Hz bandpass filter applied. Surface layer and Moho Pds phases are the only clearly observable phases. Figure 10. View largeDownload slide Receiver functions computed for ARU station with 0.02–0.5 Hz bandpass filter applied. Surface layer and Moho Pds phases are the only clearly observable phases. The number of iterations included in the WID was determined by computing the sum of the standard deviation of the residual seismograms after each iteration of the WID (Fig. 11), and the WID was stopped when the sum of the standard deviation changed by very little between iterations. In our case, 120 iterations appeared to be a sufficient number of iterations. Figure 11. View largeDownload slide Sum of the standard deviations of all receiver function traces versus number of iteration through iterative deconvolution. Figure 11. View largeDownload slide Sum of the standard deviations of all receiver function traces versus number of iteration through iterative deconvolution. Fig. 12(a) is a ‘single phase’ H–k plot computed by stacking RFs computed from ARU data using Pds traveltime equations for a range of k values using PREM S-velocities as the reference model. The largest stacking amplitude occurs for the Moho at ∼38 km in depth with k = 1.87. The mid lithosphere at 95 km depth, 410 and 660 Pds phases are also identified at k = 1.87, 1.88 and 1.89, respectively. The k values found for ARU data are different than those found for the PREM synthetics. PREM velocities are not an exact match to the velocity model appropriate for ARU so we would not expect accurate k values. However, k is only used here as a discriminant to identify phase types; any reasonable 1-D velocity model works for WID. Figure 12. View largeDownload slide Single phase H--k plot using Pds moveout curve for station ARU receiver function (a) and the H–k plot using PXds moveout curve (b). Pds phases from Moho, mid-lithosphere, 410 and 660 discontinuities are labelled; Moho PPds and PSds phases are also labelled. Figure 12. View largeDownload slide Single phase H--k plot using Pds moveout curve for station ARU receiver function (a) and the H–k plot using PXds moveout curve (b). Pds phases from Moho, mid-lithosphere, 410 and 660 discontinuities are labelled; Moho PPds and PSds phases are also labelled. To verify that the WID is not dependent on variations in a reasonable velocity model, we preformed WID on the synthetics using a modified version of PREM (changing PREM by ±10 per cent) and found that stacked phase specific RFs were virtually identical regardless of which velocity model was used. Fig. 12(b) is an H–k plot computed by stacking RFs using the PXds traveltime equation (eq. 4). We find the maximum amplitude for the Moho PPds phase is at k = 1.5, and the peak amplitude for Moho PSds occurs at k = 0.95 (Figs 12b, 14b and c). These H–k plots appear to have more phases than expected for the number of layers in our velocity model. These extra phases are stacked higher order reverberations whose amplitude are smaller than expected phases. However, the expected PPds and PSds phases are clearly observable on the phase specific RFs. Fig. 13 shows the PdsRFs (a), PPdsRFs (b) and PSdsRFs (c) produced by WID for the ARU data. We are able to identify Pds phases from a shallow surface layer, the Moho, a mid-lithosphere discontinuity (at 95 km depth), the 410 discontinuity and the 660 km discontinuity in the PdsRFs (Fig. 13a). All these phases are clearly delayed as ray parameter increases, which is expected for the Pds phase. The unlabelled phases at around 20 s are higher order Pds reverberations from shallower layers. The amplitude versus offset (AVO) relations show the expected pattern of increasing amplitude with larger ray parameter. Figure 13. View largeDownload slide Individual phase receiver function produced from ARU data after phase separation (a), (b) and (c) shows the Pds, PPds and PSds receiver function, respectively. Figure 13. View largeDownload slide Individual phase receiver function produced from ARU data after phase separation (a), (b) and (c) shows the Pds, PPds and PSds receiver function, respectively. Fig. 14 shows depths, amplitudes and k values of all the phases in the RFs produced using WID on the ARU data. Fig. 14(a) shows the H and k values used to construct the PdsRFs in Fig. 13(a). All the dots for the Psd phase in Fig. 14(a) are at depths consistent with expected discontinuities, but the k values vary, reflecting errors due to the widths of the peaks on the H–k plots, errors in the assumed velocity model, contamination by noise and interference from other phases such as higher order reverberations. The multiple k values for the same depth (H), such as the Moho depth and the surface layer, reflect that, in WID, the same depth can be revisited between iterations. It is noteworthy that largest amplitudes found for the Moho Pds (Fig. 14a), PPds (Fig. 14b) and PSds (Fig. 14c) phases are at a depth of ∼38 km. The k value of 1.88 found for the P410s and P660s phases is higher than computed from PREM (which is 1.85 for 410 and 1.84 for 660), but in WID we only use the k for phase identification and not to interpret the local geology. Errors in the assumed velocity model will leak into the errors in the estimated k. However, our goal here is to clean up RFs to produce individual phase specific PdsRFs, PPdsRFs and PSdsRFs in the time domain, not to estimate k. The PPdsRFs (Fig. 13b) appear to have strong phases from the Moho (at 15 s) and mid lithosphere (at 35 s), but, as with the synthetics, the PPdsRFs are more heavily contaminated by higher order reverberations than the PdsRFs or the PSdsRFs. The range of k values for the PPds phases (Fig. 14b) is larger than found for the Pds phases, because the PPds plot includes many higher order reverberations, which should not stack at the correct k value. We can suppress reverberations by using a narrower range of k values, but we may not want to do so, because these reverberations may be useful in waveform modelling. Figure 14. View largeDownload slide H & k value for each iteration during phase separation (a), Pds receiver function amplitudes, (b) PPds receiver function amplitudes and (c) PSds receiver function amplitudes. Figure 14. View largeDownload slide H & k value for each iteration during phase separation (a), Pds receiver function amplitudes, (b) PPds receiver function amplitudes and (c) PSds receiver function amplitudes. The Moho and the mid lithosphere phases in the PSdsRFs (Fig. 13c) are easily identified since they have reversed amplitudes compared to PdsRFs and PPdsRFs. The k values for PSds are centred around a k 1 (Fig. 14c), which is consistent with what we expected for PSds phases in setting up the PXds phase moveout eq. (4). Fig. 15 shows the Pds, PPds and PSds produced by WID and traditional RF methods. The Moho and mid lithosphere Pds phases are separated with little interference from other phases in the PdsRFs (Fig. 15a). There appears to be two phases from the Moho depth on the traditional Pds RF (Fig. 15b), but one of the two is likely a mis-stacked PPds phase from a shallower layer. In the PdsRFs (Fig. 15a), the shallower of the two phases were removed. Moho PPds and PSds are mis-stacked at depths of about 150 and 180 km, respectively, in the traditional Pds RFs (Fig. 15b), but these phases were removed from the PdsRFs (Fig. 15b) computed using WID. Figure 15. View largeDownload slide Stacks of receiver functions using all 649 events recorded at ARU for depths of up to 200 km: (a) stack of PdsRFs produced using WID, (b) stack of Pds RFs using traditional RFs, (c) stack of WID PPdsRF only including the phases with the eight largest amplitudes, (d) stack of WID PPdsRFs only including those phases identified with Vp/Vs ratios of between 1.85 to 1.89, (e) stack of WID PPdsRFs that includes phases identified with Vp/Vs ratios between 1.4 and 1.9, (f) stack of traditional PPds RFs, (g) stack of PSdsRFs produced using WID, and (h) stack of traditional PSds RFs. Note that Moho and mid-lithosphere phases are labelled since they are shown in all three stacks. Moho phases stack coherently at a depth of ∼38 km and a mid-lithosphere discontinuity is at 95 km depth. Figure 15. View largeDownload slide Stacks of receiver functions using all 649 events recorded at ARU for depths of up to 200 km: (a) stack of PdsRFs produced using WID, (b) stack of Pds RFs using traditional RFs, (c) stack of WID PPdsRF only including the phases with the eight largest amplitudes, (d) stack of WID PPdsRFs only including those phases identified with Vp/Vs ratios of between 1.85 to 1.89, (e) stack of WID PPdsRFs that includes phases identified with Vp/Vs ratios between 1.4 and 1.9, (f) stack of traditional PPds RFs, (g) stack of PSdsRFs produced using WID, and (h) stack of traditional PSds RFs. Note that Moho and mid-lithosphere phases are labelled since they are shown in all three stacks. Moho phases stack coherently at a depth of ∼38 km and a mid-lithosphere discontinuity is at 95 km depth. As discussed for the synthetic tests above, the phase separated PPdsRF stacks have more peaks than PdsRFs (Fig. 16a) and PSdsRFs (Fig. 16g) when phases from a wide range of k are included (Fig. 15c). When we further restrict the range of allowable k values (Fig. 15d) and amplitudes (Fig. 15c) in phase selection during WID, we can produce cleaner PPdsRFs. Fig. 15(c) shows the PPdsRFs with only the eight largest amplitude arrivals included, resulting in PPdsRFs with phases from the Moho, but without the mid lithosphere phase which has a small amplitude. By using a tighter range of k values (Fig. 15d), we are able retain the phases from the mid-lithosphere and an apparent double phase from the Moho depth. However, this PPdsRF has less contamination from higher order reverberations than are visible in the PPdsRF computed using a wider range of k values. It is also possible that some low amplitude phases are side lobes or artefacts leaked in during the WID. Nevertheless, the Moho PPds and PSds phases are stacked at a depth of about 38 km, which is consistent with the Moho Pds. The phase from a depth of about 95 km is apparent in all the stacks, so there is good reason to conclude that these are real phases from a mid-lithosphere discontinuity or a shallow LAB. Figure 16. View largeDownload slide (a) Stacked PdsRF produced from WID processing of all 649 events recorded at ARU. (b) Stacked PdsRF produced using WID applied to ARU data from 2014 to 2015 (55 events), (c) Stacked PdsRF produced using WID applied to ARU data from 2012 (43 events). (d) Stack of traditional RF produced all 649 RFs computed from data recorded at ARU. Figure 16. View largeDownload slide (a) Stacked PdsRF produced from WID processing of all 649 events recorded at ARU. (b) Stacked PdsRF produced using WID applied to ARU data from 2014 to 2015 (55 events), (c) Stacked PdsRF produced using WID applied to ARU data from 2012 (43 events). (d) Stack of traditional RF produced all 649 RFs computed from data recorded at ARU. To determine how well WID may work for temporary broad-band deployments, which typically last for 1–2 yr, we preformed additional tests of WID for 1 to 2 yr subsets of the ARU data. The most useful product of WID appears to be the PdsRFs, so we only focused on Pds phases in our discussion of tests using 1 and 2 yrs’ worth of data. Fig. 16(a) shows the PdsRFs using all 649 events recorded from ARU data. Figs 16(b) and (c) were produced using 2 yr of data (55 events) and 1 yr of data (43 events), respectively. Fig. 16(d) is the Pds stack of all 649 RFs from ARU using traditional RF processing. All the phase specific PdsRFs recovered Pds phases from the expected depths (Moho, Mid-lithosphere, 410 and 660). Phases from deeper layers, such as 410 and 660 km, have sufficient curvature so that they do not mis-stack as PXds phases on the raw RFs (Fig. 16b). There are more mis-stacked phases (or noise) throughout the transition zone depths on the traditional RFs than are present on all three of the PdsRFs computed using WID (Figs 16a–c). The positive Pds phase between 160 and 180 km on the PdsRFs computed from all the ARU data is not present on the PdsRFs computed from 2 yr of data and have a negative polarity on the PdsRFs produced from 1 yr worth of data. This phase is likely real since it is strong on the PdsRFs produced from all the data, but the differences on the short term RFs may reflect azimuthal variations due to differences in data distribution for the three different PdsRFs. There are phases that are unique to the PdsRFs from the 1 or 2 yr data sets that are not common to all three WID PdsRFs. The difference between these PdsRFs may reflect limitations in the distribution for the 1 and 2 yr data set that result in the inclusion of different miss-stacked reverberations or azimuthal variations in the local geology. However, the point of this paper is not to analyse ARU geology but to demonstrate that WID can produce RFs that are relatively free of higher order reverberations from shallow layers. It is clear that all the PdsRFs computed for different amounts of data have less noise and fewer mis-stacked reverberations than the traditional RFs, which demonstrates that WID can be effective for short-term deployments. To determine if there are substantial interpretable phases in the residual waveform from WID (Fig. 17), we stacked RFs made from the residual of the seismograms left after WID RFs (Fig. 18). We see there are no distinctive phases in the stacked residual data, except possibly a small residual of the Moho Pds in Fig. 18(a). This possible residual Moho amplitude reflects that the extracted amplitude from iterative process may not completely remove all amplitude from an interface. There are no obvious phases remaining in the stacks of all three residual that stand out significantly above the noise level. Figure 17. View largeDownload slide Residual receiver function for station ARU after WID. Figure 17. View largeDownload slide Residual receiver function for station ARU after WID. Figure 18. View largeDownload slide Stacks of residual receiver function (a) using Pds moveout corrections, (b) PPds moveout corrections and (c) PSds moveout corrections. Figure 18. View largeDownload slide Stacks of residual receiver function (a) using Pds moveout corrections, (b) PPds moveout corrections and (c) PSds moveout corrections. 3 CONCLUSION WID produces phase specific RFs that are easier to interpret than the stacks of traditional RFs. The WID method is the most useful for shallow phases where the pulse width is close to or larger than the amount of moveout. To get an optimal result from WID, RFs from a wide range of distances are required. PdsRFs are the most valuable of the phase specific RFs because the PdsRFs tend to have very little contamination from higher order reverberation compared to the PSdsRFs and especially the PPdsRFs. WID also works best using data with high signal-to-noise ratios. We have tested WID at many stations and found that it did not work well for many of the South American stations. However, it produced clear 3-D stacked images using California data from permanent and temporary deployments with Moho and LAB phases (Ainiwaer & Gurrola, in preparation) that are similar to recent published Sp work (Lekic et al. 2011). We also found that WID works well using transportable array data from Texas (Gurrola et al., in preparation). The Pds stacked RFs from WID (PdsRFs) are the most improved, owing to the distinct shape of the Pds moveout curve as compared to the moveout curves for the PPds and PSds phases and reverberations. The H–k analysis in WID is only used to identify a phase, but arrival times and amplitudes for specific phases on the RFs are estimated using cross-correlation of each individual seismogram, independently. Therefore, the resulting RFs can still be used to model 2-D or 3-D Earth structure. Acknowledgements We appreciate the comments from fellow students to improve the manuscript. We are grateful the IRIS data centre for providing seismological data. We thank the reviewers (G. Pavlis and A. Frederiksen) for comments and suggestions that greatly improved this manuscript. REFERENCES Aki K., Richards P.G., 2002. Quantitative Seismology . 2nd edn, pp. 144– 145, University Science Books. Andrews J., Deuss A. 2008. Detailed nature of the 660 km region of the mantle from global receiver function data, J. geophys. Res.  113( B6), B06304, doi:10.1029/2007JB005111. https://doi.org/10.1029/2007JB005111 Google Scholar CrossRef Search ADS   Ammon C.J., 1991. The isolation of receiver effects from teleseismic P waveforms, Bull. seism. Soc. Am.  81 2504– 2510. 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. https://doi.org/10.1029/JB095iB10p15303 Google Scholar CrossRef Search ADS   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( B2), doi:10.1029/2011JB008560. https://doi.org/10.1029/2011JB008560 Dueker K.G., Sheehan A.F. 1998. Mantle discontinuity structure beneath the colorado rocky mountains and high plains, J. geophys. Res.  103( B4), 7153– 7169. https://doi.org/10.1029/97JB03509 Google Scholar CrossRef Search ADS   Dugda M.T. 2005. Crustal structure in Ethiopia and Kenya from receiver function analysis: implications for rift development in eastern Africa, J. geophys. Res.  110( B1), B01303, doi:10.1029/2004JB003065. https://doi.org/10.1029/2004JB003065 Google Scholar CrossRef Search ADS   Dziewonski A.M., Anderson D.L., 1981. Preliminary reference Earth model, Phys. Earth planet. Inter.  25( 4), 297– 356. https://doi.org/10.1016/0031-9201(81)90046-7 Google Scholar CrossRef Search ADS   Gao S.S., Silver P.G., Liu K.H. & the Kaapvaal Seismic Group, 2002. Mantle discontinuities beneath Southern Africa, Geophys. Res. Lett.  29( 10), doi:10.1029/2001GL013834. https://doi.org/10.1029/2001GL013834 Gilbert H.J., Sheehan A.F., Dueker K.G., Molnar P. 2003. Receiver functions in the western United States, with implications for upper mantle structure and dynamics, J. geophys. Res.  108( B5), 2229, doi:10.1029/2001JB001194. https://doi.org/10.1029/2001JB001194 Google Scholar CrossRef Search ADS   Gurrola H., Simmons N.A. 2000. Multiple seismic discontinuities near the base of the transition zone in the Earth's mantle, Nature  405( 6786), 559– 562. https://doi.org/10.1038/35014589 Google Scholar CrossRef Search ADS PubMed  Gurrola H., Minster J.B., Owens T. 1994. The use of velocity spectrum for stacking receiver functions and imaging upper mantle discontinuities, Geophys. J. Int.  117( 2), 427– 440. https://doi.org/10.1111/j.1365-246X.1994.tb03942.x Google Scholar CrossRef Search ADS   Gurrola H., Baker G.E., Minster J.B. 1995. Simultaneous time-domain deconvolution with application to the computation of receiver functions, Geophys. J. Int.  120( 3), 537– 543. https://doi.org/10.1111/j.1365-246X.1995.tb01837.x Google Scholar CrossRef Search ADS   Gurrola H., Bernard Minster J. 1998, Thickness estimates of the upper-mantle transition zone from bootstrapped velocity spectrum stacks of receiver functions, Geophys. J. Int.  133( 1), 31– 43. https://doi.org/10.1046/j.1365-246X.1998.1331470.x Google Scholar CrossRef Search ADS   Huang X., Gurrola H., Parker J., 2016. A comparative study of popular methods to compute receiver functions applied to synthetic data, Seism. Res. Lett. , 129–121–129–4. 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. https://doi.org/10.1046/j.1365-246x.2000.00217.x Google Scholar CrossRef Search ADS   Langston C.A., 1977. The effect of planar dipping structure on source and receiver responses for constant ray parameter, Bull. seism. Soc. Am.  67 1029– 1050. Langston C.A. 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves, J. geophys. Res.  84( B9), 4749– 4762. https://doi.org/10.1029/JB084iB09p04749 Google Scholar CrossRef Search ADS   Lawrance J.F., Shearer P.M. 2006. A global study of transition zone thickness using receiver functions, J. geophys. Res.  111( B6), n/a. https://doi.org/10.1029/2005JB003973 Lekic V., French S.W., Fischer K.M. 2011. Science , 334 783– 787. https://doi.org/10.1126/science.1208898 CrossRef Search ADS PubMed  Ligorria J.P., Ammon C.J., 1999. Iterative deconvolution and receiver-function estimation, Bull. seism. Soc. Am.  89 1395– 1400. Liu K.H., Gao S.S., Silver P.G., Zhang Y. 2003. Mantle layering across central South America, J. geophys. Res.  108( B11), 2510, doi:10.1029/2002JB002208, B11. https://doi.org/10.1029/2002JB002208 Google Scholar CrossRef Search ADS   Nair S.K., Gao S.S., Liu K.H., Silver P.G., 2006. Southern African crustal evolution and composition: constraints from receiver function studies, J. geophys. Res.  111( B2), doi:10.1029/2005JB003802. https://doi.org/10.1029/2005JB003802 Owens T.J., Crosson R.S., Hendrickson M.A., 1988. Constraints on the subduction geometry beneath western Washington from broadband teleseismic waveform modeling, Bull. seism. Soc. Am.  78 1319– 1334. Plešinger A., Hellweg M., Seidl D., 1986. Interactive high-resolution polarization analysis of broadband seismograms, J. Geophys.  59 129– 139. Tauzin B., Debayle E., Wittlinger G. 2008. The mantle transition zone as seen by global Pds phases: no clear evidence for a thin transition zone beneath hotspots, J. geophys. Res.  113( B8), B08309, doi:10.1029/2007JB005364. https://doi.org/10.1029/2007JB005364 Google Scholar CrossRef Search ADS   Vinnik L.P. 1977. Detection of waves converted from P to SV in the mantle, Phys. Earth planet. Inter.  15( 1), 39– 45. https://doi.org/10.1016/0031-9201(77)90008-5 Google Scholar CrossRef Search ADS   Wang Y., Pavlis G.L. 2016. Generalized iterative deconvolution for receiver function estimation. Geophys. J. Int.  204( 2), 1086– 1099. https://doi.org/10.1093/gji/ggv503 Google Scholar CrossRef Search ADS   Zhu L., Kanamori H. 2000. Moho depth variation in southern California from teleseismic receiver functions, J. geophys. Res.  105( B2), 2969– 2980. https://doi.org/10.1029/1999JB900322 Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### DeepDyve Freelancer ### DeepDyve Pro Price FREE$49/month

\$360/year
Save searches from