Converted-wave reverse time migration imaging in subduction zone settings

SUMMARY

We use a newly developed 2-D elastic reverse time migration (RTM) imaging algorithm based on the Helmholtz decomposition to test approaches for imaging the descending slab in subduction zone regions using local earthquake sources. Our elastic RTM method is designed to reconstruct incident and scattered wavefields at depth, isolate constituent P- and S-wave components via Helmholtz decomposition, and evaluate normalized imaging functions that leverage dominant P and S signals. This method allows us to target particular converted-wave scattering geometries, for example incident S to scattered P, which may be expected to have dominant signals in any given data set. The method is intended to be applied to dense seismic array observations that adequately capture both incident and converted wavefields. We draw a direct connection between our imaging functions and the first-order contrasts in shear wave material properties across seismic discontinuities. Through tests on synthetic data using either S → P or P → S conversions, we find that our technique can successfully recover the structure of a subducting slab using data from a dense wide-angle array of surface stations. We also calculate images with a small-aperture array to test the impact of array geometry on image resolution and interpretability. Our results show that our imaging technique is capable of imaging multiple seismic discontinuities at depth, even with a small number of earthquakes, but that limitations arise when a small aperture array is considered. In this case, the presence of artefacts makes it more difficult to determine the location of seismic discontinuities.

1 INTRODUCTION

The shallow portion of subducting slabs generate Earth’s largest earthquakes, primarily through interplate earthquakes between the top of the slab and the overriding crust (Astiz et al. 1988). Seismic and tsunami hazard for near-coastal population centres depend upon the downdip limit of the locked seismogenic zone, the updip limit of near-trench rupture and placement of asperities (Petersen 2002; Hyndman 2013; Audet & Schaeffer 2018; Walton et al. 2021). Thus, accurate forecasts of hazard for near-coastal cities require knowledge of the detailed material properties of the entire megathrust. High-resolution 3-D images of Earth structure have been identified as a necessary component for improved risk management in the Cascadia subduction zone in particular (Gomberg et al. 2017).

The ideal resolution scale for imaging subduction zones is on the order of a few hundred meters to 1 km, given that this is the scale of the subduction channel (that is, the thin area between the subducting slab and overriding plate; e.g.  Nedimovic et al. 2003). However, existing techniques are not able to image structure a few hundred meters in wavelength across the entire length of the subducting slab. Active source seismic studies have the capability to illuminate material properties at high resolution, and are commonly used offshore to image the updip portion of the megathrust (e.g. Preston et al. 2003; Li et al. 2015; Shiraishi et al. 2019). But beyond the offshore and near-shore zone that may be studied via active source imaging, the resolution of structure found by seismic studies tends to degrade. Further onshore, the slab must be imaged using passive source methods, and most passive source studies rely on teleseismic waves, whose frequency content is too low to resolve finely detailed structure. These passive source methods include inversion for the location and strength of scatterers based on the Kirchhoff–Helmholtz integral representation of the full seismic wavefield (Bostock et al. 2001) and forward modelling of converted and back-scattered waves to characterize a multilayered dipping structure (Nicholson et al. 2005). These studies based on teleseismic wave analysis generally have relatively low resolution, with an upper frequency limit of approximately 0.5 Hz (e.g. Nikulin et al. 2009; Audet & Schaeffer 2018). Occasionally the receiver function method can be pushed higher, up to 2 Hz, which has lead to higher resolution imaging of Vp/Vs ratios in the subducted crust and demonstrated relationships between fluid pressure and megathrust fault slip style (Kato et al. 2010). Generally, however, as a result of the relatively low resolution of receiver function studies, it can be difficult to merge relatively low-frequency onshore data with offshore images; for instance, the National Science Foundation Cascadia Initiative was unable to combine an onshore-offshore broad-band receiver function line into a clear image of variations in physical properties from the locked to transition zone to episodic tremor and slip (ETS) zones (Janiszewski & Abers 2015; Audet & Schaeffer 2018).

High-resolution information about seismic discontinuities may be recovered using converted waves—that is, seismic phases that convert from direct S to converted P or direct P to converted S at a seismic discontinuity. Converted wave analysis has been used to image subducted slabs in many regions worldwide (e.g. Eberhart-Phillips & Reyners 1999; Horleston & Helffrich 2012; Gong & McGuire 2021). Pioneering converted wave studies used forward-wave modelling to match the delay times of observed converted waves at major discontinuities such as the intraslab oceanic Moho and the boundary between the descending and the overriding crust (James & Snoke 1994).

Seismic structure illuminated by converted waves may be imaged using reverse time migration (RTM) techniques. Methods based on back-propagation of direct and converted waves include Kirchhoff–Helmholtz integral approaches to back-propagate reflected and converted waves to their time of origin (Wapenaar & Haimé 1990; Shearer et al. 1999), and direct cross-correlation of time-reversed direct and converted wavefields measured on Earth’s surface (e.g. Shang et al. 2012; Duan & Sava 2015; Shabelansky et al. 2017; Li & Li 2022). The latter applications generally use imaging functions that cross-correlate either scalar or vector fields representing the direct and converted waves, and they are successful in refocusing energy into specular reflectors, that is one or more surfaces dividing two volumes with distinct material properties, with additional refinements addressing the issue of polarity reversal in P → S or S → P conversions (Duan & Sava 2015). These methods are particularly powerful because they are able to make use of seismic data from local earthquakes, which represent a potential underutilized source of high-frequency information about the subducting slab. There is an additional advantage to methods that cross-correlate only back-propagated wavefields, without requiring a forward calculation: these methods are source-independent and require little information a priori (Shabelansky et al. 2014).

Here, we use a merging of methods for converted wave RTM imaging. This algorithm includes a new approach for receiver wavefield extrapolation to evaluate time-reversed direct and scattered wavefields at the depths of hypothetical reflectors, based on an approximate solution of the reverse-time representation theorem. It is conceptually similar to the approaches of Duan & Sava (2015); Shabelansky et al. (2017), and Du et al. (2019) for evaluating an imaging function designed to correlate incident (P or S) waves and converted (S or P) waves that have been reconstructed at their point of origin. As in these approaches, the new method uses Helmholtz decomposition of the back-propagated receiver wavefields in order to construct the imaging function. In contrast to these approaches, our method produces images with a high signal-to-noise ratio, and is therefore capable of producing interpretable results using data from even a single earthquake. We also normalize our imaging functions to remove all dependence on the seismic wavefields from the contributing earthquake sources, so that our images reveal both the locations of seismic discontinuities and the associated ratio of converted-wave to direct-wave transmission coefficients, which is related to the difference in seismic shear velocity and density across the boundary.

