Toward Waveform‐Based Characterization of Slab & Mantle Wedge (SAM) Earthquakes

Earthquakes in subduction zones occur in the slab mantle, in the subducting crust, on the subduction plate interface, and, in some cases, in the mantle wedge–regions that are separated by strong seismic discontinuities. These discontinuities are typically imaged with techniques using teleseismic waves, while local earthquakes are located based on arrival times. While this combination of imaging and earthquake location provides a good initial overview of where the earthquakes are located, the uncertainties associated with the two approaches are too large (i.e., few kilometers) to robustly identify on which side of a discontinuity (with thickness ∼ 100 m) the earthquakes occurred. Here we investigate how the waveforms of local earthquakes, which contain secondary phases arising from wave scattering at discontinuities, can be exploited to determine the source region of subduction zone earthquakes more robustly. Our investigation involves a three‐step approach and includes an application to data from western Greece. First, to identify characteristic secondary phases, we analyzed synthetic seismograms from a generic 2‐D subduction zone. Second, to enhance the visibility of secondary phases in field data, we implemented a workflow to process three‐component seismograms. Third, to identify individual secondary phases in the data, we matched their timing to arrivals computed in a 3‐D velocity model. We identified on average two to three secondary arrivals per station. These include P‐ and S‐reflections from the plate interface which indicate hypocenters in the mantle wedge, and P‐reflections from the slab Moho which indicate hypocenters on the plate interface and in the subducting crust.

. (a) Sketch of a subduction zone with SAM earthquakes and their potential near-vertical raypaths. Depicted rays do not include possible surface multiples. Only raypaths of first-order wave interactions are plotted (i.e., one conversion, one reflection, or one combined conversion & reflection at the same interface). Phases in bold are expected to be particularly well observed (see ray geometrical consideration in Section 2 and synthetics in Section 4). (b) Seismic discontinuity model of the Peloponnese-Attica region in the Western Hellenic subduction zone. The four discontinuities generate reflected and converted arrivals in the waveforms of SAM earthquakes (green dots). Interfaces are colored as follows; blue/green shading -surface topography; light brown -continental Moho; red -subduction plate interface; blue -slab Moho (see Figure S4 for a zoomed-out view of the model). Earthquakes within the system are typically located through methods that provide absolute earthquake locations. Locations can be refined for better relative accuracy through methods such as double-difference relocation (Waldhauser & Ellsworth, 2000). Such studies have located earthquakes at the top of the slab, close to or within the LVL, and sometimes in a second slab-parallel layer 15-40 km below the slab top. The presence of two layers of earthquakes is referred to as a double Wadati Benioff Zone. While earthquakes in the lower plane clearly occur within the slab mantle, it is not always clear whether earthquakes in the upper plane occur in the slab crust or in the slab mantle (e.g., Abers et al., 2013).
Dense temporary deployments are able to record how the waveforms of SAM earthquakes change due to their location relative to the discontinuities. However, this information is typically not used in the common approach for locating earthquakes in subduction zones. To solve the earthquake location problem relative to the discontinuities, we introduce a waveform-based method that draws on common preprocessing routines from local earthquake analysis, receiver functions and reflection seismology. We apply this methodology to a high quality data set recorded across the Peloponnese peninsula in Western Greece. This location in the Western Hellenic subduction zone (WHSZ) is an excellent natural laboratory since recent studies have suggested that earthquakes in this location occur in all four source regions. Analysis of the data set has yielded high-resolution images from receiver functions Suckale et al., 2009), a 3-D velocity model from local earthquake traveltime tomography (Halpaap et al., 2018), and a catalog of precise hypocenter locations from careful relocation (Halpaap et al., 2019). Here we focus on the exact location of earthquakes relative to the seismic discontinuities, in particular whether the new methodology can support the occurrence of earthquakes in the mantle wedge.

Secondary Arrivals in Local Recordings of SAM Earthquakes
Secondary arrivals of teleseismic earthquakes are regularly used to image subduction zone structure, while locally recorded SAM earthquakes are commonly used to study the shape of the Wadati-Benioff zone. The secondary arrivals of SAM earthquakes have been exploited far less. Most existing studies focus on a subset of the secondary arrivals to better constrain subduction zone structure, for example, using guided waves (i.e., delayed high-frequency waves, e.g., Abers, 2005;Martin & Rietbrock, 2006), or S-to-P conversions from the top and bottom of the LVL (e.g., Horleston & Helffrich, 2012). A handful of studies have explored the use of seismograms from local earthquakes, using interferometry (Poliannikov et al., 2012), autocorrelation methods (D. Kim et al., 2019), and processing approaches similar to RF analysis (Bock et al., 2000). Some authors have noted a difference in waveforms depending on the source region of earthquakes (e.g., Uchida et al., 2010), or have investigated waveform fits at periods above 1-6 s for SAM earthquakes with magnitudes 5.7 (e.g., Chen et al., 2007;Okamoto et al., 2018;Song et al., 2009;Wang et al., 2014). However, there has not yet been a systematic study on how secondary arrivals from SAM earthquakes can be used to characterize the earthquake's source region.
At teleseismic distances, incoming waves can be considered planar, with virtually constant incidence angle, amplitude and polarity recorded across the stations. The P-and S-waves are well-separated in time (4-10 min), thus their converted arrivals (delays of 1-20 s) are easily identified. The incoming wave always comes from below the discontinuities, although surface-reflected multiples both near the source and the receivers may cause secondary signals. Source-side effects are canceled out during deconvolution (Bostock, 2004), while receiver-side multiples are recognized by arrival time.
At local distances from the source, waves travel as curved wavefronts from a point source (for small earthquakes) to a typical seismic deployment, where the waves are recorded at a range of incidence angles. The three-dimensional radiation pattern of the earthquake source (assumed to be a double couple) results in differing signal amplitudes and polarities for P-and S-waves between each recording station. Direct P-and S-phases are not well separated in time: P-to-S conversions may arrive at the same time as S-to-P conversions for discontinuities such as the overriding Moho (see Figure 1a). The P-and S-coda may contain overlapping wave conversions from above the source, and reflections from below the source. Reflections are expected when an earthquake occurs above a discontinuity for near-vertical observations. Both the slab Moho and plate interface may give rise to reflected arrivals for sources in the mantle wedge, while only the slab Moho HALPAAP ET AL.
10.1029/2020JB021573 4 of 28 gives rise to reflected arrivals for sources in the slab crust. Such reflected arrivals vary considerably between stations due to variations in source radiation and reflection coefficients as a function of incidence angle.
With multiple secondary phases potentially arriving at the same time, the amplitude ratio between secondary and primary arrivals determines whether an arrival may be observed. This depends on the ratio of radiated P-and S-wave energy, the source's radiation pattern, and the seismic attenuation along the raypath. A typical earthquake radiates 9-25 times more S-wave than P-wave energy (e.g., Madariaga, 2015;Prieto et al., 2004). Not considering raypath effects, this translates to a 3-5 times larger amplitude for the S-than the P-wave, owing to the quadratic relationship of energy E and amplitude A, 2 E A  . With increasing distance, particularly for raypaths through hot parts of the mantle, this effect can be negated due to stronger attenuation of S-waves. For reflected and transmitted (converted) waves, the relative amplitudes are given by the Zoeppritz equations (Zoeppritz, 1919, simplified in Schoenberg & Protazio, 1992Shuey, 1985), or can be described with scattering patterns when using the approximation of Born scattering (e.g., Beylkin & Burridge, 1990;Millet et al., 2019). These patterns tell us that, given a discontinuity with an equally large perturbation in P-and S-velocity, the P-to-S transmission coefficient is the same as the S-to-P transmission coefficient. Considering no attenuation, the secondary arrivals of the S-wave hence exhibit a larger amplitude than their P-counterparts. Once we consider attenuation, S-to-P conversions are even more likely to be observed than their P-to-S counterparts. Secondary phases from discontinuities close to the source are our main focus, since we expect these phases to differ depending on source region. Conversion near the source leads to the wavefront traveling mostly as the type of wave (P or S) that is observed at the receiver. This implies that S-to-P converted phases start out with a larger amplitude at the source and keep a larger portion of the initial amplitude than P-to-S conversions. Therefore, it should be considerably easier to observe S-to-P conversions than P-to-S conversions. Still, the individual geometry of the source and raypaths between source and receiver will be a dominant factor for the observability of secondary phases.
Earlier studies have used wave conversions from local earthquakes to look at structure (e.g., Horleston and Helffrich, 2012;Matsuzawa et al., 1990;Regnier et al., 1994;Zhang et al., 2004) but these generally focused either on wave conversions from slab earthquakes below the target discontinuity, or Moho reflections of shallow earthquakes and largely avoided cases in which the path of secondary arrivals is ambiguous.
Overall, earthquakes in the slab mantle offer the least ambiguity with regard to secondary-arrival raypaths (all raypaths with first-order interactions, i.e., only conversions, are upgoing), which makes them well-suited for mapping studies. However, not all subduction zones exhibit prominent seismicity below the slab Moho, the WHSZ being one such example (Halpaap et al., 2018). There is another class of earthquakes that offers a reduced set of multiple raypaths: Earthquakes that occur on the plate interface (Uchida et al., 2010), whose location precludes conversions or reflections from the interface itself. For that reason, interface earthquakes will prove valuable in this study.
In this study, we explore how the source region of SAM earthquakes can be determined from seismograms recorded at short distances from the source. We present two methods, one to calculate synthetic seismograms (spectral element method (SEM)), and the other to compute multiphase arrival times (fast marching method) and describe the respective models. Synthetic seismograms are used to explore how secondary arrivals can be identified based on specific criteria. We then use observed seismograms from a data set of well-located earthquakes in the WHSZ, applying a processing scheme to enhance the signal of secondary arrivals. Our discussion focuses on three questions: (a) What secondary arrivals are best suited to confirming the origin of earthquakes in the mantle wedge? (b) What challenges remain for interpreting these secondary arrivals? (c) How can the proposed methodology and our findings contribute to future studies aimed at characterizing the subduction system?

