Attenuation contrast in mantle wedge across the volcanic front of northeastern Japan that controls propagations of high-frequency S-wave later phases

Distinct later phases of waves with rich high-frequency (> 8 Hz) components were observed for intraslab earthquakes that occurred at intermediate depths, particularly at depths exceeding 100 km, in the northeastern (NE) Japan subduction zone. These high-frequency later phases (HFLPs) showed anomalously large peak-amplitude delays, up to ~ 50 s after direct S-wave arrivals at stations in the backarc region. Using a source-scanning algorithm, we investigated the locations of passing points affecting the propagation of HFLPs. The passing points were estimated to be in the forearc region in the entire NE Japan, indicating that HFLPs are scattered waves that pass through the forearc region. The propagating HFLPs seem to bypass the backarc mantle wedge, as a consequence of the distinct attenuation contrast in the mantle wedge across the volcanic front in NE Japan. These HFLP observations suggest that the high-attenuation zone in the backarc mantle wedge controls propagations of the high-frequency waves of intraslab earthquakes, in addition to the scatterers possibly locate in the forearc region.


Introduction
In northern (N) and northeastern (NE) Japan, the Pacific plate is subducting beneath the Kuril and NE Japan arcs. The Pacific plate subducts nearly perpendicular to the Japan Trench in NE Japan, whereas obliquely at the Kuril trench in N Japan (Fig. 1a). The oblique subduction of the Pacific plate is driving the southwestward migration of the Kuril arc, which collides with the NE Japan arc at the central part of Hokkaido, Japan (e.g., Kimura 1996).
Seismograms observed for earthquakes occurred at both N and NE Japan record not only direct P and S waves, but also a variety of other signals reflecting the highly heterogeneous structures related to the subduction of the Pacific plate. In the signals of intermediate-depth intraslab earthquakes occurring at depths of ~ 80-250 km, the direct S waves are usually followed by signals with frequencies higher than ~ 8 Hz. These high-frequency later phases show remarkable amplitude peaks, whereas their arrival times are often unclear. Later phases predominating in high-frequency ranges (> 8 Hz) have been observed at stations paired with earthquakes in the forearc region of N Japan (Yomogida et al. 2002) and NE Japan (Hasemi and Horiuchi 2010), suggesting that delayed arrivals of the highfrequency waves reflect the spatial variations in seismic wave attenuation beneath these two regions, along with the presence of scatterers in the forearc region. However, the relationship between the later phases of the backarc intraslab earthquakes and the heterogeneous structures has not been revealed.
In this study, we investigate the high-frequency later phases (HFLPs) of backarc intraslab earthquakes in NE Japan. Wu and Irving (2018) have already reported arrival delays of high-frequency waves in an earthquake cluster at a depth of ~ 150 km, but here we further identified the HFLPs of earthquakes over the entire Pacific slab in the study area. We analyzed the amplitude records of intraslab earthquakes with HFLPs and revealed that the HFLPs were scattered waves generated by scatterers at depths corresponding to the lower crust and the mantle wedge in the forearc region; their propagation paths are similar to those of the forearc intraslab earthquakes retrieved by Hasemi and Horiuchi (2010).
The propagation paths of HFLPs through the forearc region resulted from the S-wave attenuation contrast across the volcanic front in NE Japan (e.g., Yoshimoto et al. 2006;Takahashi 2012). Because the amplitude of the direct S wave is strongly attenuated when passing across the volcanic front (e.g., Nakajima et al. 2013), the characteristics of the S-wave attenuation property in the backarc region in NE Japan are poorly understood. Therefore, the HFLPs of backarc intraslab earthquakes may provide a new tool for investigating heterogeneous structures in the subduction zones of N and NE Japan, including S-wave attenuation beneath the backarc region.