The imaging functions that we construct are controlled by the contrast in density, Δρ, and shear velocity, Δβ, across a seismic discontinuity. Therefore, our imaging technique is highly sensitive to shear wave structure and density, which are strongly linked to the presence of fluid and other mechanical properties. This is significant because the properties of subduction zones that control rupture extent and velocity likely scale with porosity, fluid pressure, and the metamorphic state of the subducted sediments (Kato et al. 2010), all of which can contribute to variations in shear velocity and density. Downdip rupture extent is likely controlled by the detailed mechanical properties of the sediments in the subduction channel (Liu 2013; den Hartog & Spiers 2014). Furthermore, rupture velocity is related to shear velocity. The rock-sediment and sediment-rock interfaces above and below the thrust interface are expected to result in large amplitudes in our imaging functions, and these are reflective of the mechanical state of the subducted sediments to some degree. Thus, our imaging function has potential to target the subduction channel described by, for example Nedimovic et al. (2003).

We illustrate the proposed converted wave RTM imaging method with tests on a synthetic velocity structure representative of a subducting slab and surrounding continental crust and mantle. We focus on relatively high-frequency waves (3 Hz) that could be generated by local seismic events and observed by a dense receiver array at the Earth’s surface. Application of our methodology to synthetic data sets demonstrates the power of the method to recover major seismic discontinuities in the downdip end of seismogenic zones, specifically in the 15–40 km depth range onshore where other methods have difficulty. We also test the limitations of our imaging method when the receiver array aperture is small and/or data from few seismic sources are available.

2 METHODS AND TOOLS

2.1 Theoretical treatment of imaging conditions

Seismic waves emanating from an earthquake will interact with discontinuities at depth, producing transmitted and converted waves. We focus here on the converted waves, which are characterized to first order as direct P to scattered S, denoted P → S, Ps or $\acute{P}\acute{S}$⁠; and direct S to scattered P conversions, denoted S → P, Sp or $\acute{S}\acute{P}$⁠. We assume that the scattered waves are generated at specular discontinuities. We further assume that both the direct waves and converted waves are observed at a receiver array at Earth’s surface.

As described in Appendix A, we characterize the strength of specular discontinuities through imaging functions designed to recover normalized $\acute{P}\acute{S}$ or $\acute{S}\acute{P}$ wave conversion coefficients. This is accomplished by: (1) approximating the direct and converted wavefields at depth through suitably defined receiver wavefields (i.e. from time reversal of the wavefield recorded by seismic receivers at Earth’s surface), (2) isolating P and S contributions to these receiver wavefields through Helmholtz decomposition, (3) computing a zero lag cross-correlation on the extrapolated direct and converted wavefields to evaluate single-source un-normalized imaging functions and (4) stacking and normalizing these imaging functions over available sources.

For any wavefield u(r, t), we extract the P and S components through Helmholtz decomposition as (eqs A12 and A13)

Let u^‡ be the time reversed receiver wavefield. (u^‡)_P and (u^‡)_S then represent the reconstructed P and S wavefields at depth. Following common convention, we will subsequently write the reconstructed converted wavefields as (δu^‡)_S or (δu^‡)_P in order to emphasize that they include scattered wavefields.

At every imaging point the dip of the discontinuity is assumed a priori. At a potential location of wave conversion ξ, the stacked imaging conditions for normalized $\acute{P}\acute{S}$ and $\acute{S}\acute{P}$ (eqs A27 and A28) are, respectively,

In these expressions, j indexes the source; p(ξ, j) is the ray parameter associated with the direct wavefield from source # j at ξ, calculated via eq. (A23); α and β are P and S velocities on a reference structure and $\hat{\mathbf{n}}$ and $\hat{\mathbf{t}}$ are local normal and tangent vectors to the hypothetical discontinuity at ξ. $\hat{\mathbf{n}}$ and $\hat{\mathbf{t}}$ are determined by the specified dip of the discontinuity via eqs (A14) and (A24). θ is the angle of the ray parameter with respect to the dipping discontinuity rather than a horizontal axis. An illustration of these parameters is shown in Fig. 1. In our implementation, we do lateral averaging of the denominator ${\int }_0^T{\left|{({\mathbf{u}}_j^{‡}(\xi ,t)\cdot \hat{\mathbf{n}})}_P\right|}^2\mathrm{d}t$ or ${\int }_0^T{\left|{({\mathbf{u}}_j^{‡}(\xi ,t)\cdot \hat{\mathbf{t}})}_S\right|}^2\mathrm{d}t$ to ensure numerical stability. This is achieved by calculating the value of the denominator at all points in the mesh and then computing the average of this value over all points in an X–Y plane for each point at depth (Z) in the mesh.

There are three primary characteristics that differentiate our technique from similar methods explored by others, for example Shabelansky et al. (2017) and Du et al. (2019). First, we do not backpropagate the full seismograms. Instead, we use only the direct wavefield and the scattered wavefield, (δu^‡)_S or (δu^‡)_P. To accomplish this, we window all seismograms in the receiver gather to remove phases other than those of interest and backpropagate the windowed seismograms. This removes many artefacts that would otherwise be produced by non-relevant phases that interfere with one another. For the case of incident P producing scattered S, we eliminate the direct S wave and the S-wave coda, which can dominate real datasets and often results from body-to-surface wave scattering at shallow depths which could be incorrectly mapped to deeper structure. For the case of incident S producing scattered P, we eliminate the direct P wave.

Secondly, we enhance the magnitude of the imaging function kernel by considering only the components of the P and S wavefields that are expected to have the largest amplitude in our cross-correlation. Previous studies have taken a dot product of the P and S wavefields as the imaging function kernel (Shabelansky et al. 2017), but this results in an imaging function that has relatively low amplitude because large components of P are multiplied by small components of S, and vice versa. Our imaging kernel, shown in eq. (A17), multiplies the component of S tangential to the discontinuity plane by the component of P normal to this plane. These are the components that are expected to be the largest, so their product has a higher magnitude than a dot product would.

Lastly, we use a stable sum of ray parameters p(ξ, j) to normalize our stacked imaging functions. This results in an imaging function that does not depend on the wavefield (specifically, the incidence angles of source wavefields with respect to the discontinuity plane), so we are able to deduce the relative amplitudes of contrast in seismic material parameters across discontinuities at various depths in the target structure from our imaging function.

2.2 Code implementation

We use the spectral element code SEIS2.5D (Pollitz & Langer 2023) designed to simulate seismic wave propagation in 2.5-D. It is based on the quasi-static deformation code of Pollitz (2014) but extended to seismic wave propagation by adding inertial terms and seismic wave attenuation to the governing equations. It has been previously used to characterize regional seismic wave propagation (Pollitz & Mooney 2015) and to explore synthetic RTM imaging (Pollitz 2019).

In a spherical coordinate system, a 2-D model domain is defined about a pole of symmetry, varying with colatitude and depth but not with azimuth. 3-D seismic wave propagation on this 2-D structure is simulated in the frequency domain by expanding the total displacement field into a series of independent azimuthal components, the amplitudes of which are determined by solving a series of inverse problems in the 2-D domain, one for each azimuthal order number m. To achieve 2-D seismic wave propagation in the present study, we retain only the azimuthally symmetric (m = 0) component of the solution. We use SEIS2.5D to generate synthetic source wavefields using moment tensor point sources, and to evaluate receiver wavefields using three-component forces at the receivers.