Methodology
To constrain the types of secondary arrivals that can be observed in records of SAM earthquakes, we used a two-step approach that simulates waveforms and calculates multiphase traveltimes. First, we calculated and compared synthetic waveforms from a generic subduction zone model ( Figure 2). In principle, the waveform carries information on all discontinuities at which a conversion or reflection took place. When . Velocity-attenuation structure and wavefield simulations of a generic subduction zone (simplified from Hacker et al., 2003). Attenuation is only considered in a subset of simulations described in Table S2. The panels (a-d) show simulations for earthquakes located in different source regions. Red, green and turquoise wavefronts show the wavefield at three snap shots in time, 1.5, 7.5, and 13.5 s, respectively. Red triangles mark the location of three seismic stations (updip, above and downdip from the earthquakes) for which we analyzed seismograms. Different shades of gray represent seismic velocities, which are specified for each region in (a). The inset in (d) shows an example of the spectral element mesh within the corner of the mantle wedge. A movie of the four simulation runs is available at doi. org/10.5281/zenodo.4334041 (see also caption for Movie S1).
we consider multiple interactions with discontinuities (i.e., conversions, reflections, or combined reflection plus conversion of P-and S-waves), the number of all potential secondary arrivals quickly increases beyond what can be realistically observed on the records (e.g., 80 potential arrivals within 10 s for wave interactions with two out of the three discontinuities). Synthetic waveforms can tell us which secondary arrivals we can expect to observe, based on the prominence of each phase. Hence, in a first step, we model synthetic waveforms to analyze which secondary arrivals we can observe. But computing synthetic seismograms is too prohibitively expensive to simulate the high-frequency wave propagation for several hundred earthquakes in a 3-D model. Moreover, the method does not track which interfaces a specific wavefront has interacted with. We tackled these limitations in the second step, where we calculated only the traveltimes of direct and prominent reflected/converted phases for a 3-D velocity-discontinuity model of the WHSZ.

Waveform Modeling Based on Spectral-Element Method
In this first step, we employed the SEM as implemented in the Specfem2D software, which accurately simulates the complete wavefield in 2-D heterogeneous media and allows synthetic seismograms to be extracted at any given point . The SEM is widely used to simulate the wave propagation at global (e.g., Bozdağ et al., 2016), regional (e.g., Zhu et al., 2015), and local scales (e.g., Tape et al., 2009;Tong et al., 2014). At its core, SEM solves the weak form of the seismic wave equation on a finite element mesh in the time domain. It allows complexities such as topography at the free surface, internal discontinuities, anisotropy, (an-/poro-) elasticity, and any lateral variations of elastic parameters and density. The SEM is highly accurate owing to the implementation of high-degree polynomial basis functions. The mesh consisting of individual spectral elements can be distorted to a certain degree, which allows element boundaries to coincide with sharp interfaces. This approach does not require material properties to be interpolated between elements, and it avoids spurious diffractions from staircase-like interfaces .
For modeling wave propagation in a realistic subduction zone, there are three main restrictions regarding the finite element mesh that we needed to address: (a) The mesh needs to represent at least three sharp seismic discontinuities (overriding Moho, plate interface, slab Moho), (b) The mesh needs to limit the spectral elements' degree of distortion to avoid unstable simulations, and (c) The mesh needs to balance the size of individual elements throughout depth to maintain near constant number of points per wavelength and hence optimize simulation run times (cf. , and see text in the Supporting Information). In light of these restrictions, we created a mesh for a generic subduction zone with an average element size of 260 m (maximum: 500 m, see inset in Figure 2d), which allows us to model frequencies up to 6.4 Hz at a time step of 3 ms (see Table S1 for simulation parameters). At this resolution, we were able to adequately model the propagation of a source time function with the shape of a first-order Gaussian-derivative and a dominant frequency of 2.9 Hz. The subduction zone model includes an overriding continental crust and a subducting oceanic crust with low seismic velocities (see Figure 2a for values of Vp, Vs, and density), while the mantle wedge and slab mantle have higher velocities. The subducting crust dips at a 12 angle down to 60 km depth, where it gradually steepens to a dip of 26. At a depth of 80 km, the seismic velocity in the subducting crust increases such that there is only a small remaining contrast relative to the velocities within mantle wedge and slab mantle.