Characteristics of HFLPs
The direct S-wave signals of intraslab earthquakes at approximate depths of 100-250 km in NE Japan are commonly followed by an HFLP. The HFLPs are characterized by a predominance of frequencies higher than ~ 8 Hz, whereas direct S waves are rich in lower frequency signals (Figs. 2 and 3).
The precise arrival times of HFLPs are difficult to determine because the amplitudes of these signals gradually increase. However, the increasing amplitude of HFLPs is clearly seen in records, and the peak amplitude is discernible. As seen in Fig. 3, the peak amplitude shortly follows the arrival of direct S waves at stations in the forearc region, but is delayed by up to ~ 50 s in the backarc region. The HFLP peaks clearly propagate from the forearc region toward the backarc region with an apparent velocity of ~ 5 km/s. The characteristics of HFLPs can be summarized as follows: (1) a distinct high-frequency (8-32 Hz) signal following the direct S-wave arrival; (2) a remarkable peak amplitude, especially at 8-16 Hz, with long duration; (3) longer delay of peak amplitude in the backarc region than in the forearc region.  (Nakajima et al. 2009). The contour interval is 50 km Two plausible mechanisms can explain the observed delay of the peak amplitude of HFLP in the high-frequency range of intraslab earthquakes: multiple-forward scattering phenomena along the ray path of direct S waves (Obara and Sato 1995;Kennett 2005, 2008;Takahashi et al. 2009), and the presence of multi-passing waves (Yomogida et al. 2002;Hasemi and Horiuchi 2010). In their study of intraslab earthquakes in the forearc region, Hasemi and Horiuchi (2010) reported a later phase with similar characteristics to the HFLPs studied here [characteristics (1) and (2) above]. They attributed the later phase to scattered waves bypassing the high-attenuation mantle wedge. Like HFLPs, the later phases identified for forearc intraslab earthquakes were reported over the entire NE Japan (Hasemi and Horiuchi 2010). Therefore, we interpret the present HFLPs as scattered waves, as discussed in Hasemi and Horiuchi (2010).

Data
We selected stations and intraslab earthquakes located within 10 km of each of five cross-sections (A-E) almost perpendicular to the NE Japan (Fig. 1a). These intraslab earthquakes with clear HFLP signals were collected from the Japan Meteorological Agency (JMA) unified catalog over the 2003-2015 period. The magnitude range of the selected events was 2.5-4.5. Velocity seismograms were collected from the nationwide seismograph network operated by JMA, National Research Institute for Earth Science and Disaster Resilience (NIED) [Hi-net (Okada et al. 2004;NIED 2019)], and Tohoku University. The direct S-wave arrival times were manually picked.  Fig. 1b. a, c Seismograms bandpass-filtered at 1-32 Hz, and b, d their MS envelopes in different frequency ranges: 1-8 Hz (filled gray lines) and 8-32 Hz (thick black lines). Arrows indicate the identified peak amplitudes of the HFLPs. White inverted triangles show the expected delay times of the waves scattered at the passing points (their temporal ranges correspond to the location uncertainties displayed in Fig. 5). Gray diamonds indicate the expected delay times at the passing points in Fig. 7. Gray curves denote the arrival times of the direct P and S waves expected for the JMA2001 one-dimensional velocity model (Ueno et al. 2002). Dashed lines represent the locations of the arc volcanic front (VF) In this study, intraslab earthquakes at depths of about 50-100 km and 100-250 km were selected as the forearc and backarc events, respectively (Fig. 1b). Table 1 lists the numbers of backarc events and stations along each cross-section. Over cross-sections A, B, C, D, and E, we analyzed 50, 101, 49, 43, and 97 waveforms of backarc events, respectively.