The forward simulations include the effect of attenuation. The strongest attenuation is assigned to the crust where Q_μ = 60 and Q_κ = 150, consistent with commonly used crustal attenuation values (e.g. Montagner & Kennett 1996). A lossless medium is used for the computation of the time-reversed wavefields ${\mathit{u}}^{‡}$ in order to avoid instability. This procedure results in a small underestimation of the amplitude of the reconstructed wavefield at depth, affecting back-propagated S more than P. Relative to a case in which attenuation is not included in the forward calculation, this will produce a modest reduction in imaging amplitude in the P → S case and an increase in imaging amplitude in the S → P case.

The Helmholtz decomposition defined by eqs (1) and (2) requires evaluating second spatial derivatives of seismic displacement throughout the model domain, and this requires greater spatial resolution of the grid than is needed to evaluate only the displacements. The 2d grid has elements of varying size that are on the order of 50–60 per cent of the S-wavelength at the dominant frequency of 3 Hz. With 6 GLL points per 1-D element, we find this to be sufficient to produce images mostly free of numerical noise with a high signal to noise ratio. If we were calculating displacements only, we could use a mesh spacing of 90–100 per cent of the S-wavelength, nearly double what we need for the Helmholtz decomposition. The greater spatial resolution required for the Helmholtz decomposition can result in large data sets that would be difficult to visualize. This problem is fairly minimal for our 2-D examples, but it could become a greater challenge in a 3-D implementation.

2.3 Workflow

First, synthetic seismograms at a receiver array are generated using a forward model computed on a true (i.e. target) Earth structure that contains seismic discontinuities. Next, the resulting seismic record sections must be processed. Real data, or synthetic data produced via a time-domain wave propagation code, must be processed with a low-pass filter to remove very high frequency information, since high-frequency signals tend to exacerbate the artefacts from the Helmholtz decomposition calculations unless the grid is made more dense. We are using a code that does calculations in the frequency domain, so only perform our initial calculations up to 3 Hz and filtering is not required. We have not added noise to the synthetic seismograms in this study. A smaller signal to noise ratio is expected to reduce the amplitude of imaging functions at discontinuities relative to surrounding artefacts. This is one reason we have chosen a 3 Hz filter; in real seismic data, lower frequencies may be contaminated with ocean noise (Agnew & Berger 1978).

As we have previously discussed, the data in the receiver gather must be windowed to isolate phases of interest before it is backpropagated. This is accomplished by determining the P-wave arrival time via an STA/LTA algorithm. Then, for the case of incident S converting to P, we edit out everything before 5 s after the P arrival time in order to remove the direct P wave. This allows isolation of both the back-propagated converted wavefield ${(\delta {\mathbf{u}}_j^{‡}(\xi ,t)\cdot \hat{\mathbf{n}})}_P$ and back-propagated direct wavefield ${({\mathbf{u}}_j^{‡}(\xi ,t)\cdot \hat{\mathbf{t}})}_S$ required to calculate I_{S → P} via eq. (4). For the case of incident P converting to S, we determine the S-wave arrival time by multiplying the P-wave arrival time by 1.73 and edit out the S wave and everything after it, including the S-wave coda. This allows isolation of both the back-propagated converted wavefield ${(\delta {\mathbf{u}}_j^{‡}(\xi ,t)\cdot \hat{\mathbf{n}})}_S$ and back-propagated direct wavefield ${({\mathbf{u}}_j^{‡}(\xi ,t)\cdot \hat{\mathbf{t}})}_P$ required to calculate I_{P → S} via eq. (3).

After the seismic data are processed, we backpropagate them into a domain with a smoothed version of the true model (Fig. 2). During the backpropagation run, we perform the Helmholtz decomposition at each time step in each location in the mesh and write out the result. The imaging function is then calculated using an assumed slab dip, the decomposed wavefields and the terms inside the sums of eqs (3) or (4). At this stage, the ray parameter is also computed via eq. (A23). These steps must be performed separately for each earthquake source in the stack. The imaging functions from the separate earthquake sources are then stacked to obtain the complete result of eqs (3) or (4).

3 EXAMPLES

We use the synthetic subduction zone velocity model shown in the left panel of Fig. 2 as the true model for all of the examples in this study. This model has a Vp/Vs ratio of 1.73 throughout the domain. Although our code does not require that the mesh honours boundaries between different layers of the velocity structure, it is beneficial to construct a mesh that does honour these boundaries because if it does not, the velocity structure will be resampled, and this results in a form of smoothing that leads to less-sharp converted phases. The smoothed model for time reversed wavefield calculations is constructed by using a Guassian filter to smooth the true model. The true model has regular spacing so that the smoothing will not result in distortion of the layers. Our smoothed model is shown in the right panel of Fig. 2. The mesh used for this model has 66 cells in the z-direction and 63 cells in the x-direction. In SEIS2.5D, E is aligned with the x-direction. N, which is aligned with the y-direction, is not relevant to us because our simulations are 2-D.

The earthquake sources used for all of our examples are shown as stars in Fig. 2. The earthquakes are all located at a depth of approximately 40 km because we want a spacing of at least a few km between the source and the lowermost discontinuity so that the converted wave front will show up clearly on the record section. There is a slight variation in the depth of each earthquake to avoid potential additive source effects in the stacked imaging functions that might result in amplified source artefacts. Earthquakes have a shear moment of 1 × 10¹⁴ N m.

3.1 Isotropic source

The first example that we present is a simplified case intended to illustrate how the converted-wave imaging method works. We use 8 isotropic seismic sources located about 2 km apart at approximately 40 km depth. This example uses the full array of 129 stations.

We first do a forward run for each of the sources to produce synthetic seismic data. To illustrate the results of the Helmholtz decomposition, we have plotted snapshots of the decomposed forward wavefield from one of the sources in Fig. 3. A full movie of this wavefield simulation can be viewed in the Supporting Information. Note that we generally perform Helmholtz decomposition of only the time-reversed receiver wavefield, not the forward wavefield. This example shows how the decomposition cleanly separates the P and S wavefields. The primary wave front visible in the P wavefield, in the left-hand panel of Fig. 3, is the direct P wave which is the only direct wave produced by an isotropic source. The right-hand panel of Fig. 3 shows the S wavefield, which contains three primary wave fronts, each corresponding to a P → S conversion at one of the discontinuities. Note that we do not see clearly separated wave fronts for the two uppermost discontinuities because they are very close together; there is only 1–2 km separation between these boundaries. We would be able to distinguish between these discontinuities if we used a higher frequency, and hence smaller wavelength, in our calculations, but this would increase the computational cost.

Fig. 4 shows record sections for this source. The seismograms are consistent with the decomposed wavefield shown in Fig. 3. The Z-component, shown in the right panel, has a single strong arrival—the direct P wave—along with low-amplitude arrivals later in the record. These low-amplitude arrivals are due to scattering from Ps phases converting back to P waves at later (shallower) discontinuities. The E-component, shown in the left-hand panel of Fig. 4, has three strong arrivals after the P wave, corresponding to the three converted wave fronts seen in Fig. 3. We window out all data after the S-wave arrival time prior to backpropagation to reduce interference from the scattered waves that appear towards the end of the record section.

