Seismic detection of the X‐discontinuity beneath the Ryukyu subduction zone from the SdP conversion phase

The X‐discontinuity, which appears at the depth of approximately 300 km, is an important seismic interface with positive velocity contrasts in the upper mantle. Detecting its presence and topography can be useful to understand phase transformations of relevant mantle minerals under the high‐temperature and high‐pressure circumstance of the Earth's interior. In this study, we detect the X‐discontinuity beneath the Ryukyu subduction zone using five intermediate‐depth events recorded by the dense Alaska Regional Network (AK). The X‐discontinuity is successfully revealed from the robust slant stacking of the secondary down‐going and converting SdP phases. From the depth distribution of conversion points, we find that the X‐discontinuity's depth ranges between 269 km and 313 km, with an average depth of 295 km. All the conversion points are located beneath the down‐dipping side of the Philippine Sea slab. From energy comparisons in vespagrams for observed and synthetic seismograms, the strong converted energy is more likely from a thin high‐velocity layer, and the S‐wave velocity jumps across the X‐discontinuity are up to 5% to 8% with an average of 6.0%. According to previous petrological and seismological studies, the X‐discontinuity we detected can be interpreted as the phase transformation of coesite to stishovite in eclogitic materials within the oceanic crust.


Introduction
The X-discontinuity in the upper mantle has been described from previous seismological observations as a seismic interface with a velocity contrast of 1.5%-4.8% for P-and S-waves and a gradient zone less than 5 km (e.g., Revenaugh and Jordan, 1991;Deuss and Woodhouse, 2002;Bagley and Revenaugh, 2008;Shen et al., 2014;Chen T et al., 2015). The X-discontinuity is found at depths ranging from 250 to 350 km; it is sometimes termed "the 300 km discontinuity" (e.g., Zhang and Lay, 1993;Williams and Revenaug, 2005;Xie CX et al., 2012;Schmerr et al., 2013). The characteristics of the X-discontinuity are sporadically observed in some regions, such as the Eurasia, Australia and North America continents (e.g., Hales et al., 1980;Li AB et al., 2002;Williams and Revenaug, 2005), the subduction zones and beneath oceanic plates (e.g., Revenaugh and Jordan, 1991;Zhang Z and Lay, 1993;Deuss and Woodhouse, 2002;Bagley and Revenaugh, 2008), and a few hot spots (Courtier et al., 2007;Bagley et al., 2009). Shen XZ et al. (2014) proposed that the X-discontinuity could be ubiquitous, using global stacks of P-wave receiver functions. However, it re-mains uncertain whether the X-discontinuity is a global feature.
The mineral physics mechanism for the X-discontinuity is also debatable (e.g., Williams and Revenaug, 2005). So far, several mechanisms have been proposed to interpret it, such as the formation of hydrous phase A (Mg 7 Si 2 O 8 (OH) 6 ) in subduction zones (Revenaugh and Jordan, 1991;Ohtani, 2005), the transition of orthoenstatite to high-pressure clinoenstatite (e.g., Woodland, 1998;Xie CX et al., 2012), and the transformation from coesite to stishovite in subducted eclogitic materials (e.g., Williams and Revenaug, 2005;Chen T et al., 2015). The phase A is limited in cold subducted slabs, and not stable above 1000 °C at 300 km depth (Komabayashi et al., 2005;Ohtani, 2005;Chen T et al., 2015). The impedance contrast caused by the transition of orthoenstatite to high-pressure clinoenstatite in a pyrolite composition is too small (1.0%-1.4% for S-wave and 0.7%-0.9% for P-wave) (e.g., Williams and Revenaug, 2005;Chen T et al., 2015). The phase transformation from coesite to stishovite can produce large P-and S-wave impedance increases (3%-8%) proportional to the content of free silica (e.g., Chen T et al., 2015), which could come from the eclogite assemblage. (e.g., Woodland, 1998;Jin ZM et al., 2001). The eclogitic materials in the upper mantle are likely remnant oceanic crust, consisting primarily of subducted mid-ocean ridge basalt (MORB) ( e.g., Williams and Revenaug, 2005;Schmerr et al., 2013;Woodland et al., 2014). To some extent, the coesite to stishovite transformation can be thought of as a probe for subducted oceanic materials in the upper mantle (e.g., Chen T et al., 2015) The Ryukyu subduction zone is located in the western margin of the Philippine Sea ( Figure 1a). The Philippine Sea Plate has been subducted beneath the Eurasian Plate along the Ryukyu Trench (e.g., Hall et al., 1995;Zang SX et al., 2002;Fukao and Obayashi, 2013;Huang ZC et al., 2013), since approximately 50 Ma, when the Pacific Plate changed its movement direction from NNW to NWW (Park et al., 1998;Sharp and Clague, 2006). The present convergence rate between the Philippine Sea Plate and Eurasian Plate is as much as 5-7 cm/year (Park et al., 1998). Seismic tomographic studies revealed the Philippine Sea Plate as the relatively high-velocity zone, with the thickness of 30-40 km (e.g., Zhao DP et al., 2002;Huang and Zhao, 2006;Li C et al., 2008a). The intermediatedepth earthquakes occur at depths down to about 280 km. Abundant earthquakes in the Ryukyu area provide reliable and high-quality waveform data from global seismic networks. Seismic detections of mantle discontinuities in this area are useful to understand the phase transformations of relevant minerals under the circumstances of the Earth's interior.
In this study, we use data from the dense permanent Alaska Regional Network (hereafter, AK) to detect the X-discontinuity beneath the Ryukyu subduction zone (Figure 1c). The slant stack method is applied to extract the secondary SdP phase, which is the conversion of a down-going S wave to a down-going P wave at a discontinuity with the depth "d" (Figure 1a). Our detections imply possible phase transformations in the oceanic crustal fragments, and are helpful for understanding the regional geodynamics about slab subduction activities.