Method
The propagation paths of the HFLPs were investigated using a source-scanning algorithm (SSA) (Kao and Shan 2004). The conventional SSA evaluates the location of an energy radiation source based on spatiotemporal changes in brightness values computed from the amplitude records of seismic waves. Using this SSA feature, we searched for passing points representing the distribution of scatterers explaining the observed peak delays of the HFLPs.
Following Hasemi and Horiuchi (2010), we assumed that each HFLP is a single scattered wave. The brightness, b ijk , of an earthquake-station pair is defined as where E ij is the mean-square (MS) envelope amplitude of the i th earthquake observed at the j th station, and τ ijk is the delay time of the scattered wave passed at the k th grid point from direct S-wave arrival. The delay time is given by τ ijk = τ L − τ S , where τ S and τ L are the travel times of the direct S wave and the scattered wave, respectively. The brightness, B k , at the k th grid point is calculated from the values of b all of the available earthquake-station pairs as follows: where b ijk is the brightness value normalized by the maximum b for each earthquake-station pair, and N is the number of available pairs. Normalizing the b value reduces the bias due to source and site-dependent effects in the SSA analysis (e.g., Kao and Shang 2004). The grid on which B k is maximized represents the common passing point that can explain the peak-amplitude delays of the HFLPs. ( Applying the SSA technique, Kosuga (2014) found the source location of the reflected waves observed for local inland earthquakes in NE Japan. In that study, the background coda-wave amplitude was subtracted from the observed records by the synthesized theoretical amplitude based on the single scattering model (Sato 1977;Sato et al. 2012). This procedure efficiently highlights the amplitude of the targeted later phases, but requires spatial homogeneity of the seismic wave attenuation. Therefore, it is inapplicable to the NE Japan subduction zone, in which attenuation properties are spatially variable (e.g., Tsumura et al. 2000;Takahashi 2012;Nakajima et al. 2013). Rather than Kosuga's approach, we considered the spatial changes of the velocity and attenuation structures in the brightness evaluations, as detailed below.
The present study considered the effects of both geometrical spreading and path attenuation of the scattered waves. In the single scattering model (Sato 1977;Sato et al. 2012), the background coda-wave energies are obtained by summing the wave energies of the singly scattered waves from all scatterers having a common delay time. The scatterers distributed on one isochrone of delay time form one scattering shell. As the scattering shell expands in time, the number of scatterers on the corresponding shell increases and the amplitude of the scattered wave from each scatterer is reduced by geometrical spreading. Thus, the contributions of the individual scatterers on the shell to the background codawave amplitude should decrease over time. To express the expansion effect of the scattering shell relative to that of a direct S wave, we define the following weight, w T : We also define the path attenuation effect, w Q , on each scattered wave relative to that of the direct S wave: where t * S and t * L denotes the path attenuations of the direct S wave and scattered wave, respectively. The path attenuation, t * , is obtained by integrating the travel time, τ , and the quality factor, Q , of the attenuation along a given ray path: where V is the seismic velocity. Multiplying w T and w Q by the MS amplitude, the brightness of an earthquakestation pair given by Eq. (1) is modified to Using Eq.
(2) together with Eq. (6) in the SSA, we estimated the locations of passing points of HFLPs in the attenuating medium. The inclusion of w Q often yielded unrealistic estimates of B k , because the attenuation structure was oversimplified as described below. This ill-posed behavior was improved by limiting w Q to the 0.1-10.0 range.
Characteristic structures, such as the subducting Pacific slab and upwelling mantle flow, are commonly known to be distributed along NE Japan (e.g., Nakajima et al. 2001Nakajima et al. , 2013Takahashi 2012). Therefore, the present SSA analysis was conducted on two-dimensional vertical cross-sections across the NE Japan arc, labeled as cross-sections A-E in Fig. 1. The model space was divided into four layers by three structural boundaries: the Conrad and Moho discontinuities (Katsumata 2010), and the upper boundary of the Pacific slab (Nakajima et al. 2009). The S-wave velocity in each layer was taken from recent tomographic models (e.g., Nakajima et al. 2001;Matsubara et al. 2017). We also incorporated the low-velocity subducting oceanic crust in the uppermost part of the Pacific slab (e.g., Tsuji et al. 2008;Shiina et al. 2013). In NE Japan, a high S-wave attenuation ( Q S ) zone in the mantle wedge has been constrained to the depth of about 60 km (e.g., Takahashi et al. 2009;Takahashi 2012). This high-attenuation zone has also been detected as a P-wave attenuation ( Q P ), and is continuously distributed into a deeper position of the backarc mantle wedge (e.g., Tsumura et al. 2000;Nakajima et al. 2013). Based on the results of P waves and assuming that Q P /Q S = 2 (e.g., Shito and Shibutani 2003;Nakajima et al. 2013), we adopted the high-attenuation mantle wedge extending below the backarc region. The assumed attenuation structure reproduces the path-averaged properties of direct S waves, as confirmed in the Discussion Section. For example, the velocity and attenuation structures over cross-section D are shown in Fig. 4.
The travel times and path attenuations of the direct S and scattered waves through the assumed velocity and attenuation structures were computed using the ray-tracing method of Zhao et al. (1992). Grid points were set at 2.5-km intervals in the horizontal and vertical directions.
The brightness values in the frequency range of 8-16 Hz were determined following Hasemi and Horiuchi (2010), who reported delayed high-frequency waves of forearc intraslab earthquakes in NE Japan. From the two horizontal components of the seismograms, the MS amplitudes were synthesized with a 1-s moving average. The noise amplitude was calculated as the 5-s average of the MS amplitudes before the P-wave arrival. The MS envelopes were visually checked, and then accepted if their amplitude at 2.5 times the lapse time from the direct S-wave travel time was twice the noise amplitude. Figure 5 shows the estimated locations of the passing points together with the brightness values B . The passing points in the five cross-sections were estimated at horizontal positions of 35-55 km and depths of 30-80 km. The passing points in the forearc region reproduced the peak-amplitude delays of the HFLPs, and the increasing delay times towards the backarc region (see Fig. 3b, d).

Results
The brightness values depended on the assuming velocity and attenuation models. Representative results in cross-sections D and E are shown in Fig. 6, while the passing points in the other cross-sections were determined at almost identical in Fig. 5. When the attenuation effects were ignored, the passing points appeared in the backarc mantle wedge and the forearc region (see Fig. 6b, f, respectively). The passing point in cross-section D induced no long delay times at the backarc stations, contradicting the observed HFLP features (Fig. 3b). These examinations of different models firmly support the robustness of the passing point locations in the forearc region and indicate that scatterers in the backarc region a b Fig. 4 a S-wave velocity and b attenuation models including the high-attenuation zone in the mantle wedge (HMW) along cross-section D (Fig. 1). b Values in white boxes are the Q S values at 8-16 Hz. The three black lines represent the structural boundaries of the Conrad and Moho discontinuities (Katusmata 2010), and the upper boundary of the Pacific slab (Nakajima et al. 2009) do not affect the observed peak amplitudes in the records with HFLPs even if they exist there.
Because the spatial distributions of the brightness values do not necessarily represent the location uncertainty in the passing points and the intensity of scatterers, the passing-point uncertainties were determined in bootstrapping tests. In these tests, earthquake-station pairs were randomly selected, and w T and w Q were perturbed by randomly generating additional values of τ L and t * L , respectively. The additional values followed Gaussian distributions with zero average. The standard deviations of τ L and t * L were assumed as 1.0 s and 0.4 s, respectively, based on the residuals of the travel times and path attenuations in the reference structures (e.g., Nakajima et al. 2001Nakajima et al. , 2013Matsubara et al. 2017). After 10,000 trials, we obtained the 95% confidence interval of the passing point location.
The 95% confidence intervals were mainly distributed in the forearc region (Fig. 5). According to the results, the scattered waves related to HFLPs propagated outside of the high-attenuation zones in the backarc mantle wedge, as discussed in previous studies (Hasemi and Horiuchi 2010;Wu and Irving 2018). Note that our dataset may not precisely constrain the depths of the passing points, because the confidence intervals were vertically elongated. In some cross-sections, the confidence intervals were split over two areas at horizontal positions of around 0-20 km. Therefore, although the scatterers along each cross-section were difficult to determine exactly, the passing points definitely lay on the forearc side in all five cross-sections.

Discussion
We demonstrated that a region of the passing points in the forearc region can explain the observed peak delays of HFLPs. Hasemi and Horiuchi (2010), who similarly reported scattering regions in NE Japan, located the scattering sources in the lower crust based on the arrival times of high-frequency waves of forearc intraslab earthquakes. However, even after accounting for the location uncertainties, the passing points estimated from the a b c d e Fig. 5 Estimated locations of passing points from the backarc intraslab earthquakes. White crosses and black contours indicate the passing points defined by the maximum brightness values, B , and the 95% confidence intervals, respectively. Stars and inverted triangles denote the earthquakes and stations, respectively. Gray-shaded areas represent the high-attenuation zone in the mantle wedge. Black lines depict the structural boundaries shown in Fig. 4 a b e f c g d h Fig. 6 Estimated locations of passing points at cross-sections D and E in four models: a, e the HH2010 model, a layered model with velocities used in Hasemi and Horiuchi (2010). b, f The WoQ model, excluding attenuation effects. c, g The HighQ model: assuming weak attenuation ( Q S = 300 ) in the backarc mantle wedge (gray-shaded area). d, h The LowQ model, assuming strong attenuation ( Q S = 100 ) in the backarc mantle wedge. The attenuation structure in the HH2010 model and the velocity structure in the WoQ, HighQ, and LowQ models are the same those shown in Fig. 4. Other symbols are explained in the caption of Fig. 5 backarc events in the present study appear to be deeper than the lower crust in cross-sections D and E (Fig. 5d, e).
To further discuss the possible depths at which scattered waves are excited, we estimated the scattering regions of forearc intraslab earthquakes, in the same manner estimating passing points of the HFLPs of the backarc intraslab earthquakes. The results are shown in Fig. 7, and the numbers of forearc events in each crosssection are listed in Table 1. Most of the passing points related to high-frequency waves of the forearc events were located within or just below the lower crust (Fig. 7), as previously reported by Hasemi and Horiuchi (2010). In cross-section B, the estimated depth of passing point was ~ 90 km. This result may be anomalous, because there is large vertical uncertainty in the lower-crust mapping. Figure 8 compares the expected delay times of the scattered waves at the passing points in Figs. 5 and 7. Although the passing-point depths in cross-section D differ by ~ 40 km between Figs. 5 and 7, the delay times of the scattered waves were similar in both scenarios (Fig. 8). In cross-section E, the expected delay times at the passing point were faster in Fig. 5 than in Fig. 7 (Fig. 8). However, when the location uncertainties in the passing points are translated into ranges of delay times, the delay   (Fig. 7). Gray-filled symbols are the expected delay times in the cross-sections A-C. White symbols are the expected delay times in cross-sections D and E, and their temporal ranges corresponding to the 95% confidence intervals in the forearc region in Figs. 5 and 7 times at the passing points are not significantly different in Figs. 5 and 7. That is, the obtained differences in passing-point depths produce no significant differences in the delay times of the scattered waves in the backarc events. Therefore, we consider that scattering sources in the lower crust (Hasemi and Horiuchi 2010) can explain the generation and propagation of the HFLPs. Hasemi and Horiuchi (2010) also suggested the highintensity scattering sources, with a reflection coefficient of 0.56. Such a high reflection coefficient was necessary to reproduce the relative amplitude of later phases following the direct S waves in forearc events. If abnormally intense scattering sources were presented, seismic records of the forearc stations would become complicated by excited reflected and converted waves. However, signals other than direct waves were not recognized in the seismograms recorded at the forearc stations (Fig. 3). Hasemi and Horiuchi (2010) assumed a homogeneous attenuation structure when deriving the high-intensity scattering sources, whereas the attenuation properties in NE Japan are spatially variable (e.g., Yoshimoto et al. 2006;Takahashi 2012;Nakajima et al. 2013). Thus, when analyzing the scattering properties affecting HFLPs, one must consider the heterogeneity of the attenuation structure.
The S-wave attenuation properties in NE Japan have been constrained down to the depth of about 60 km (e.g., Takahashi et al. 2009;Takahashi 2012). At depths exceeding 60 km, the contrasting attenuation properties between the forearc and backarc mantle wedge affect the seismic wave propagations of backarc events. To confirm the S-wave attenuation contrast in the mantle wedge during backarc events, we applied the coda-normalization method (e.g., Aki 1980;Yoshimoto et al. 1993) to seismic records containing HFLPs. The scattering properties affecting coda-wave amplitudes can be regarded as almost homogeneous when the lapse time is sufficiently spaced from the direct S-wave arrival. Therefore, the Q S values obtained by the coda-normalization method probably characterize the spatial variations in attenuation property along the propagation paths of direct S waves across the volcanic front.
The direct S-and coda-wave amplitudes were averaged over 2.56-s and 5.12-s time windows, respectively, starting from the S-wave arrival and 100-s lapse time. The envelopes were selected by the criterion described in the Method section. In addition, an envelope was eliminated if the signal-to-noise ratio of its coda-wave amplitude was less than 3 or if the travel time of the direct S wave was longer than 50 s. The Q S was then estimated by the method following Yoshimoto et al. (1993).
During Event A, the lapse times of the direct S waves and HFLPs in the backarc events were 30-40 s and 40-60 s, respectively (Fig. 3b). At such earthquake-station pairs, the direct S waves cross in the backarc region, whereas the HFLPs propagate through the forearc region, as revealed in this study. Thus, the Q −1 S values obtained in the backarc and forearc regions likely represent the path-averaged attenuation properties of direct S waves and HFLPs, respectively. Assuming Q −1 S values of ~ 0.002 and ~ 0.0008 for direct S waves and HFLPs, respectively, a reflection coefficient of about 0.05-0.10 is required to explain the relative amplitudes of the high-frequency phases discussed by Hasemi and Horiuchi (2010). That is, the relative amplitudes of the HFLPs can be reproduced by scattering sources without an anomalously high intensity. This magnitude of intensity is a value when considering an isolated scattering source. If multiple scattering sources are present, the intensities of the individual scatterers will be further reduced. The duration of the HFLPs often exceeded the delay times of about 10-30 s that correspond to the spatial uncertainties in the passing points (Fig. 3b, d). The location uncertainties were 30-80 km in the horizontal direction and 50-100 km in the depths (Fig. 5). Thus, scatterers other than those discussed as the passing points are broadly distributed in the forearc region, and can attribute to the observed long duration of HFLPs.
Finally, we considered the across-arc variations in the attenuation properties of NE Japan. The Q −1 S values determined at the forearc stations (Fig. 9a) are comparable to previously estimated Q S values of forearc events (e.g., Hoshiba 1993;Takahashi et al. 2005;Takahashi 2012). This consistency suggests weak S-wave attenuation in the Pacific slab, similar to that discussed for P waves (e.g., Shiina et al. 2018). In other words, the low-Q −1 S values at the forearc stations of NE Japan support the localization of high-attenuation zones in the backarc mantle wedge (e.g., Tsumura et al. 2000;Takahashi 2012;Nakajima et al. 2013;Liu et al. 2014). When determining the Q −1 S values, we included the MS envelopes even if then direct S waves were apparently hidden by P-coda waves. Additionally, coda-wave amplitudes decay more rapidly in the backarc region than in the forearc region (Yoshimoto et al. 2006). These factors would underestimate the Q −1 S values in the backarc region, suggesting that the Q −1 S contrast in the forearc and backarc regions of NE Japan may differ by a greater extent than estimated here.
We conclude that two fundamental factors affect the observations of HFLPs of backarc intraslab earthquakes in NE Japan: the presence of scatterers in the forearc region and strong attenuation contrast across the volcanic front in the mantle wedge. In particular, the highattenuation zone seems to highlight the HFLPs in seismic records by attenuating the amplitudes of the direct S wave. Unfortunately, the actual mechanisms of the scattering phenomena in the forearc region, including scattering intensities and their spatial variations, that affect the HFLP observations were not retrieved in this study. Takemura and Yoshimoto (2014) discussed that smallscale heterogeneity localized around the lower crust and the uppermost part of the Pacific slab causes peak delay and envelop broadening in high-frequency ranges. Additionally, other studies have proposed scattering sources in the overlying plate (e.g., Yomogida et al. 2002;Katayama et al. 2012;Zhao et al. 2015) and/or within the subducting plate (e.g., Furumura and Kennett 2005). Numerical approaches (e.g., Furumura and Kennett 2005;Przybilla et al. 2009;Takemura et al. 2015a, b) have quantitatively constrained the high-frequency wave propagation and its relations to heterogeneous structures. By numerically simulating HFLP propagation, we expect to reveal the detailed properties of high-frequency wave scattering in the NE Japan subduction zone.

Conclusions
Distinct later phases with rich high-frequency (> 8 Hz) components were observed in the seismograms of intermediate-depth intraslab earthquakes in the backarc region of NE Japan. Applying the SSA, we estimated the propagation paths of the HFLPs from their amplitude records. The passing points were estimated to lie between the lower crust and the uppermost part of the Pacific slab in the forearc region, indicating that HFLPs are scattered waves generated in the forearc region (Fig. 10). The passing-point depths ranged between 30 and 80 km, but the uncertainties in depths were enlarged by the limited