The stacked imaging function calculated with all eight seismic sources is shown in Fig. 5. All of the significant discontinuities have been recovered with a high degree of resolution. As we expect, the second boundary from the top has not been clearly recovered. This is consistent with our observation that we do not obtain a clean separate converted wave front from this boundary because it is so close to the uppermost discontinuity relative to the wavelength implied by our 3 Hz calculation.

Eq. (A25) indicates that the sign of the imaging function primarily corresponds to the sign of Δβ, defined as β_shallow − β_deep. In this case, we have a P → S imaging function, so we expect that when a discontinuity separates an upper higher-velocity structure from a deeper low-velocity structure, the sign of the imaging function should be negative. When a discontinuity separates an upper lower-velocity structure from a deeper higher-velocity structure, we expect the imaging function to have a positive sign. These expectations are consistent with our findings for all significant boundaries. The uppermost topographic discontinuity separates an upper higher-velocity structure from a deeper lower-velocity structure, and it has a negative imaging function with positive side lobes. The two deeper discontinuities separate upper lower-velocity structure from deeper higher-velocity structure, and they have positive imaging functions with negative side lobes. The opposite sign side lobes are caused by the narrow band nature of the time-reversed waveforms. These artefacts can be greatly reduced through the use of a higher frequency band, but this would result in the need for a denser mesh and hence significantly higher computational costs. We are able to recover the shallowest boundaries most strongly, and to a greater spatial extent, but all significant boundaries have been imaged well in this case.

The interpretation of the imaging function would not change if a different source time function was used. The frequency content of the source time function will affect the location and amplitudes of the sidelobes, but the amplitude and location of the central peak located on the discontinuity is expected to remain the same. This is because we will be filtering below the corner frequency (about 5 Hz) for the magnitude 2–4 earthquakes that we would use.

3.2 Non-isotropic source: wide-angle dense array

We now investigate both P → S and S → P imaging functions with moderately dipping non-isotropic double couple earthquake sources in our subduction zone structure. Because a non-isotropic source produces both direct S and direct P waves, there will be a greater quantity of scattered waves, since both direct S and direct P will produce converted, transmitted and reflected waves at each boundary. Therefore, windowing of the seismograms is even more critical in these cases.

Each of these examples uses the full array of 129 stations across the top of the domain. We choose our sources based on the type of converted waves we are targeting. Because our simulations are 2-D, all sources have a strike of 0° or 180°. For the P → S case, we use eight 45° dipping sources at the locations of the blue stars in Fig. 2(left-hand panel). This type of earthquake is expected to send significant direct P-wave energy upwards, thus producing strong converted Ps phases at the discontinuities above the source. For the S → P case, we use six +75° dipping sources at the locations of the red stars and six −75° dipping sources at the locations of the yellow stars. Additional sources were used in the S → P case because we do not have as much energy directed upwards; the optimal sources for S → P conversions have much of their energy directed towards the edges of the domain, so we used extra sources to enlarge the area that we are sensitive to and thereby enlarge the area of the discontinuity that we recover. If real data were used, the researcher should choose whether to compute the P → S or S → P imaging function based on their knowledge of the seismic source and the type of converted waves that may be observed.

The P → S and S → P stacked imaging functions are shown in the left- and right-hand panels, respectively, of Fig. 6. The P → S imaging function looks similar to the one we produced with the isotropic source, shown above in Fig. 5. But the S → P image both overlaps and recovers different portions of the discontinuities than the P → S image. This arises because Ps waves tend to refract to lower incidence angles (with respect to the dipping discontinuities) as P converts to S, making them more likely to be recorded by the receiver array, while Sp waves tend to have ray paths that diverge away from the far ends of the receiver array as S converts to P. Hence converted waves sample different portions of the discontinuities in the two cases. The lowermost discontinuity is well imaged in the S → P image, but it is not well imaged in the P → S image, where it is poorly resolved and slightly mislocated. This is because there are very few conversion points on the lower boundary due to geometry of P to S scattering at that depth. Both imaging functions have some slight artefacts around the sources and near the bottom boundary. The source artefacts, however, are greatly reduced relative to their appearance in a non-stacked single earthquake imaging function (see Fig. 9). Source artefacts are more prominent in S → P images because it is difficult to completely window out the direct P wave without also removing the Sp phase from the lowermost discontinuity.

While the imaging functions from the P → S and S → P cases are similar in magnitude, they are not identical. Eqs (A26) and (A29) indicate that the imaging conditions should differ only by sign. The magnitude difference arises due to our treatment of attenuation in the simulations. We do include attenuation in the forward runs used to produce synthetic seismograms because attenuation is present in the real Earth and we would like our seismograms to be as realistic as possible. However, we do not implement attenuation in the time-reversed simulations because time-reversing attenuation leads to instability. Attenuation is most pronounced in the continental crust, which is on average h = 23 km thick and has Qs = 60 and Qp = 90, based on the assigned Q_μ = 60, Q_κ = 150 and the relations given by eqs (9.59)–(9.60) of Dahlen & Tromp (1998). The attenuation of P and S waves through the crustal layer are $e^{(-Q{p}^{-1}\pi hf/\alpha )}$ and $e^{(-Q{s}^{-1}\pi hf/\beta )}$⁠, respectively, where f = 3 Hz is the dominant frequency. Inspection of eqs (A27) and (A28) then implies that I_{P → S} and I_{S → P} will be too low or too high, respectively, by a factor of e^{( − πhf*[1/(α Qp) − 1/(β Qs)])}. With α = 5.8 km s⁻¹ and β = 3.35 km s⁻¹, this factor is 1.9. Thus, our S → P imaging functions are expected to be biased low by a factor of 1.9 and P → S imaging functions systematically high by the same factor. This is consistent with the differing magnitudes of the imaging functions shown in Fig. 6.

To investigate the S → P imaging case further, we plot snapshots of the decomposed seismic wavefield from the forward run in Fig. 7. A full movie of this wavefield simulation can be viewed in the Supporting Information. We see that the direct P wave, the direct S wave, and the converted Sp waves all have significant amplitude. The direct P wave does produce converted S waves as well, but because the P-wave energy is lower for this source, the Ps conversions have a lower amplitude and do not show up well in these plots. Both wavefields have artefacts around the lower boundary, possibly as a result of calculating the Helmholtz decomposition in the absorbing PML layers; this is the likely cause of the lower boundary artefact seen in the imaging functions.

Fig. 8 shows the record sections for this earthquake source. The seismograms are consistent with the decomposed wavefields shown in Fig. 7; they show strong direct P and S waves, and coherent Sp arrivals on the Z-component. Between the direct P and direct S waves on the E-component, we see a very large number of scattered phases. This is due to the fact that both the direct P wave and the direct S wave are producing converted phases at each of the boundaries, and although the Ps conversions have lower amplitude, as we have noted, they do show up on the record section. We window out the direct P wave to reduce interference, but it would not be possible to disentangle the Sp phases of interest from the other scattered waves. However, because the Ps phases have low amplitude, they produce only a small amount of noise in the S → P imaging function shown in Fig. 6(right-hand panel). When real data are used, it would be necessary for the researcher to determine whether a P → S or S → P imaging function should be calculated, based on the amplitudes of the Sp and Ps arrivals in the data. The preceding examples would suggest, however, that since shear sources produce relatively strong direct S, Sp waves are generally expected to dominate over Ps.