Multiphase Arrival Time Calculation Based on Multistage Fast Marching Method
In this second step, we calculated phase arrival times using the multistage fast marching method that Rawlinson and Sambridge (2004) and de Kool et al. (2006) implemented in the FM3D software. The choice of a second modeling approach is necessary for two reasons: first, while the SEM computes accurate seismograms, it does not provide a means to identify the boundaries at which secondary phases originate; second, to identify arrivals in the data, we need to take into account the 3-D distribution of earthquakes, stations, and geometry of a real subduction zone. Considering these constraints, the multistage fast marching method is optimally suited to compute accurate arrival times for a large number of phase combinations of different events, stations, and interface interactions. The method computes arrival times for multiple phases with any chosen number of reflections and/or transmissions (including conversions) by solving the eikonal equation in a layered 3-D domain with heterogeneous velocity distribution. Its advantages are computational speed and numerical stability. As a method based on the eikonal equation, it yields diffracted phases to all points in a region, even in the presence of LVLs or ray shadow zones (de Kool et al., 2006;Rawlinson & Sambridge, 2004).
For multiphase traveltime calculations in an accurate model of the WHSZ, we constructed a grid comprising four seismic discontinuities embedded in a 3-D background P-and S-wave velocity model (see Figure 1b). The four discontinuities are the surface topography/bathymetry, the overriding Moho, subduction plate interface, and the slab Moho. We defined these discontinuities from the following data sources: topobathymetric data are from the Shuttle Radar Topography Mission (Farr et al., 2007) and the Smith & Sandwell bathymetry database , the overriding Moho model is from CRUST1.0 (Laske et al., 2013, see Figure S4 for an explanation of these choices), and the subduction plate interface is from Halpaap et al. (2019). We defined the slab Moho relative to the plate interface, at the bottom of an 8 km thick slab crust [thickness according to Bohnhoff et al. (2001); Kokinou et al. (2005Kokinou et al. ( , 2006; Pearce et al. (2012)]. For velocities we used the 3-D models from Halpaap et al. (2018).
The accuracy of calculated traveltimes depends on the complexity of the velocity model and increases with grid resolution. Here we chose to sample the velocity model on a spherical grid with 0.02 node spacing horizontally, and 2 km vertically. The grid extends along 401 × 401 × 221 nodes between 35.5  N to 42.5  N and 18.5  E to 25.5  E, padded by another 0.5 in each direction, and from 0 to 200 km depth, padded by another 10 km at the top and bottom. This dense sampling balances computational speed with stability and accuracy for computing the raypaths. Unstable calculations may otherwise abort the traveltime calculation. For 0.02 grid spacing, 6% of ray paths were invalid and <0.03% were missing, while for 0.05 grid spacing, 4% were invalid and 30% missing. Based on the accuracy estimates of Rawlinson and Sambridge (2004, for a 2-D case) and de Kool et al. (2006, 3-D), we expect traveltimes on our 2 km-spaced grid to be accurate to about 10-200 ms (best case vs. worst case) for a typical observation in a subduction zone, at a distance of 100 km.

Analyzing Secondary Arrivals in Synthetic Seismograms
We simulated the wave propagation for two source types-explosion and a set of realistic earthquake mechanisms-at four different locations within the SEM-mesh: mantle wedge, interface, slab crust, and slab mantle (see Figure 2). As earthquake mechanisms, we used four different focal mechanisms that commonly occur in subduction zone settings (e.g., in Greece Halpaap et al., 2019): a normal fault mechanism, a low-angle thrust mechanism, a 45 thrust mechanism (all dipping in direction of subduction), and a perpendicular thrust mechanism. For the 45 we ran simulations both with and without seismic attenuation to analyze possible effects that attenuation may have on the observed waveforms, resulting in a total of 20 simulations (see Table S2 and Movies S1-S3). It is important to note here that we did not seek to fully predict observed seismograms, for which additional source effects (e.g., exact fault orientation, earthquake magnitude, nearsource heterogeneity) play an important role beyond the resolvable level of complexity. In fact, there are three factors that restrict our ability to predict seismograms: (a) the lack of a reference velocity-discontinuity model appropriate for waveform simulations at up to 10 Hz for the highly curved and structure-rich WHSZ (such models exist only in few places such as Japan: JVSIM-model Koketsu et al., 2012;Okamoto et al., 2018), (b) limited knowledge of the moment tensors for individual earthquakes, (c) the events' small magnitudes, which bias the recorded signals toward high frequencies, and (d) computational requirements, which guided our choice of 2-D simulations (rather than 3-D, which would be needed to model arbitrary source mechanisms).
We ran our SEM-simulations for various source types and positions, extracted synthetic displacement seismograms at three stations, and plotted these as station gathers in Figure 3. Here we present a reduced set of 48 seismograms for simulations with the explosion source and the 45 thrust mechanism, a compromise which allows us to demonstrate which phases may realistically be observed, and how these phases differ between source regions. For an exhaustive survey of the results from all 20 simulations, we refer the reader to Figures S1-S3. Each row in Figure 3 (and Figures S1-S3) shows seismograms for one of the two source types (explosion and thrust mechanism, four types in Figures S1-S3), and for the four different source regions (mantle wedge source at the top). The different columns represent seismograms recorded at one of three possible station locations: updip, vertically above, and downdip (see Figure 2). The sources lie straight below the central station, and are vertically separated by 5 km. To highlight coherent phases, we aligned seismograms on the direct P-wave and show the vertical (Z) and horizontal (X) components independently, sorted by event depth (shallowest at the top). We show both explosion (upper panel) and earthquake sources (lower panel), because it facilitates the identification of observed phases: As the explosion source radiates only P-wave energy, only P-to-P and P-to-S interactions appear on the plot (P-S-P conversions are usually negligible). We identified the phases by tracking the wavefronts on snapshots of the wavefields, taken every 0.3 s, that we combined into three movies (available at doi.org/10.5281/zenodo.4334041 and doi. org/10.5281/zenodo.4944977; see captions for Movies S1-S3).
Station gathers of aligned seismograms show phases with a range of moveouts, including (a) positive (delay increases with depth), (b) neutral, and (c) negative moveout. These phases mainly correspond, respectively, to (a) direct S-arrivals and S-to-P conversions, (b) direct P-arrivals and P-to-S conversions (specifically at the overriding Moho), and (c) reflections from either the plate interface or the slab Moho. In Figure 3c, we mark these arrivals, naming them according to the successive interactions (either reflection or transmission) they have undergone. "P" or "S" identifies the wavetype before and after an interaction, denoted by a letter in subscript: m -subducting oceanic moho; t -slab top (=plate interface); M -overriding plate Moho; T -surface Topography (see also Figure 1a).
Visible phases change across stations and depending on the source region of the earthquake, but there are usually only a few prominent secondary arrivals. In most setups, direct P-and S-waves exhibit the largest amplitude, but there are exceptions. If a station is located such that a vanishing amount of energy is radiated along the direct raypath, then there are always secondary arrivals with substantial amplitude (see Figure S7). This is due to raypaths interacting with dipping slab discontinuities (i.e., not the overriding Moho) taking off from the source at an angle that differs considerably from the take-off angle of the direct ray. Generally, the highest amplitudes for secondary arrivals occur for (a) reflections rather than conversions, (b) interactions with a discontinuity close to the source, and (c) S-to-S or S-to-P interactions rather than P-to-P or P-to-S interactions. Typically, the secondary arrivals that are easiest to spot on single traces are reflections off the plate interface or the slab Moho. Conversions, for their part, are better observed on station gathers as they tend to exhibit a strong coherence from trace to trace. These general observations are supported by the analysis of all 20 simulations, which include other types of source mechanisms. Different earthquake sources lead to individual secondary phases becoming stronger for some incidence angles, but they do not change the general patterns.
Individual waveforms vary considerably as a function of source region. Below, we list the main waveform attributes for each source region.
Mantle wedge earthquakes show: 1. Shortest S-P time, that is, shorter than any slab earthquake when the slab dip is shallow ( 25) and source and receiver raypaths are similar, 2. Longest S-coda, due to reflections off the interface and the slab Moho, 3. Reverse polarity for the earliest secondary arrivals (P t P in this case) compared to the direct P-arrival (this is unique to this source region).
Interface earthquakes show: 1. Longer "quiet"-time after the direct arrivals (P and S) and after the earliest secondary arrivals, 2. Equal polarity for the earliest secondary arrivals compared to the direct P-arrival.
Slab crust earthquakes show: HALPAAP ET AL.

Observing Secondary Arrivals in Field Data
We now shift our attention to assessing how secondary arrivals can be used to deduce the source region of SAM earthquakes from field data. Here, our primary targets are earthquakes at depths between 40 to 90 km -for two reasons: First, this depth range is best covered by our test data set from the WHSZ; second, in this depth range, the subducting crust is not yet being transformed to eclogite (cf. Figure 3 in Halpaap et al., 2019). The latter implies that strong velocity contrasts exist at both the bottom and top of the subducting crust in this depth range. These contrasts are a prerequisite, without which we would not expect any difference between waveforms of earthquakes above and below these boundaries.

Data Set
The data set consists of seismograms recorded across the Peloponnese peninsula in the WHSZ (see Figure 4). This location is well-suited to this study because it is well-instrumented and there is strong evidence that earthquakes occur in the various source regions of the subduction system. The seismograms were recorded by temporary stations of the Medusa and Egelados experiments across the Peloponnese-Attica region in Halpaap et al. (2019) provided precise hypocenters for these earthquakes from relocation using the double-difference method in a 3-D velocity model (see linked data sources in data and materials availability statement). The relocation yielded hypocenter locations that are accurate both relative to one another (average relative error: 0.04 km) and in an absolute sense (average absolute error: 1.7 km). Halpaap et al. (2018) also remigrated the image from scattered teleseismic waves with a 1-D velocity model based on the 3-D model used for relocation. The remigration ensured structure imaged with teleseismic waves matches the local seismicity patterns. Together with a model of the top of the slab (Halpaap et al., 2019), the accurate hypocenters allowed them to calculate the distance between each earthquake and the plate interface to within 1 km. The seismicity catalog contains more than 200 P-picks at each one of 11 stations situated above the earthquake cluster (on average 375 picks per station, from a total of 468 events). These 11 stations (S009-S016, PE02, PE05, PE07, black inverted triangles in Figures 4c and 5) show varying noise levels, but together they were able to record events ranging in local magnitude M L between −0.5 to 2.9 (see Figure S6 for noise analysis). Another three stations (PE04, S017, LTK) above the earthquake cluster yielded too few high-quality recordings to be used here.

Waveform Processing
Raw seismic data recorded in the field contain a combination of signals, including signals from primary and secondary phases (the phases of interest), incoherent signals from local scattering at small-scale heterogeneities, and signals from sources other than earthquakes. The two latter signals, often considered noise, HALPAAP ET AL.
10.1029/2020JB021573 12 of 28 mask the phases of interest and the key characteristics we see in the synthetics. To make secondary phases visible, we processed three-component seismograms following a seven-step workflow. We relied on the functionality of the Gismo seismic data analysis toolbox (Reyes & West, 2011;Thompson & Reyes, 2017) and implemented new processing and plotting schemes. The workflow consists of the following steps, illustrated in Figure 6 on the field data (see also Figures S1-S3 for processed synthetic seismograms): 1. Bandpass filtering between 1.5 and 10 Hz 2. Alignment of waveforms based on arrival picks and cross correlation 3. Rejection of recordings with SNR <2.5 4. Rotation of channels from ZNE (vertical, N-S, E-W) to ZRT (vertical, radial, transverse) 5. Polarization filtering on three-component data 6. Automatic gain control on three-component data 7. Sorting of waveforms according to distance from the plate interface and plotting of waveform envelopes

Bandpass Filter
The frequency content of microseismic events (M < 3) is limited depending on the event's magnitude. Seismograms of such small events are strongly affected by noise, which also varies with frequency. Our analysis of the signal-to-noise ratio (see Figure S6) showed that noise level is highest in the microseism frequency band below 0.4 Hz, and that it decreases up to a frequency of 5 Hz, where it flattens out toward higher frequencies. Larger earthquakes have a higher SNR and longer source time function, which allows them to be observed down to lower frequencies. While noise dominates at low frequencies, high frequencies are likely dominated by wave scattering. To allow major converted/reflected arrivals to stand out, we filter records within one universal pass band that best suits all data, between 1.5 and 10 Hz.
We considered deconvolving the seismic traces to avoid misidentifying source time function-generated reverberations as secondary arrivals. However, where feasible in this study (for M  2.5), deconvolving waveforms yields no advantage that would allow secondary arrivals to be recognized in our study's data set. In the data set, earthquakes range in magnitude M L between −0.5 and 2.9. Even with conservative empirical estimates of the source time function scaling with magnitude (i.e., very low stress drop and above-average source time function length), all earthquake's source time functions here are shorter than the minimum time difference between arrivals that we considered (considering the scalings from Hiramatsu et al., 2002;Prieto et al., 2004, see Text S3 for a detailed explanation).

Rotation
To best separate the energy of incoming P-and S-waves, we rotated the seismograph components from ZNE (vertical, north-south, east-west) to ZRT (vertical, radial, transverse). Rotation to ZRT accounts for the backazimuth between station and event, but not the incidence angle. It does not yield complete separation of P-and S-arrivals for rays incoming at oblique paths, but achieves good separation for subvertical paths. To achieve complete separation, we would have to rotate to the P-SV-SH system (or at least LQT, see Rondenay, 2009). Such a rotation was impractical here, because the incidence angle changes substantially between direct and converted/reflected arrivals: For the direct rays, the average absolute incidence angle is 35. For secondary arrivals from the same source, the incidence angle is on average 6 smaller, with average minimum and maximum differences of −20 and 11, respectively. We also considered the horizontally polarized SH-waves generated by earthquakes, as these become coupled to the P-SV system in the presence of conversion/scattering at 2-D and 3-D structures and anisotropic materials (e.g., an layer of serpentine near the plate interface inferred by Olive et al., 2014). With the slab discontinuities dipping at about 21 , with some anisotropic layers, and with sources and receivers distributed in 3-D, we expected conversions from and to SH-waves to influence all components.

Polarization Filter
Scattering at small-scale heterogeneities is a major source of signal disturbing the identification of secondary arrivals. Compared to scattering and random noise, secondary arrivals generally show a particle motion more rectilinear rather than elliptical (Montabetti & Kanasewich, 1970). The particle motion of a signal can be taken into account with three-component seismograms, through application of data-adaptive multi-component filters, including the product of two components (Regnier et al., 1994) or a polarization filter Figure 6. Stepwise demonstration of processing scheme. A selection of seismograms from mantle wedge earthquakes recorded at station S010 are shown. Rows 1-7 represent processing steps 1-7, with alignment step applied twice: after bandpass-filtering and again on rotated seismograms after polarization-filtering. The first three rows display ZNEcomponents, while the following rows display ZRT-components. In row 7, annotated phase names and colored lines mark arrivals computed from fast marching method (green: direct P-, S-phases, blue: purely converted phases, red: purely reflected phases, purple: reflected and converted phases). (Montabetti & Kanasewich, 1970;Samson & Olson, 1981). Both methods increase signal energy when the particle-motion rectilinearity is high. The product between pairs of channels is a useful tool to concentrate information from multiple input-channels onto one output channel, and has been employed in the identification of local converted phases in the Chilean subduction zone (Regnier et al., 1994). A polarization filter is better suited to partition signal energy between an equal number of input-and output channels, which allow differentiation of P-, SV-and SH-arrivals. Here we implemented a polarization filter as suggested by Montabetti and Kanasewich (1970).  and Horleston and Helffrich (2012) successfully applied such a filter to local recordings of slab earthquakes. Below, we give a short summary of the filtering strategy and chosen parameters.
The polarization filter iterates through a three-component signal and calculates two gain functions: for rectilinearity and directivity (Montabetti & Kanasewich, 1970). Here, rectilinearity is a scalar measure of polarization degree, while directivity is a three-component measure indicating which direction (in reference to the components) a polarized signal dominates. Both gain functions are calculated from the eigenvalues and eigenvectors of the signal's covariance matrix, which itself is computed for a small time window of the three-component signal. The rectilinearity gain function is calculated as the ratio between largest and intermediate eigenvalue of the covariance matrix, while the directivity gain function is obtained from the largest eigenvector. Montabetti and Kanasewich (1970) introduced three parameters (n, J, K) to weight the eigenvalue ratio, the rectilinearity, and the directivity, which we chose according to their suggestion (n = 0.5, J = 1, K = 2). Alternative choices appeared to have negligible effect. Before application, the gain functions are smoothed within a window of a chosen length -a crucial user-defined parameter. If the chosen window is too short, random maxima in the gain functions may occur; if it is too long, a long wave train around secondary arrivals is boosted. For the window length, Montabetti and Kanasewich (1970) suggest twice the dominant period. With the aim of consistent processing, the dominant period is not an obvious measure for a data set with earthquakes of varying magnitude. Through trial-and error, we found that a length of 0.5 s is a generally sensible choice for our data set, independent of the bandpass filter frequencies.
A window length of 0.5 s is close to twice our 1.5-10 Hz-passband's central period of 0.258 s, which would agree with the suggestion of Montabetti and Kanasewich (1970). However, when changing the passband, a value of 0.5 s remained optimal. Particularly for higher frequencies, it appeared that the signal is dominated by scattering, for which a smaller window length leads to undesired amplification.

Waveform Alignment
The better the alignment between waveforms in a record section, the easier it is to distinguish coherent arrivals across traces. While we estimated our manual arrival picks to be accurate within 0.1-0.2 s for lownoise records (likely larger for high-noise records), even such small errors may have led to an offset of multiple cycles between traces. Through cross-correlation of traces, we aimed to align waveforms to within 0.05 s difference of their P-arrival. We performed the alignment twice, based on the P-arrival on the Z-component waveform-first on the envelope of wider-band (1.5-15 Hz) waveforms, and second on the envelope of the narrower-band (1.5-10 Hz), polarization-filtered waveforms. With this two-stage strategy, we aimed to avoid a series of common issues that may arise when comparing different events: (a) Alignment should be on the P-arrival onset rather than the maximum of the source time function (hence a wider passband at first), (b) the alignment should be independent of the first motion polarity (hence envelopes), and (c) the alignment should not be affected by pre-arrival ringing introduced by a finite impulse response filter (hence repeated alignment after boosting rectilinear signal). We cross-correlated every waveform within a gather, over a window ranging from −0.2 to 0.5 s relative to the picked P-arrival, to find an optimal lag time between each waveform pair. We then adjusted the P-arrival times of each waveform relative to the waveform with the smallest average lag, but we did not allow arrival shifts when correlation yielded lag values of more than half the central period (i.e., 0.258 s).

Automatic Gain Control on Three-Component Data
When secondary phases are invisible purely because of their weak amplitude and not because they overlap with other phases, then an amplitude correction can make them visible. In reflection seismology, an automatic gain control (AGC) is often used to balance reflection amplitudes, with the trade-off that relative amplitude information is then lost. An AGC normalizes amplitudes in a moving window across a trace. Here we tested both a single-channel AGC and a three-component AGC. The latter uses the mean amplitude of all traces to compute a correction factor at each time step, which is then used as a weight on the three components. The tests showed that, with a window length of 2 s, the three-component AGC accomplished two desired effects: first, it allowed secondary arrivals to be better observed between the direct arrivals and at long time offsets after the S-wave; second, it emphasized the component (Z, R, or T) on which a secondary phase is dominant by reducing the phase's relative amplitude on the other two components.

Final Event Selection
After completing the preprocessing steps described above, out of an original set of 4,120 observations (station-event pairs with a P-pick), 2,885 three-component seismograms were retained. For the 11 selected stations, the number of retained three-component seismograms ranged between 123 (station S015) and 359 (station S010). The processing steps that removed seismograms from the data set are: rotation, polarization filtering, and the check for SNR. While the former two steps required that three-component recordings were complete and without gaps, the latter step rejected recordings with high noise (i.e., SNR <2.5 on all three channels, computed from two windows of 2-s length each, one preceding the P-arrival by 0.5 s and one starting at the P-arrival). For computed arrivals, we chose to plot all arrivals arising from one single interaction with a discontinuity (i.e., reflection, conversion, or combined reflection plus conversion). This restriction is warranted by the synthetic modeling (see guidelines in Section 4) and by the actual observations.

Phase Identification
In panel (b), an additional plot indicates two amplitude ratios between the wavelets of the direct arrivals-we show one ratio between the S-arrival on the R-component and the P-arrival on the Z-component (approximately the SV/P-ratio), and another one between S-arrivals on the R-and T-components (approximately the SV/SH-ratio). We computed these ratios from the maximum envelope amplitudes in a window −0.3 to 0.7 s around the arrivals. For the P-wavelet we used the arrival time from the relative seismogram alignment, while for the S-wavelet we used the picked S-minus-P arrival when the misfit between picked and computed times was <0.3 s, otherwise we used the computed arrival. As these measurements are very sensitive to noise (for small events) and imperfect component separation, we applied a running median filter across 4 traces to better extract systematic changes across the station gather. Four stations (S012, S013, PE05, S016) show clear amplitude ratio differences across the major discontinuities above/below which earthquakes occur-that is, the plate interface ( st d  0 km) and the slab Moho ( st d  −8 km).
To identify secondary arrivals in the processed station gathers, we visually search for peaks in the seismogram envelopes that appear consistently across at least 10 traces, and match computed arrival times of secondary phases 0.5 s. We rate each identified arrival on a three-level quality scale: A-rated phases are clearly visible and their identification is without reasonable alternative; B-rated phases are visible and identified with their most probable arrival name, and C-rated phases are least well-observed, being either only faintly visible and/or their identification is unreliable (see Figure 8 for a selection of seismograms where we identified A and B-rated phases). The secondary arrivals and information regarding their identification are summarized in Table 1. Additional major phases that we cannot explain from any predicted arrival (e.g., station S010, at 0.8 s) are described in the comment column. For each identified phase, the table also lists the potential source region of the associated earthquake.

Results -Observed Phases
For each of the 13 analyzed three-component station gathers (Figures 7 and S7-S27), we find on average two to three clear secondary arrivals that are listed in Table 1. As an example, for station S010 we clearly observe on the Z-component an S M P-conversion at the overriding Moho that we identify with high certainty. Even though this arrival appears at a similar time as the S m P-reflection from the slab Moho for earthquakes on the plate interface, it exhibits a characteristic S M P moveout across traces for earthquakes in the slab crust (especially 1-5 km below the plate interface). Compounded with the fact that no apparent S m P are observed for earthquakes in the slab crust, this points to a S M P origin. We observe three more secondary arrivals on the other two components which we identify with a "most probable" quality level as: (a) an S M P conversion at the overriding Moho (visible on T-component, in particular for interface earthquakes), (b) an S t S-reflection off the interface (on T-component, for mantle wedge earthquakes), (c) an S m P-reflection off the slab Moho (on T-component, for slab crust earthquakes). With faint visibility we observe another secondary arrival on the Z-component that we identify as likely being an S m S-reflection of the slab Moho (for earthquakes on the interface). As listed in the table's comments, the Z-component seismograms also show a secondary arrival shortly (<1 s) after the P-arrival that does not match any secondary arrivals that we expect.
Generally, the clearest secondary arrivals are the S-to-P conversion at the overriding Moho (S M P, e.g., stations S010, PE02), and a P-reflection from the slab Moho (P m P, e.g., stations S015, PE02), both observed predominantly on the Z-component. In optimal conditions, these arrivals are consistent with the predicted times. They appear primarily on the expected component (Z for incoming P-waves), and show the characteristic moveout compared to the P-arrival-that is, either negative (reflections) or positive (S-to-P conversions). The latter of the two phases, the P m P-reflection, is particularly well-observed because it appears highly coherent and with little shift across multiple traces. The traces that show this arrival (e.g., stations S015 and PE02) all originate from earthquakes we previously located within 1 km normal distance of the subduction interface.
HALPAAP ET AL.   . We identified phases marked in (b, d, and e) with high confidence, and the phase in (c) with medium confidence (see Table 1). Theoretical arrivals of other potential phases are also marked, but not labeled (see legend for color coding, and Figures 7 and S7-S27 for complete station gathers with labels).
Other secondary arrivals require a certain level of seismogram interpretation, either because their amplitude varies between traces, or because they do not align as well as others. These include S t S-reflections off the plate interface (e.g., station S009, Figure 8b) particularly late S m S-reflection off the slab Moho (e.g., station S014, Figure 8c). The arrivals with the highest uncertainty are those with a combined reflection and conversion, which often overlap with other arrivals in the P-and S-coda, and which are particularly difficult to identify in the S-coda. These include a P t S-reflection off the plate interface (stations S013, S016, PE07) and a P m S-reflection of the slab Moho (stations S012, PE02, S013, S016, PE07), for which we rate most identifications with the lowest quality rating C. In addition to the identified arrivals, on some station gathers we observe clear secondary arrivals that do not match any expected secondary arrivals. Such unexpected arrivals are particularly visible at stations S010 (0.8 s after P, see Figure 7, arrow 1) and S013 (1 s after S, Figure S19).
The amplitude ratios between direct P-and S-arrival wavelets show no clear trends at nine out of 13 stations, but at the remaining four stations, we observe clear contrasts between sets of earthquakes. These four stations are PE02, S012, S016, and S013, with the latter ( Figure S19) showing a very clear amplitude ratio contrast of 50:1 between earthquakes that we previously located within <1 km distance from the plate interface, and earthquakes located farther below the plate interface.

Discussion
The focus of this study is to assess whether the analysis of seismic waveforms can help support the existence of mantle wedge earthquakes. While it is generally accepted that subduction earthquakes occur in the slab crust, in the slab mantle and sometimes even on the interface down to 90 km depth (Abers et al., 2013;Uchida et al., 2016), earthquakes occurring within the mantle wedge have proven more controversial. Generally, the controversy stems from the fact that mantle wedge seismicity is rare. Only a small, but growing set of targeted high-resolution studies in the last decade have been able to distinguish mantle wedge earthquakes from those in the slab (Bie et al., 2019;Chang et al., 2017Chang et al., , 2019Davey & Ristau, 2011;Halpaap et al., 2019;Laigle et al., 2013;Malusà et al., 2016;Nakajima & Uchida, 2018;Nakata et al., 2019;Paulatto et al., 2017;Uchida et al., 2010;White et al., 2019). In the specific case of western Greece, the putative mantle wedge earthquakes overlie an aseismic portion of the slab, raising questions as to whether these are, in fact, mislocated slab earthquakes (i.e., upward shifted). An apparent alignment of earthquakes at a dip of 45 (much steeper than the imaged LVL dip of 21), raised further uncertainty. Previous hypocenter solutions yielded small absolute errors (<1.7 km), on a scale similar to the resolution of the images from scattered teleseismic waves. However, imaged discontinuities and earthquake locations relied on independent data-that is, teleseismic waveforms and arrival time picks from local waveforms. As we have seen from the synthetic models, these locally recorded waveforms carry secondary arrivals which may confirm where an earthquake has occurred in relation to the discontinuity.

Criteria to Identify Mantle Wedge Earthquakes
Both the synthetic and observed (processed) seismograms show promising signals that help determine the source region of subduction earthquakes. The key observations extracted from the waveforms are the amplitude ratios and multiple secondary arrivals. Together, these observations support the case of earthquakes in the mantle wedge below Tripoli.
Amplitude ratios of direct P-, SV-, and SH-arrivals, measured at different azimuths and takeoff angles, largely depend on an earthquake's radiation pattern and hence its focal mechanism. For a confined cluster such as below Tripoli, azimuths and takeoff angles vary between stations, but are nearly constant for all clustered earthquakes observed at a single station. Halpaap et al. (2019) showed that earthquakes in the mantle wedge, on the interface and in the slab have systematically different focal mechanisms, with normal faulting dominating in the wedge, slab-parallel thrusting dominating on the interface, and trench-parallel thrusting dominating in the slab. While the focal mechanism analysis included a larger data set spanning 11 years, for the earthquakes selected here we were only able to retrieve focal mechanisms for 12 out of 468 earthquakes from the 2006-2007 period. Only these 12 earthquakes yielded focal mechanism solutions because they were particularly well sampled (>15 distributed stations) and of a magnitude M L >1.6 that allowed first arrival polarities to be picked. Although focal mechanisms for smaller magnitude events could not be obtained, the amplitude ratios of all events in the data set are consistent with three systematically different source mechanisms. Stations PE02, S012, S013, and S016 exhibit a clear contrast in amplitude ratio between at least two of the three possible source regions defined by their varying focal mechanisms-that is, the wedge, interface and slab. Compared to other stations in the network, these four stations appear to be favorably located to record large variations in amplitude ratios. Such extreme amplitude ratio variations occur where ray paths take off from the source, either close to the nodal planes or in the center of the compressional or tensile regions, for at least one of the three prevalent focal mechanisms (see Figure S5 and Text S4 for details).
The secondary phase observation that best supports the existence of mantle wedge earthquakes is an S-reflection from the interface (unique to these earthquakes), combined with a very late S-reflection from the slab Moho (cf. Figures 8b and 8c). Reflections from the slab Moho are common in our observations, especially P m P but also to a lesser degree S m S. For an earthquake in the mantle wedge, when a reflection from the Moho is observable, then a reflection from the interface (S t S) is also usually observable (see, e.g., Figure 8b), since both phases take off from the source at similar angles, resulting in minimal difference in the amount of energy radiated from the source toward the two raypaths. For these earthquakes, the S m S wave also becomes the most delayed secondary arrival after the direct P. To cause such a late arrival, the earthquake needs to be far above the deepest discontinuity, that is, in the mantle wedge. Shallower earthquakes in the overriding plate could also be responsible. However, based on our experience from arrival picking, these are, these are generally easy to distinguish from deeper earthquakes based on scattering-dominated waveforms and different first-arrival moveout. Arrivals with even larger delay times can only be caused by multiple reflections in a layer such as the crust, or by reflections from a discontinuity below the slab Moho. As previously discussed, the latter option is unlikely because such reflections have never been observed.

Mismatches Between Modeled Arrivals and Observations
Our systematic search for secondary arrivals largely confirms the expected earthquake locations and the assumed velocity-discontinuity model for the WHSZ, but it also points out some model-data mismatches. It is generally difficult to verify phase identifications across a network as observed arrivals vary between stations. This variability arises from the station geometry of the Medusa network, which follows a SW-NE line with stations spreading out in a triangular shape toward the back arc. Designed for receiver function imaging, this spread-out configuration was chosen because the dip direction of the slab was not well known prior to the deployment. While the average station spacing along the profile is 8.5 km in the forearc, the profile-normal spacing is up to 40 km. With such separation, the recorded seismograms are expected to vary considerably from station to station due to 3-D complexity of source and discontinuities. We expect that studies such as ours may benefit from short inter-station distances along a narrow line, aided by a few off-line stations for 3-D coverage.
In this study, we tackled the challenge posed by the network by taking a 3-D discontinuity-velocity model into account. The least well-defined interface in the model is the overriding Moho, because we could not deduce its 3-D structure from previous analyses (2-D receiver function migration and 3-D local earthquake tomography) and had to rely on a global model. Of all conversions, those from the overriding Moho may hence have the highest uncertainty in the arrival time. Station PE07 shows Moho conversions particularly well. It shows that the observed S M P-conversions arrive earlier than predicted, while the P M S-conversions arrive later than predicted. Together, these indicate that the true location of the Moho is shallower than in the model. Even if predicted phases match the secondary arrivals well, there is always a possible tradeoff between the location of a discontinuity and the velocity model. One example where we suspect such a tradeoff to be important is in the low velocity layer. P m P-reflections from the slab Moho, which we consider as the discontinuity at the bottom of the LVL, are particularly well-observed at stations S015 and PE02. The reflected arrivals are a good match to the predicted arrival times, which we computed for an 8 km thick slab crust with velocities imaged in our tomographic study. However, in tomographic studies, small features with strong velocity contrasts always suffer from smearing and smoothing. It is thus possible that the ve-locity in the slab crust is lower than imaged, which in turn would require a thinner LVL to fit the observed reflections (cf. Shiina et al., 2017;Song et al., 2009).

Unknown Discontinuities
Clear secondary phases that we could not match to expected arrivals (e.g., station S010, 0.8 s after P, and S013 (1 s after S) suggest that there is at least one discontinuity that has not yet been resolved in imaging studies of the WHSZ. These early secondary arrivals are best observed for earthquakes close to the interface. The clustered interface earthquakes extend along a plane that measures 8 km in the dip direction and 5 km in the trench-parallel direction. As all these earthquakes show the arrival at a similar time, it appears that: (a) The earthquakes actually occur on one very narrow fault plane, which is narrower than the hypocenter solutions would suggest (1 km), (b) the arrival stems from a discontinuity (hereafter called "X") that parallels the interface close to the top of the slab (cf. Bostock, 2013). Unexpected arrivals may also be produced by structure below a seismograph station. In that case, all SAM earthquakes recorded at that station should show the arrival at the same relative time, and the arrival should appear as a P-to-S conversion on a horizontal component. Conversions at shallow structure should be uniform across channels, as the difference in ray paths from the clustered SAM seismicity is minimal along the shallowest part of the ray. However, as the early arrival at station S010 appears predominantly on the Z-component and is only observable for earthquakes close to the interface, we consider this early arrival to be caused by a deep structure.
The X-arrival may be due to layering within the oceanic crust, as has been documented in wide-angle studies in the oceans (i.e., Christeson et al., 2007) and in the Alaskan subduction zone (D. Kim et al., 2019). This layering typically includes the following units: layer 1, consisting of sediments; layer 2, consisting of hydrated pillow basalts and sheeted dikes; and layer 3, consisting of low-porosity gabbros (see e.g., Grevemeyer et al., 2018). The thickness of layer 1 in the slab (subducted sediments) ranges from 0.2 to 1.4 km, while layer 2 thickness ranges from 1.5 to 2.0 km (Syracuse et al., 2010). Layer 2 is further divided into layers 2A and 2B, which are separated by a discontinuity representing either the boundary between lava and dykes, or an alteration boundary (Christeson et al., 2007). Seismological studies on the structure of the slab can often resolve only one composite LVL, as in Figure 5. Competing models attribute the low-velocity signature to different parts of the system, including (a) the complete oceanic crust, (b) a fluid-saturated, overpressured layer confined to the upper crust, (c) a distributed plate interface shear zone, (d), a layer of underthrusted sediments, or (e) a layer of serpentine at the top of the slab (Bostock, 2013;Chuang et al., 2017).
With only station S010 clearly showing an arrival at 0.8 s, we cannot tell with certainty whether it is an S X P-conversion, a P X P-reflection, or an S X P converted reflection. In any case, the boundary is located no farther than 3 km from the cluster of interface earthquakes when considering conservatively high velocities ( P V  7.5 km/s, S V  4.3 km/s). For lower velocities, that is, as low as P V  6.5 km/s (Shiina et al., 2017) and S V  2.0 km/s (Song et al., 2009), the layer could be as thin as 1.3 km if the secondary arrival stems from an S X P-reflection. Assuming that the earthquakes occur at the top of the subducted sediments, a layer of 1.3-3 km thickness below may consist of: (a) only layer 1, (b) layer 1 and layer 2A, (c) layer 1 and layer 2. Obtaining better constraints on the properties of this thin layer will be key to better understanding the occurrence of earthquakes on the interface and in the mantle wedge, and to coupling these findings to fluid flow and slab structure. These tasks were, however, beyond the scope of this paper. Future studies will have to identify the X-arrival on multiple stations, which would allow for a thorough analysis of its moveout and thus better estimates of the upper layer's thickness, seismic properties, and position relative to earthquakes.

Pitfalls of Using Secondary Arrivals
Our analysis exposes some possible pitfalls of using secondary arrivals from local SAM earthquake data to identify the source region of these earthquakes. For all source regions that offer several converted/reflected raypaths between source and receiver, there is often a risk of misidentifying phases. The risk is minimal when analyses are restricted to earthquakes that occur well below the slab Moho (i.e., slab mantle) -a strategy adopted by many authors (e.g., Horleston & Helffrich, 2012;Matsuzawa et al., 1990;. Conversely, the risk is large especially for earthquakes that occur in the slab crust, for which reflected S m P phases off the subducting Moho may be mistaken as S M P conversions from the overriding Moho. When it is not feasible to consider only earthquakes from the slab mantle, then additional constraints need to be taken into account, such as the polarity of the phase (e.g., Horleston & Helffrich, 2012) and/or the energy separation between components. The data in our study showed that the polarity of secondary phases could not be viably measured for earthquakes with 3.0 M  . For the events studied here, the energy separation between components did not always match expectations. It appears that the complexity in structure (e.g., dipping interfaces) and anisotropy strongly affect records, leading to what is an expected P-to-SV-converted phase sometimes arriving on the T-component (where only SH-energy is expected) rather than the R-component.
To avoid such pitfalls and misidentification of phases, it is important to understand the timing and prominence of different secondary arrivals. While the fast marching method can accurately predict arrival times for 3-D models of velocity structure and discontinuities, the method does not resolve which few out of about 80 considered secondary phases are strong enough to be observed. One solution is modeling the wavefields for a set of earthquakes and analyzing which secondary arrivals are observable for a specific setup. This setup may differ in receiver locations, source parameters (location, mechanism), and the structure of the subduction zone, including the geometry of discontinuities. But subduction zone structure is commonly not known well enough to directly compare modeled and observed seismograms at the high frequencies relevant to earthquakes with M < 3. At these frequencies, modeling the wavefields is also computationally expensive, hence some simplification is necessary. Here, we have shown that a suitable simplification is to employ the 2-D spectral-element method to model the wavefield in only two dimensions, and for a set of only four earthquake locations with the same mechanism. This strategy allowed us to conclude that we should expect as observable secondary phases mainly S-to-P conversions, P-reflections, S-reflections, and S-to-P converted reflections, all with an interaction at only one discontinuity. This is an approximately 10fold reduction in secondary phases in comparison to the 80 considered phases, without which the identification of observed phases would not be possible.

Conclusions and Outlook
In this paper we used reflected and converted arrivals in the seismograms of SAM earthquakes to distinguish between earthquakes that occur in the mantle wedge, on the interface, and in the subducting crust. We computed synthetic high-frequency seismograms to show how common station gathers of P-aligned seismograms from multiple events can be used to recognize secondary phases based on their timing and moveout. We found positive, neutral, and negative moveouts (where positive equates a traveltime increase with earthquake depth) that are characteristic for the following arrivals: (a) positive: S-arrivals and S-to-P conversions; (b) neutral: direct P-arrivals and P-to-S conversions; (c) negative: reflections from the plate interface or the slab Moho. For observed three-component seismograms we developed a seven-step processing sequence to enhance the signal-to-noise ratio of secondary arrivals: (1) rotation, (2) bandpass filtering, (3) selection of high-quality records, (4) polarization-filtering, (5) cross-correlation based alignment, (6) three-component gain control, and (7) plotting of the waveform envelopes together with computed phase arrivals in common station gathers, sorted by the source's relative distance from a discontinuity (the plate interface). To calculate theoretic arrival times of secondary phases, the fast marching method proved to be an efficient tool that can take a 3-D velocity-discontinuity model into account. We showed that the workflow succeeds in characterizing earthquakes of varying magnitude and signal-to-noise ratio, but that the observation of a secondary phase depends critically on the source's radiation pattern. Applying the workflow to a cluster of deep seismicity in Greece and analyzing the resulting station gathers, provided the following information: 1. We identified on average two to three secondary arrivals per station. 2. Reflections from the plate interface and the slab Moho provide strong evidence that there are indeed earthquakes occurring in the mantle wedge beneath the Peloponnese. 3. Amplitude ratios of direct P-and S-waves showed that focal mechanisms differ systematically for earthquakes located above, on, and below the plate interface. 4. There is a previously unidentified discontinuity 1.3-3.0 km below the top of the slab, which may represent the boundary at the base of either the subducted sediments, or oceanic crust layers 2A/2B.
Motivated by these observations, we propose to exploit secondary arrivals for mapping the overriding Moho, as well as the LVL's location, thickness, and internal structure. Improved constraints on the internal properties of the slab and the mantle wedge would allow us to better understand the distribution of pore fluid in the LVL. Prime targets for this type of study are sites such as the earthquake cluster below Tripoli, where fluid escapes the subducting crust through a vent on the interface and causes earthquakes in the mantle wedge. Compared to other subduction zones with pronounced double seismic zones (e.g., New Zealand, Chile, and Japan), in Greece the task would be complicated by the lack of earthquakes deep within the subducting plate. Based on traditional converted-wave methods (i.e., Regnier et al., 1994;, in this case, secondary arrivals would primarily allow mapping of the overriding Moho [cf. Figure 8e, similar to Uchida et al., 2010 in Japan]. However, with the reflected and converted secondary arrivals that we have newly identified in this study, adequate ray tracing and/or migration of reflections would also afford mapping of the plate interface and the subducting Moho, and inversion for seismic velocities. Ultimately, resolving these parameters would give us insight into the processes that promote fluid migration and the associated generation of earthquakes in subduction zones.  Figure 1b were downloaded from the United States Geological Survey at https://dds.cr.usgs.gov/srtm/. Bathymetry data for that model are from the Smith and Sandwell bathymetry database at the Scripps Institution of Oceanography, University of California at San Diego, at http://topex.ucsd.edu/marine_topo/. Supplementary movies S1-S3 show wave propagation simulations and are available at doi.org/10.5281/zenodo.4334041 and doi.org/10.5281/ zenodo.4944977. We performed wave propagation simulations with the SpecFem2D code that is, available from http://github.com/geodynamics/specfem2d. We computed multiphase arrival times with the FM3D code that is, available from http://www.iearth.org.au/codes/3Dfastmarching/. Maps were created with the m_map Matlab toolbox (http://www.eoas.ubc.ca/∼rich/map.html). The GISMO code for seismic data processing in Matlab is available from Zenodo under http://doi.org/10.5281/zenodo.1404723, and the newest version of GISMO is available at https://github.com/geoscience-community-codes/GISMO. A modified version of the GISMO code, including the scripts and selected data needed to reproduce the seismogram processing and plotting can be accessed under http://doi.org/10.5281/zenodo.4332370.

Acknowledgments
We thank the Associate Editor and two anonymous reviewers for thoughtful and constructive reviews, which helped improve the initial version of the manuscript.