Data Collection
We collected data from the intermediate-depth earthquakes that occurred from 2005 to 2017 beneath the Ryukyu subduction zone. Because our goal was to detect the X-discontinuity around 300 km, we chose earthquakes with focal depth range between 130 and 220 km, in order to avoid the interference of direct P and pP phases on the weak SdP. Furthermore, we chose only earthquakes with simple source time series and source duration less than 5.0 s, because complicated source time functions may cause a misreading of the multiple phase (Zhou YZ et al., 2012); event magnitudes should be strong enough that observed seismograms would have clear P and pP phases. Finally, five intermediate-depth events were selected (Figure 1b). The basic source parameters of Ev.1-5 (Table 1) Figure 1. The source to receiver geometry along the great circle path (a), epicenters of five events used in this study (b) and station locations of the Alaska Regional Network (AK) (c). The epicentral distance contours are shown in dark gray dashed lines in (a). The insert map in (a) describes the schematic illustration of the secondary SdP conversion phases. The SdP phase means that a down-going S wave converts to a down-going P wave at the discontinuity with the depth labeled as "d".
In this study, we use data from the dense AK network, which has suitable epicentral distances (62°-75°). The AK network has been operated by the Alaska Earthquake Center since 1987, and consists of 126 broadband stations. To ensure the stability of slant stacking, we select only 94 stations with higher density and smaller distance ranges-less than 7° (Figure 1c). On the basis of the grid search method, we choose the station with the minimum differential distance to other stations as the central station.
After retrieving the vertical component seismograms of five events from the Data Management Center of Incorporated Research Institutions for Seismology (IRIS DMC, http://ds.iris.edu/), we manually pick seismograms with clear P phases, and the signal-to-noise ratios (SNR) more than 3.0. For each event-to-network pair, 60 seismograms with top SNR Cui et al., 2017) are used, if sufficient trace numbers are available. In the record sections for five event-to-network pairs, the SdP phases are transparent for Ev.3 (Figure 3), and visually weak for other events ( Figure S1).

Data Processing
The SdP conversion phases are usually weak in observed seismograms, and subject to possible background noise interference. To suppress the noise level and enhance consistent seismic signals, a stacking method is commonly applied in the vertical-component network/array data processing (e.g., Vidale and Benz, 1992;Kawakatsu and Niu, 1994;Schimmel and Paulssen, 1997;Rost and Thomas, 2002;Tibi and Wiens, 2005;Niu FL, 2014). To extract the secondary SdP phase effectively, the N-th root slant stack method (e.g., Kawakatsu and Niu FL, 1994;Vinnik et al., 1998;Rost and Thomas, 2009;Cui QH et al., 2018) is used to stack seismograms in the domain of the travel-time and slowness differences. To make the results more convincing, we use both N=4 and N=1 (also called linear stack) to jointly identify the SdP phases in vespagrams and constrain their relative slowness.
The procedure of seismic data processing can be described as follows: (a) Removing the means, trends, and instrument responses from raw seismograms to velocity records; (b) Filtering seismograms at 0.2-1.0 Hz (e.g., Castle and Creager, 2000;Zang SX et al., 2006) with a 2-pole Butterworth filter ( Figure  2), then aligning and normalizing them with reference P phases; (c) Stacking observed seismograms by the 4-th root and linear slant stack methods, then using Hilbert transform to get their envelopes, and finally taking the logarithm to get the vespagrams; (d) Calculating the travel-time and slowness of P and SdP phases with the IASP91 model (Kennett and Engdahl, 1991) using the TauP package (Crotwell et al., 1999), then identifying the SdP phases and inverting their depths; (e) Calculating synthetic seismograms using the modified IASP91 models adding a thin layer (thickness=7 km) under the conversion depths in (d), and then linearly stacking them to the vespagrams; (f) Determining the S-wave velocity jumps from the energy comparisons between vespagrams in (c) and (e), and then inverting their horizontal locations with the modified IASP91 models.
Based on an analysis of the time-frequency spectrums and instrument responses, we filter the original seismograms with the frequency band of 0.2-1.0 Hz, which is sensitive to mantle discontinuities and is used in the slant stacking (e.g., Castle and Creager, 2000;Zang SX et al., 2006;Chen J et al., 2014). As shown in Figure 2, background noise in the original velocity seismograms is effectively suppressed after the band-pass filtering; the P and SdP phases become clearer in the filtered seismograms.
To better extract the SdP conversion phases, we align the observed seismograms with P-wave peaks, which denote the maximum energy of P phases. Taking the record section of Ev.3 as an example (Figure 3a), all the traces have maximum amplitudes at 0 s. The travel-time differences between the P-wave peaks and predicted travel-times concentrate at 2.3-4.1 s for the record section.
An obvious secondary phase (S269P) becomes visible at ~16.0 s and has positive polarity with the P phase ( Figure 3a). If the SdP conversion phase occurs at a sharp discontinuity, then an S-wave velocity jump of 15% is required to generate the energy of S269P (-18 dB) from synthetic seismograms ( Figure S2) calculated with the reflectivity method (Wang RJ, 1999). Such a large velocity jump is unusual in the upper mantle, and the strong energy may imply conversion through the scatterers, rather than a sharp discontinuity. For computation efficiency, in this study we use a 1-D thin layer model to approximately constrain the parameters of scatters. Based on the IASP91 model, we adopt a high-velocity layer with the thickness of 7 km, proposed in Niu FL (2014) and Yang ZT and He XB (2015), and calculate the synthetic seismograms ( Figure 3b). The S-wave velocity jump of 8% yields energy comparable to that of the observed one in vespagrams ( Figure S2). From previous studies on the X-discontinuity (e.g., Zhang Z and Lay, 1993;Williams and Revenaug, 2005;Shen XZ et al., 2014;Chen T et al., 2015), such a velocity change (ΔV s /V s =8%) strongly suggests the coesite to stishovite transformation in a thin layer.

Stacking Results
We apply the 4-th root and linear slant stackings to analyze five event-to-network pairs, and then obtain the corresponding vespagrams ( Figure 4). In each sub-plot, the P phase is located at travel-time 0 s and has maximum beam power (0 dB); purple lines show the theoretical arrivals with travel-time and differential slowness for SdP, and provide the essential reference for conversion phase identification. Only the SdP phase, which is near the theoretical trends with slowness deviation within -0.1 to +0.1 s·deg -1 , can be used for studying mantle discontinuities. We identify those secondary phases in vespagrams, both from 4-th root and from linear slant stackings, as the SdP phases and mark them with the dark gray arrows in Figure 4a-e.
The S300P phases are all visible for five event-to-network pairs. The secondary phases in Figure 4b show the complexities around the theoretical SdP at 15-25 s, and we regard the stronger phase as the S300P phase among the possible SdP phases in the vespagrams. After the time-to-depth conversion, the conversion depths of the identified SdP phases are estimated at 269 to 313 km with the average depth of 295 km. In Figure 4e, there is a weak secondary phase at ~22 s labeled as S412P, which is in agreement with the vespagram for synthetic seismograms ( Figure S3). The S412P in Figure 4e corresponds to the 410 km discontinuity, which is inferred to be the phase transformation from olivine to wadsleyite (α→β) (e.g., Revenaugh and Jordan, 1991;Ringwood, 1991;Deuss and Woodhouse, 2002). The theoretical S660P phases for five event-to-networks pairs are located near the strong pP or between the pP and sP phases (Figure 4), and they would be hard to be distinguished.
In this study, we mainly focus on the S300P phase and the strength of the X-discontinuity. To effectively constrain the slowness of the S300P phase, we manually pick the observed slowness in vespagrams and list them in Table 1. The differential slowness to the theoretical one is 0.09 s·deg -1 for Ev.1 and Ev.2, and 0.00-0.04 s·deg -1 for Ev.3-5 (Table 1). Those slowness differences imply ray path deviations of SdP phases relative to the great circle plane of P phases. We applied the beam-forming analysis (also called f-k analysis (e.g., Rost and Thomas, 2002) in a time window around P (-2 to +2 s) and S292P (-1 to +1 s) for Ev.1 as an example. The relative slowness of S292P is -0.10 s·deg -1 (Figure 5), which is consistent with the observed one (-0.11 s·deg -1 ) in the vespagrams (Figure 4a). The relative back azimuth of S292P is 0.14° to P phase ( Figure 5), so that the secondary phase is slightly out of the great circle plane of P wave. The ray path and slowness differences between SdP and P phases may be caused by the pos- sible slightly dipping geometry of scatterers beneath the subduction zone (e.g., Niu FL, 2014).
We attempt to constrain the velocity contrasts across the X-discontinuity: First, for linear stacking we manually pick the maximum beam power for S300P phases in the vespagrams ( Figure 4). Second, we calculate the synthetic seismograms with gCMT focal mechanisms (Table S1) and S-/P-wave velocity and density contrasts. The slope relationships between velocity and density perturbations are taken from Chen T et al. (2015). We linearly stack the synthetic seismograms into vespagrams. Third, we make the energy comparisons between both vespagrams. The Swave velocity jumps across the X-discontinuity for Ev.1, 3 and 5 are obtained ( Figure S2 and S3); they range from 5% to 8% with an average of 6%; the jumps for Ev.2 and 4 are not well constrained, due to the extremely weak SdP in the synthetic seismograms.

Depth Uncertainty Estimation
We estimate that the depth uncertainty of conversion points is about 11 km, based on these three considerations: (a) Phase picking errors. The arrival time picking errors are usually less than 0.5 s (Collier and Helffrich, 1997), which introduces up to 5 km depth error, calculated with the IASP91 model (Kennett and Engdahl, 1991); (b) Mantle heterogeneities. We extract the depth sections from the global S-wave S40RTS model (Ritsema et al., 2011) and P-wave GAP_P4 model . The tomographic depth sections along ray paths from event to network are shown in Figure 6. We use the fast marching method (de Kool et al., 2006) to calculate the travel-times of P and SdP phases across tomographic sections. We find that mantle heterogeneities contribute up to 1 km depth error and have little influence on the time-todepth conversion; (c) Focal depth error. This error is generally considered to be 5 km (e.g., Collier and Helffrich, 1997;Chen J et al., 2014).

Discussions
As projected on the high-resolution bathymetry of the GTOPO30 model (https://lta.cr.usgs.gov/GTOPO30) (Figure 7a), the conversion points concentrate in the northern and southern parts of the Ryukyu subduction zone. The conversion depths range between 269 and 313 km with an average of 291 km in the northern part, and between 292 and 303 km with an average of 298 km in the southern part. From our results, the average depth of the X-discontinuity is 7 km deeper in the southern part than that in the northern part.
We calculate the distances between the conversion points, events, and the Ryukyu Trench, then project them in the vertical depth  (Figure 7b). The seismicity within the range in Figure 7a is obtained from ISC-EHB Bulletins (1960-2014) (http://www.isc.ac. uk/isc-ehb/), and also projected in the depth section. The approximate geometry of the subducted Philippine Sea slab is depicted with dashed lines in Figure 7b. All the conversion points are located beneath the down-dipping side of the Philippine Sea slab, which has been continuously subducting beneath the Eurasian Plate from the Ryukyu Trench (e.g., Park et al., 1998;Zang et al., 2002;Huang JL and Zhao DP, 2006).
Previous studies revealed that the X-discontinuity (250-350 km) has a velocity contrast of 1.5%-4.8% (P-and S-waves) and a small gradient zone of less than 5 km (e.g., Revenaugh and Jordan, 1991;Deuss and Woodhouse, 2002;Bagley and Revenaugh, 2008;Shen XZ et al., 2014;Chen T et al., 2015). Revenaugh and Jordan (1991) reported an X-discontinuity near 300 km depth beneath the northwestern Pacific with reflectivity coefficients of 1.1%-2.9%. Zhang Z and Lay (1993) found that the X-discontinuity beneath the western Pacific regions (Izu-Bonin, North Korea and North Okhotsk Sea) is at 325 to 335 km, with average S-wave impedance contrast of 9%, based on stacking of sSH precursors. Bagley and Revenaugh (2008) reported X-discontinuity near 300 km depth beneath the Pacific Ocean, with an S-wave impedance increase of 3%-8%, based on ScS reverberations. Xie CX et al. (2012) found that the X-discontinuity beneath the Tonga subduction zone is at 301 to 328 km, with an average depth of 314 km. The depth of the X-discontinuity we detected ranges between 269 km and 313 km with an average depth of 295 km, and the variation (44 km) is a little larger than that in previous Pacific Ocean results (e.g., Revenaugh and Jordan, 1991;Zhang Z and Lay, 1993;Bagley and Revenaugh, 2008;Xie CX et al., 2012). From the energy comparisons in the vespagrams for observed and synthetic seismograms, the S-wave velocity jumps across the X-discontinuity are up to 5% to 8% with an average of 6%, and the cor-   Figure S2. The vespagrams for synthetic seismograms of Ev.3 calculated with a layer model (thickness=7 km, ΔV s /V s =8%) (a) and a sharp discontinuity (V S /V S =15%) (b). The upper panels show the vespagrams derived from the linear slant stacking. Black and green circles denote the S300P and S410P predicted by the IASP91 model ( Kennett and Engdahl, 1991). Purple lines denote the trends of travel-time and slowness differences of SdP to P phase. The lower panels show the synthetic seismograms.  Figure S3. The vespagrams for synthetic seismograms of Ev.1 (a) and Ev.5 (b). Black, green and blue circles denote the S300P, S410P and S660P phases predicted by IASP91 model (Kennett and Engdahl, 1991). Purple lines denote the trends of travel-time and slowness differences of SdP to P phase. The lower panels show the synthetic seismograms.
responding S-wave impedance contrasts are up to 6.6% to 10.6% with an average of 8.0%. These depth and velocity contrasts seem to be consistent with previous X-discontinuity results in other regions (e.g., Zhang Z and Lay, 1993;Bagley and Revenaugh, 2008). Such a velocity or impedance contrast can be well explained by the coesite to stishovite phase transformations with free silica, which could come from the eclogite assemblage formed by the subducted oceanic crust (e.g., Jin ZM et al., 2001;Williams and Revenaug, 2005;Schmerr et al., 2013;Woodland et al., 2014).
As shown in Figure 7b, the conversion points are close to the subducted Philippine Sea Plate. On the basis of the waveform synthetics, the strong scattered energy is likely converted from a thin high-velocity layer. The thin layer at this depth range may be fragments of oceanic crust, where transparent coesite phase transformation in the eclogitic materials can occur (e.g., Woodland et al., 2014).
Numerical modeling experiments revealed that slab detachment can take place over a few million years in the subduction zone, and it will be accelerated by non-Newtonian strain-rate softening and focused thermal erosion (e.g., Gerya et al., 2004;Wortel and Spakman, 2000); the detached slabs would sink into deep mantle with some coherent rotations, for their density is higher than that of the ambient mantle (Gerya et al., 2004). Therefore, we can infer that the oceanic crustal fragments we detected possibly detached from the Philippine Sea slab in its earlier subduction period, after that the Philippine Sea slab continued subducting and moving forward to its present location. The detached crustal material has a higher density than the ambient mantle after the phase transformation, and would thus be sinking deeper.
The free silica content plays an importance role in property changes, such as the S-wave, P-wave impedances and density contrasts, in the coesite to stishovite phase transformation (e.g., Ricard et al., 2005;Williams and Revenaug, 2005;Woodland et al., 2014). Richard et al. (2005) proposed that the S-and P-wave velocity contrasts are about 7.3% and 5.8% respectively, for the MORB mineralogy with 14 wt% silica. Chen T et al. (2015) suggested that 4-8 wt% SiO 2 is required to reproduce the seismically observed Pand S-wave impedance contrasts of 3%-8% and density increase of 2.0%. Taking the linear dependence slope (d(lnV s )/d(SiO 2 )-0.60 (wt%) -1 ) between S-wave perturbation and SiO 2 content proposed in Chen T et al. (2015), the average S-wave velocity jump of 6% corresponds to the 10 wt% SiO 2 in the MORB.

Conclusion
In this paper, we study the X-discontinuity beneath the Ryukyu subduction zone using five intermediate-depth events recorded by the dense AK network. The presence and depths of the X-discontinuity are clearly revealed by the SdP conversion phases, through the 4-th root and linear slant stack methods. The robust SdP conversion points appear from 269 to 313 km with an average depth of 295 km. The conversion points are all located beneath the down-dipping side of the Philippine Sea slab. From energy comparisons in the vespagrams for observed and synthetic seismograms, the S-wave velocity jumps across the X-discontinuity are up to 5% to 8% with an average of 6%. Based on the wave-form synthetics, we infer that the strong converted energy is likely from a thin high-velocity layer. According to previous studies, the X-discontinuity we detected can be interpreted as the phase transformation of coesite to stishovite in eclogitic materials within the oceanic crust.