3.3 Non-isotropic source: sparse event coverage

The examples that we have shown thus far use a moderate number of earthquakes in each stack: eight sources for the P → S image, and 12 sources for the S → P image. While this number of events could be observed in many subduction zones within a reasonable time frame, the Cascadia subduction zone has exceptionally low seismicity, so it may take a few years to observe that many local earthquakes with a seismic array in Cascadia. Therefore, we would like to investigate whether a single earthquake would be sufficient to recover meaningful high-resolution information about a subduction zone.

We select a single seismic source from those that contributed to each of the stacked images shown in Fig. 6, and plot the imaging function for that single event without normalization. We avoid normalization in this case because when fewer data are used, it is likely that the ray parameter will have values close to zero, making normalization unstable. We chose the seismic source located at approximately (−20, −40) km for both images. The results are shown in Fig. 9. Imaging functions for all seismic sources contributing to the stacked images shown in Figs 5 and 6 may be found in the Supporting Information.

Both single-source imaging functions recover at least part of each of the significant discontinuities. In the P → S image, shown in the left panel of Fig. 9, nearly the entire topographic region of the uppermost boundary is recovered, while the lower boundaries are only partially recovered, with the horizontal range of the imaged area decreasing with depth. In the S → P image, shown in the right-hand panel of Fig. 9, we see a similar phenomenon, although a smaller portion of the uppermost discontinuity is recovered. This is consistent with our observations for the stacked S → P image using the full array (Fig. 6, right-hand panel): Sp phases tend to diffract away from the edges of our array, so fewer converted phases are observed and the boundaries are not recovered to the same extent as in the P → S case.

Both the P → S and S → P single source images have less imaging power and stronger artefacts than the corresponding images with the wide-aperture array (Fig. 6), similar to the results of Du et al. (2019). Nonetheless, these examples suggest that our converted-wave imaging technique can recover seismic discontinuities to considerable depth using data from only a single suitably located earthquake using Ps or Sp converted waves.

3.4 Small-aperture array

The images shown in Fig. 6 were produced from backpropagated seismograms from 129 stations spread over about 35 km. Although this array geometry is not unrealistic, in practice nodal seismic arrays often cover a smaller area due to permitting issues, vegetation and other challenges. Therefore, we would like to determine how well imaging functions may be recovered using data from a smaller aperture array.

If only a small number of stations are used, we do not have sensitivity to the entire domain because the ray paths from our seismic sources to the limited number of receivers only covers a small area. For a single source and receiver, the recorded wavefield is generally sensitive to only a small volume surrounding the direct ray path based on finite-frequency amplitude and traveltime kernels (e.g. Marquering et al. 1999). Therefore, the imaging function calculation would not have any meaning outside the vicinity of the direct ray paths from our sources to the stations in the mini array. To reduce noise in the image and ensure that we do not mistakenly interpret areas that are not meaningful, we calculate the direct ray path from each source to the centre of our mini array and mask out the imaging function for that source in the area outside a zone near the direct ray path, roughly one Fresnel zone of the highest contributing frequency (∼3 Hz).

In our first example with a single mini array, we selected the 16 stations out of the original 129 that lie between approximately −16 and −20 km. All eight 45° dipping seismic sources shown as blue stars in Fig. 2 were used to form the stacked imaging functions. For each source, we backpropagated the windowed seismograms from the stations in the mini array, and performed the Helmholtz decomposition. Although calculating the imaging function for a given source, we also computed the direct ray path between the source and the centre of the array, and only saved the imaging function inside the area near the direct ray path. This results in the stacked imaging function in the left panel of Fig. 10. It recovers a similar pattern as the wide-angle array (Fig. 6) on all of the seismic discontinuities but over a limited area defined by the masking. As in the wide-angle array image, we expect the shallow topographic discontinuity, which separates an upper higher-velocity region from a deeper lower-velocity region, to be marked by a positive (blue) imaging function with negative (red) side lobes. Similarly, we expect the two deeper discontinuities to be marked by negative imaging functions with positive side lobes, consistent with observations. We do not normalize the imaging function in this case as well because normalization with the ray parameter is not stable due to the small amount of data used here.

Our second example has an additional mini-array consisting of 16 stations between −26 and −30 km. Again, all eight 45° dipping sources are used. For each source, we backpropagate the windowed seismograms from all 32 stations together and compute the Helmholtz decomposition. Then the imaging functions are calculated, and a direct ray path is computed from the source to the centre of each of the mini arrays. We save the imaging function from each source only in the vicinity of these direct ray paths. The stacked imaging function is shown in Fig. 10(right-hand panel). In addition to the areas that were imaged by the single mini array, we have successfully imaged a broader swath of the shallow topographic discontinuity as a positive region with negative side lobes, though the imaging function is distorted and more difficult to interpret. The deepest discontinuity has been imaged successfully in areas beneath each of the mini arrays as a bright negative spot with positive side lobes, but in areas not directly beneath the mini arrays, artefacts dominate the image. The middle discontinuity is perhaps imaged best in this case as a negative region with positive side lobes. Although the two mini arrays allow us to image a broader area than the single mini array does, we also have more noise in this image, so there is a significant trade-off. It may be better in practice to backpropagate data from mini arrays separately in order to aid in interpretability.

The masked area is significantly wider for the second mini array than it is for the first. This is due to ray paths bending around the thickest part of the low-velocity layer. For the first mini array, the ray paths are pinched as they tend to avoid the topographic low-velocity zone. For the second mini array, the ray paths travel through the low-velocity zone with less distortion and the resulting ray paths are distributed more widely.

4 VELOCITY MODEL ERRORS

Errors in the assumed velocity structure are likely to have an impact on calculated imaging functions. We explore this concern with an inaccurate velocity model for time-reversed wavefield calculations. The velocity structure shown in Fig. 11 is highly smoothed and has had the topographic layer on the surface of the slab removed. With this model, we investigate whether our technique can recover discontinuities at high resolution when available velocity models have low resolution, and whether we can recover the structure of a slab correctly when the available model contains inaccurate structural information.

The right panel of Fig. 11 shows a stacked S → P image recovered by back-propagating wavefields in this modified velocity structure. The workflow followed for this image is otherwise identical to workflow for the S → P image shown in Fig. 6. The shape of the topographic surface is recovered very well by this image, even though it is not present in the smoothed model. The boundaries recovered by this image are high resolution despite the high degree of smoothing in the velocity structure. However, directly beneath the topographic bulge on the slab, the lower boundaries are mislocated; the image infers that they are deeper than in the target model. This is because the back-propagated wavefields travel too quickly through the model due to the thinner low-velocity layer at the top of the slab.

5 DISCUSSION

In this study, we have presented examples of S → P and P → S converted wave imaging separately. For a given seismic source, a researcher must determine the dominant type of converted waves in order to choose which of these methods to use. If multiple sources of different types are used, the P → S and S → P imaging functions from the different sources may be combined by multiplying imaging functions of one type by −1 prior to stacking, in accordance with eqs (A26) and (A29). However, as we noted in Section 3.2, the lack of attenuation in time-reversed simulations also leads to a magnitude difference between S → P and P → S imaging functions, due to the fact that S waves are more strongly attenuated in the forward simulations (and in the real Earth). Therefore, one should also correct for this discrepancy prior to stacking P → S with S → P imaging functions. The correction factor that would account for propagation through an attenuating upper crust of thickness approximately 23 km in our examples, as presented in Section 3.2, is about 1.9. Our imaging functions are source independent, so stacking would not be affected by differences in the focal mechanism of earthquakes contributing to the stack. Focal mechanisms would matter in an imaging technique that requires forward simulations to generate the incident wavefield.

The imaging functions that we have calculated allow us to infer the relative strengths of seismic discontinuities in the structure. We can also interpret imaging functions in terms of the quantitative material contrast across a boundary. eq. (A25) states that imaging functions are proportional to Δρ and Δβ. In most cases, we would need to account for both attenuation and uneven sensitivity prior to making quantitative inferences about material contrasts from imaging functions. We can, however, make a quantitative comparison between imaged I_{P → S} and theoretical I_{P → S} for the case of isotropic sources, where the effect of uneven sensitivity is minimized. Table 1 presents this comparison at the locations of the three principal discontinuities at X = −17.2 km in Fig. 5, where the discontinuity depths are 21.80 km (continental crust–ocean sediment interface), 26.78 km (mid-oceanic crust discontinuity) and 33.41 km (oceanic Moho). Theoretical I_{P → S} compares reasonably well with the imaged quantity 1.9 × I_{P → S}. Note that in the real world, where the true model is unknown, one could often interpret images by neglecting Δρ and using eq. (A25) to estimate the contrast in S-wave speed, Δβ, across a given discontinuity, since the Δβ term is usually dominant.

For all of the examples in this study, we have used an ordinary Poisson’s ratio of 0.25. However, Poisson’s ratio can vary in subduction zones. This leads to large variations in Δβ, so imaging functions are likely to have large amplitudes due to this strong contrast in subduction zone regions.

While the imaging functions we have presented are relatively free of high-frequency artefacts, some low frequency artefacts remain in the images (e.g. Figs 5 and 6). Note that the imaging functions I_{P → S} and I_{S → P} are designed to recover the location and amplitude of specular reflectors, so they differ from kernels designed to image volumetric scatterers (Zhu et al. 2009; Luo et al. 2013). The use of the Helmholtz operators shown in eqs (1) and (2) acts as a high-pass temporal and spatial filter, helping to remove low frequency artefacts that are common in conventional reverse time migration with cross-correlation imaging conditions. The remaining low frequency artefacts could be addressed with least squares RTM (Dai et al. 2012). They could also be addressed through consideration of multiple frequency bands, as I_{P → S} and I_{S → P} are expected to be independent of frequency at the correct locations of seismic discontinuities but depend on frequency elsewhere.

When real data are used, windowing would be especially critical because of the likelihood that the imaging function could be contaminated by undesired interactions between other seismic phases. In particular, everything after the direct S wave should be windowed out to avoid added noise from surface waves. This also results in a less expensive backpropagation simulation, because the record sections can be cut after the direct S-wave arrival time. In some cases, body to surface wave scattering could also contaminate the interval between the P and S direct waves. However, in many datasets, converted phases are the largest arrivals between P and S so a strong signal-to-noise ratio is still likely. We also expect that most mini-arrays should be able to isolate the arrivals with near-vertical incidence, so the converted phases of interest could be identified and other phases removed prior to backpropagaton.

Seismic data should also be filtered with a low-pass filter to avoid the introduction of artefacts when the Helmholtz decomposition is performed. If one desires to resolve very small-scale structure, a higher cut-off frequency could be used, but the mesh must be more dense to avoid artefacts. This would result in a greater computational cost.

Our imaging method requires that the dip of the slab be known a priori.  Pollitz (2019) found that an incorrect assumed dip has only a modest effect on imaging function amplitudes, provided that the direction of the dip is known, which is a given in subduction zone settings; in fact, generally the large scale dip of the slab is well known to within a few degrees at most. At the wavelengths implied by a 3 Hz low-pass filter, this is what dominates the orientation of $\hat{\mathbf{n}}$⁠. Initial values would come from synthesis studies like Slab2.0 (Hayes et al. 2018), recorded seismicity, or active source studies if relevant.

The workflow that we have presented in this study is applicable to seismic sources that are located underneath the subducting slab. However, our method for converted-wave imaging can be extended to other wave types that involve interaction of a specific wave type with a seismic discontinuity and can be approached via wavefield reconstruction at depth. This could include imaging with reflected waves or twice reflected waves (Shiraishi & Watanabe 2022).

The imaging function derivation presented in Appendix A applies in 3-D as well as in 2-D. For a 3-D imaging calculation, all input quantities, such as seismograms, $\hat{\mathbf{n}}$ and $\hat{\mathbf{t}}$ should be specified in three dimensions. The workflow and calculations are the same. However, as we have previously mentioned, we expect that a 3-D implementation would be more challenging because the decomposed 3-D P and S wavefields require very large memory allocations. An efficient workflow would be required to visualize the decomposed wavefields and compute the 3-D imaging function.

6 CONCLUSION

In this study, we have presented a workflow for converted-wave elastic reverse time migration imaging that enables us to recover high-resolution images of subduction zone-like structure. The imaging functions that we use are sensitive to shear wave velocity structure. Our synthetic examples demonstrate that when data from a wide-angle, dense seismic array are available, this technique can recover broad swaths of multiple discontinuities at high resolution (e.g. Fig. 6). Multiple discontinuities at depth may be recovered even if data from only a single earthquake are used (Fig. 9), though a smaller horizontal range of the discontinuity is imaged. Similarly, if a small-aperture seismic array is present, this technique can recover a small segment of each significant discontinuity (Fig. 10), though significant artefacts are present which make interpretation of the images more difficult. Because our imaging method can recover seismic discontinuities using only a small number of seismic sources, we intend for it to be used to image regions that have previously been poorly imaged due to low seismicity and data sparsity, such as the Cascadia Subduction Zone.

SUPPORTING INFORMATION

suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

We thank Ettore Biondi, Elizabeth Cochran, Christina Morency and an anonymous reviewer for thoughtful comments which helped us improve the paper. Any use of trade, firm or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

DATA AVAILABILITY

No new data were generated or analysed in support of this research. The SEIS2.5D code is available from https://code.usgs.gov/esc/2.5D_Seismic_Wave_Code (Pollitz & Langer 2023).

References

APPENDIX A: CONVERTED-WAVE IMAGING FUNCTIONS

A1 Scattered and source wavefields

A1.1 Elastodynamic Green’s function

To use index notation we work in a x₁ − x₂ − x₃ Cartesian coordinate system. Following Aki & Richards (1980), we define the elastodynamic Green’s function $G_{in}^{(0)}(\mathbf{x},t,\xi ,\tau )$ as the i −th displacement component at position x and time t resulting from application of a unit point force at ξ in the ${\hat{\mathbf{x}}}_n$ direction at time τ on a reference model. Specifically, $G_{in}^{(0)}={\hat{\mathbf{x}}}_i\cdot \mathbf{u}(\mathbf{x},t,\xi ,\tau )$⁠, where displacement u satisfies

where ρ is density, $\mathit{\sigma }$ is the stress tensor, and δ denotes the Dirac delta function.

A2 Wavefield reconstruction at depth

The receiver wavefield is used to reconstruct both the scattered seismic wavefield and source wavefield at the locations of scattering in Earth’s interior. Receiver wavefield extrapolation is a common problem in seismic imaging (Masson et al. 2014; Yuan et al. 2017). In our approach, we isolate the scattered and source wavefields from the recorded wavefield, then extrapolate them back in time.

Assume a local x₁ − x₂ − x₃ Cartesian coordinate system, where x₃ points upward (Fig. A1). Suppose that recordings of the total wavefield are available at the ‘acquisition surface’ S^free, the free surface of the Earth at x₃ = 0, which we suppose to be a heterogeneous elastic half-space in x₃ ≤ 0. In addition to S^free, let S^int be an interior surface such that S^free∪S^int surrounds a large volume including the seismic source, and let Σ be a small surface surrounding the seismic source. The time-reversed representation theorem (e.g. eq. 2 of Masson et al. 2014) provides the solution for the time-reversed wavefield, which we denote with $\hat{\mathbf{u}}$⁠:

where c_ijkl are components of the elasticity tensor, n_j are components of the outward unit normal vector to the surface, and $\mathit{\sigma }$ is the stress tensor associated with u. T is an arbitrary end time of a seismogram. Note that eq. (A2) is applicable to any wavefield (u and $\mathit{\sigma }$⁠) that satisfies the elastodynamic wave equation with sources that are not inside the volume defined by S^free∪S^int∪Σ. We decompose the total wavefield u and $\mathit{\sigma }$ into their respective contributions from the incident (i.e. upgoing) wavefield and reflected (i.e. downgoing wavefield) as

Conceptually, we regard the incident wavefield as that which satisfies the elastodynamic wave equation with the empty half-space x₃ > 0 replaced with a homogeneous elastic half-space with material properties identical to those at x₃ = 0; hence eq. (A2) is applicable to u^inc and ${\mathit{\sigma }}^{\mathrm{i}\mathrm{n}\mathrm{c}}$ regardless of the boundary conditions imposed on G⁽⁰⁾. As we wish to reconstruct the incident (i.e. upgoing) wavefield, we thus apply eq. (A2) under the assumption that reflected waves are neglected, i.e., ${\mathbf{u}}^{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{l}}=\mathbf{0}$ and ${\mathit{\sigma }}^{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{l}}=\mathbf{0}$⁠. We suppose that both the interior surface and the source region are sufficiently distant from the region of interest that contributions of the initially radiated wavefield from the integrals over S^int and Σ may be neglected. For the remaining integral over S^free, the leading term of the integrand is zero because the Green’s functions are traction-free on S, hence

In Appendix B, we determine ${\mathit{\sigma }}_i^{\mathrm{i}\mathrm{n}\mathrm{c}}$ at x₃ = 0 for the P–SV system on a uniform isotropic elastic half-space, assuming small incidence angle and hence small ray parameter p. There it is shown that for either incident P or SV, with error of order (αp)² or (βp)², the incident-wavefield traction at a point on x₃ = 0 is given by

where ρ, α and β refer to the material properties at x₃ = 0 and v is the total velocity.

Eqs (A4) and (A5) yield the solution

Although eq. (A5) was derived for plane-wave conversion at the free surface on a uniform elastic half-space, the Green’s function G⁽⁰⁾ on the heterogeneous model may be used as long as the dominant wavelength of wave propagation is much shorter then the spatial scale of near-surface variations in material properties. This would extend to the free surface of a spherical model as long as its curvature is much larger than the dominant wavelength. We assume that eq. (A6) is applied with a suitable smooth heterogeneous reference model that is absent of discontinuities.

Let the seismic wavefield be observed at a collection of receivers {r_r|r = 1, ⋅⋅⋅, N} generated by a number of sources. Let u^obs denote the observed wavefield and v^obs = ∂_tu^obs. Since we have only a finite sample of the wavefield at Earth’s surface, we approximate eq. (A6) with a sum over the receivers:

If we make the substitutions τ → T − τ and t → T − t, then

It is convenient to define

An integration of eq. (A8) by parts yields

Note that $u_n^{‡}(\xi ,\tau )$ is distinct from the commonly used receiver wavefield defined by

The salient difference is in the weighting of horizontal and vertical force components, which are in the ratio α: β for $u_n^{‡}$ and 1: 1 for $u_n^{†}$

A3 Imaging functions for ratios of reflection/transmission coefficients

Let u⁰ denote the hypothetical wavefield on the reference model. Eqs (A9) and (A10) may be applied to both the observed wavefield u^obs and residual wavefield δu = u^obs − u⁰ in order to reconstruct the source and scattered wavefields at depth. Knowledge of both the reconstructed scattered wavefield and source wavefield permit a straightforward way to characterize seismic discontinuities through evaluation of their constructive interference at the locations of such discontinuities. Here we define imaging functions that describe the scattering of incident P to converted S or incident S to converted P that use wavefield-extrapolation estimates of both wavefields. They are designed to target contrasts in material properties across a hypothetic discontinuity at ξ, with $\hat{\mathbf{n}}(\xi )$ being the unit normal to the discontinuity surface at ξ (Fig. A1).

For incident body waves, which are furnished by the back-propagated waves contained in u^‡, we assume that u^‡ represents a ray that depends upon propagation distance x with dependence ∼(x − vt), where v is the phase velocity.

are unit vectors pointing in the direction of P-wave or S-wave propagation.

Let $\hat{\mathbf{t}}$ be the unit tangent vector to the hypothetical plane that is parallel to the projection of the incident wave polarization onto that plane,

The reconstructed wavefield u^‡ may be separated into P and S wavefields with Helmholtz decomposition, for example eq. (4.1) of Aki & Richards (1980):

and similarly for δu^‡. For the case of incident P producing scattered S we define the kernel

The kernel for the case of incident S producing scattered P is defined as

Finally, we define normalization factors

The functions ${\tilde{K}}_{P\to S}$ and ${\tilde{K}}_{S\to P}$ can be related to reflection or transmission coefficients that describe wave conversion across a specular reflector. We consider here the cases of incident upgoing P or S through a lower medium producing transmitted upgoing P or S through an upper medium. The wavefields u and δu are supposed to be provided by recordings from a receiver array on Earth’s surface. Direct inspection of the transmission coefficients for converted waves $\acute{P}\acute{S}$⁠, $\acute{P}\acute{P}$⁠, $\acute{S}\acute{P}$ and $\acute{S}\acute{S}$ given in table 5.4 of Aki & Richards (1980) then leads directly to the following relations:

where α and β are understood to be the density and P and S velocities on the reference model. Here p(ξ) = sin θ/c is the ray parameter defined with respect to the dipping discontinuity plane at ξ, where c is phase velocity (=α or β) and θ is the incidence angle of P or S with the boundary at an imaging point ξ (Fig. A1). One estimate of p is through weighted averages

In eq. (A23) we have dropped the argument (ξ, t) from (u^‡)_P, (u^‡)_S, ${\hat{\gamma }}_P$ and ${\hat{\gamma }}_S$⁠. As an alternative to numerical estimates of γ_P (or γ_S) and p, one may use an a priori time-independent estimate of γ_P or γ_S from knowledge of the source and path and use simply $p={\alpha }^{-1}{\hat{\gamma }}_P\cdot \hat{\mathbf{t}}$ or $p={\beta }^{-1}{\hat{\gamma }}_S\cdot \hat{\mathbf{t}}$ in place of eq. (A23). In any case, with the definition of $\hat{\mathbf{t}}$ in eq. (A14), we have p ≥ 0. Information on the geometry of incidence is embodied instead in $\hat{\mathbf{t}}$⁠, as is evident from eq. (A14) or the alternative expression

The transmission coefficients are provided in chapter 5 of Aki & Richards (1980). Eqs (A21) and (A22) are in principle applicable to imaging the coefficient ratios on non-planar and/or dipping interfaces, as long as p is referenced to the local dip of the interface. Because the coefficients $\acute{P}\acute{S}$ and $\acute{S}\acute{P}$ are proportional to p, the sign of their product with p does not change as incidence angle changes sign, a property shared by the proposed converted-wave imaging functions of Duan & Sava (2015) and Shabelansky et al. (2017). In the case of imaging contrasts in shear modulus Δμ or bulk modulus Δκ, the right-hand sides of eqs (A21) and (A22) are both positive if modulus values in the upper medium are less than those in the lower medium, and negative in the converse case. Hence the kernel ratios provide information on both the size and sign of hypothetical discontinuities at imaging points.

Inspection of the formulas for the transmission coefficients shows that with error of relative order (αp)² or (βp)², the quantities $p^{-1}\frac{\acute{P}\acute{S}}{\acute{P}\acute{P}}$ and $p^{-1}\frac{\acute{S}\acute{P}}{\acute{S}\acute{S}}$ depend upon the material property contrast and not upon p (Chapter 5 of Aki & Richards 1980). Specifically, if Δρ, Δβ and Δα are the material property contrasts across the discontinuity (the lower-medium value minus the upper-medium value), then with error of order (αp)² or (βp)²,

With information from several seismic sources available, we use eqs (A17)–20) and eqs (A21) and (A22) to obtain stacked estimates of these quantities, which are the imaging functions we wish to calculate:

where j indexes the source.

Note that eq. (A26) implies that

The converted-wave imaging functions given by eqs (9) and (10) of Shabelansky et al. (2017) which correspond to our ${\tilde{K}}_{\mathrm{\Delta }\mu }(\xi {)}_{P\to S}$ and ${\tilde{K}}_{\mathrm{\Delta }\mu }(\xi {)}_{S\to P}$ involve time integrals of the vector dot product [u^†(ξ, t)]_S · [u^†(ξ, t)]_P, where u^† is the conventional time-reversed receiver wavefield of eq. (A10). Note that eqs (A27) and (A28) are similar to their imaging functions with three chief differences. First, our source wavefield is provided by u^‡, that is eq. (A9), rather than u^†. Secondly, terms expected to be dominant in practice, such as ${[{\mathbf{u}}^{‡}\cdot \hat{\mathbf{t}}]}_S$ and ${[{\mathbf{u}}^{‡}\cdot \hat{\mathbf{n}}]}_P$ in eq. (A27), multiply each other rather than quantities such as ${[{\mathbf{u}}^{‡}\cdot \hat{\mathbf{n}}]}_S$ or ${[{\mathbf{u}}^{‡}\cdot \hat{\mathbf{t}}]}_P$ that are expected to be smaller (of relative order αp or βp). Thirdly, the manner of normalization, involving a sum of ray parameters over contributing sources, is different. In particular, since p(ξ, j) ≥ 0 owing to the definition of p in eq. (A14), the normalization factor 1/∑_jp(ξ, j) in eqs (A27) and (A28) is stable if available sources sample a diversity of incident angles.

The products in the integrand of eq. (A17) for the P → S kernel resembles a receiver function migrated to depth, with the term ${[{\mathbf{u}}^{‡}(\xi ,t)\cdot \hat{\mathbf{n}}]}_P$ playing the role of a vertical-component incident P and ${[{\mathbf{u}}^{‡}(\xi ,t)\cdot \hat{\mathbf{t}}]}_S$ the role of a radial-component converted S wave. eq. (A17) is a generalization that accounts for the dip of the hypothetical interface and which propagates the incident and converted wave components back to the depth of origin using the elastodynamic wave equation. The same interpretation applies to the S → P case in eq. (A19).

APPENDIX B: TRACTION ASSOCIATED WITH THE INCIDENT WAVEFIELD

Here we address the question of the traction associated with a plane wavefield incident on the free surface of uniform elastic half-space. We derive it for the case of P-wave incidence in the P–SV system. The task is, with only knowledge of the recorded total displacement (or velocity) wavefield, to separate the incident wavefield contribution from the reflected wavefield contributions to the traction at the free surface; while the sum of these components is zero, the contribution of each separate component is non-zero.

As we follow the approach of section 5.2.2 of Aki & Richards (1980) for characterizing the incident and reflected wavefields, we adopt their notation and conventions. Let x and y denote the horizontal coordinates and z the vertical coordinate measured positive downward. The density and P and S wave speeds are uniform half-space are denotes ρ, α and β, respectively. Let i be the angle of incident and reflected P waves, and j the angle of the reflected S waves. Incident and reflected wavefield potentials for plane P waves propagating in the x direction are given by

and the reflected S-wave potential by

The total displacement field is given by

The ray parameter is given by p = sin i/α = sin j/β. The incident P and reflected P and S velocity wavefields are given by

With error of order p², the incident and reflected traction is given by

The boundary conditions

yield, with error of order p²,

Since v = v^inc + v^refl, eqs (B5) through (B10) lead to

In the case of incident SV waves, we parametrize incident and reflected wavefield potentials for plane S waves propagating in the x direction by

and the reflected P-wave potential by

The total displacement field is given by

The incident S and reflected S and P velocity wavefields are given by

With error of order p², the incident and reflected traction is given by

The free-surface boundary conditions of eq. (B9) yield

Eqs (B16)–(B20) lead again to the relation of eq. (B11) between the traction at the free surface and the total velocity field.