Seismic structure of subducted oceanic crust near the slow-earthquake source region in the southern Ryukyu arc

Abstract

Seismic tomography and receiver function analysis were carried out to investigate the relation between the slab structure in the southern Ryukyu arc region and the occurrence of slow-slip events that repeat biannually. For calculation of the receiver function, 212 teleseismic earthquakes with magnitudes larger than 6.0 were selected, and teleseismic waveforms were observed using two short-period seismometers and one broadband seismometer. Assuming that each later phase in a receiver function was a wave converted from P to S at depth, the time-domain receiver function was transformed into depth domain along each ray path using a reference velocity model. P- and S-wave arrival times, selected manually by the Japan Meteorological Agency (JMA), were used for the analysis of seismic tomography. In all, 6,750 earthquakes from January 2002 to March 2014 were used. The results showed that a cluster of slab earthquakes is distributed about 10 km below the plate interface. This suggests that the slab earthquakes occur in the lower part of the oceanic crust near the oceanic Moho within the slab. Moreover, the fault of the slow-slip events corresponds with the plate interface. The results also showed that a low Vp and high Vp/Vs area is distributed within the subducted oceanic crust underlying the fault associated with the slow-slip events. These structures are similar to those in the Tokai district, which suggests that the thermal condition of the southern Ryukyu arc is similar to the case of a hot-slab area where slow-slip events occur.

Findings

Introduction

Slow-slip events (SSEs), events with rupture velocities and durations that are slower and longer, respectively, than those typically associated with earthquakes, have been detected widely in subduction zones (e.g., Obara 2002; Ozawa et al. 2002, 2013; Rogers and Dragert 2003; Kostoglodov et al. 2003; McCaffrey et al. 2008). High-pressure fluids within the subducted oceanic plate affect the occurrence of SSEs (Shelly et al. 2006; Liu and Rice 2007; Audet et al. 2009; Song et al. 2009, Kato et al. 2010), and their repeated occurrence results in stress accumulation in the up-dip part of the fault, which could potentially induce megathrust earthquakes.

Heki and Kataoka (2008) reported the detection of SSEs in the southern Ryukyu Trench, where the Philippine Sea Plate (PHS) is subducting in a northwesterly direction beneath the Eurasian Plate (Figure 1). The SSEs occur repeatedly at intervals of 6 to 7 months. The released moments of the SSEs correspond to Mw 6.6 earthquake faulting, and their fault depths range from 32 to 49 km (Figure 1d).

Figure 1

Seismotectonic setting of the southern Ryukyu and locations of seismic stations. (a) Index map for the studied region. EU Eurasian Plate, PHS Philippine Sea Plate, RA Ryukyu arc, RT Ryukyu Trench. (b) Seismic stations from the Incorporated Research Institute for Seismology (open triangles) and the Japan Meteorological Agency (solid triangles) networks. Circles show the nodes used in the 3-D seismic tomography. IR Iriomote Island. (c) Distribution of relocated hypocenters (circles). Stations IGK, KUROSH, and HATERS are those employed for the receiver function analysis. Black triangles and red triangles denote short-period and broadband seismic stations, respectively. The shaded area shows the fault of repeated SSEs. (d) North-south vertical cross section of relocated hypocenters in the rectangular area denoted as X-Y in (c).

The Ryukyu Trench is a typical subduction zone at which large interplate earthquakes are rare. Both seismic coupling and interplate coupling are assumed to be weak at the Ryukyu Trench; the seismic coupling estimated from large earthquakes is 5% (Peterson and Seno 1984). Measurements by the Global Navigation Satellite System (GNSS) network in the Ryukyu Islands have revealed movement of the stations in the direction of the Ryukyu Trench resulting from rifting of the Okinawa Trough, but these cannot detect the compressional strain accumulation induced by interplate coupling (Nakamura 2004; Nishimura et al. 2004; Watanabe and Tabei 2004). However, a large tsunami struck the southern Ryukyu region in 1771, and the source of this event has been postulated to be an M 8-class thrust-type earthquake along the Ryukyu Trench (Nakamura 2009). Additionally, interplate coupling with the potential to cause M 7.5 to 8.5 earthquakes has been estimated from GNSS observations in the southwestern part of the Ryukyu Trench (Hsu et al. 2012).

In addition, fault depths of the SSEs are shallower than those of slab earthquakes (Figure 1d). The depths of the SSEs have been estimated using the grid search method to realize the smallest post-fit residuals of the GNSS data (Heki and Kataoka 2008). As SSEs usually occur on the plate interface, accurate assessment of the depth of the plate interface is important for both modeling the SSEs and establishing the temperature-pressure relation of the plate interface beneath the southern Ryukyu arc.

This paper investigates the seismic structure of the plate interface and underlying oceanic crust to clarify the relationship between the occurrence of SSEs and fluid in the oceanic crust.

The depth of the plate interface has been estimated to range from 30 to 45 km using the travel times of reflected S-waves that can be observed via the waveforms of local crustal earthquakes (Nakamura 2001). However, this estimation includes uncertainties related to the mis-picking of reflected S-waves (identification between the S-to-S reflected phase and other converted or reflected phases). Thus, receiver function analysis was employed to determine the depth of the plate interface, and the 3-D velocity structure was determined to estimate the velocity structure in and around the SSE faults.

Data and methods

Receiver function

Receiver functions were calculated using teleseismic waveforms recorded by one broadband seismograph station (IGK) of F-net, operated by the National Research Institute for Earth Science, and two short-period stations (KUROSH, HATERS, with a natural period of 1 s) of the Japan Meteorological Agency (JMA) (Figure 1c). In all, 212 earthquakes (M > 6.0) that occurred from January 1, 2000, to March 30, 2014, were used (Figure 2a). The epicentral distances of the earthquakes used ranged from 30° to 90°, and the depths of the earthquakes ranged from 10 to 644 km.

Figure 2

Distribution of epicenters used for receiver function analysis and velocity structure. (a) Epicenter distribution of earthquakes (red circles) used in the receiver function analysis. Triangle shows the position of the seismic station. (b) 1-D velocity structure. Solid and dotted lines show the initial and final models, respectively.

Receiver functions effectively distinguish P-S waves from among all the observed phases that have been converted at seismic velocity discontinuities just beneath a station. Waveforms with durations of 95 s extending from 5 s before to 90 s after the onset of P-waves were used to construct the receiver functions. A cosine taper was applied 5 s from the beginning and end of the waveform, and the instrument response was removed. The water level method (Langston 1977, 1979), in which the value of the water level is defined as 0.1 and 0.01, was applied to the original waveforms in the frequency domain, and a Gaussian filter was applied for 1.6-Hz high-cut filtering of the data.

Assuming that each later phase in the receiver functions represents the wave converted from P to S at depth, the time-domain receiver function was transformed into a depth domain along each ray path using a reference velocity model. To construct the reference velocity model, a local 1-D structure of the area was determined using VELEST software (Kissling et al. 1994) (Figure 2b). The JMA2001 (Ueno et al., 2002) velocity model, which is the standard velocity structure in the Japan region, was employed as the initial model of the velocity structure in VELEST (Figure 2b). For the calculation of the 1-D velocity structure, P and S arrival data of 1,950 earthquakes (M ≥ 2.5), which were picked and listed by the JMA, were used. These earthquakes occurred between January 1, 2002, and March 15, 2014, and their epicenters were distributed within an area defined by 23.7° to 24.8° N and 122.9° to 124.6° E. We used 164 traces to image the receiver functions, selecting only noise-free seismograms after visual inspection. We categorized the area under investigation into blocks with a horizontal scale of 1.5 km and depth of 1.0 km. The receiver function amplitude shown at a certain depth was attributed to the block at the same depth through which the rays passed. When two or more rays passed through the same block, the average amplitude of all the rays was calculated. The results obtained from the receiver functions were compared with those of the 3-D seismic tomography and relocated hypocenters.

Seismic tomography

We used SIMUL2000 software (Thurber and Eberhart-Phillips 1999), which allows for relocation of earthquakes with a 3-D P-wave velocity and Vp/Vs model. To obtain the data, 12 seismic stations with short-period seismometers of the Japan Meteorological Agency, 2 broadband seismometers of F-net, and 5 Taiwan seismic stations of the Incorporated Research Institute for Seismology were used. The station distribution is limited to the Ryukyu Islands and Taiwan, which means that their azimuthal coverage is not advantageous for the tomography. The average interval between stations is approximately 30 km in the Iriomote Island and Ishigaki Island area, but it extends to more than 50 km outside of this area. In all, 6,750 earthquakes (M ≥ 1.2) that occurred between January 1, 2002, and March 15, 2014, were used. The epicenters were distributed within the area of 23.7° to 24.8° N and 122.9° to 124.6° E, and the depths of the earthquakes ranged from 0 to 152 km.

The numbers of used P- and S-wave arrival times, which were picked and listed by the JMA, were 43,898 and 41,134, respectively. The picking accuracy of P-wave arrival times was estimated to be 0.1 s, but the uncertainty of S-wave arrival times was slightly larger, at 0.2 s. The grid intervals in the horizontal direction were set to 10 to 20 km, and the vertical grid intervals were set to 7.5 km for the depth range of 0 to 75 km and 25 to 50 km for greater depths. For the initial 3-D tomography model, a 1-D velocity structure developed using VELEST (Kissling et al. 1994) was used (Figure 2b). The earthquake catalogue used for the 1-D velocity structure was the same as that used for the receiver functions. The P and S travel times were then inverted for Vp and Vp/ Vs perturbations from the reference 1-D velocity model. After four iterations, the average root mean square residual decreased from 0.237 to 0.173 s for the P-phase and from 0.379 to 0.286 s for the S-phase.

Checkerboard resolution tests were then conducted to confirm the reliability of the obtained tomographic images. To make a checkerboard, alternating positive and negative velocity anomalies of 5% were assigned to all the 3-D grid nodes for two grid intervals. Random errors with a standard deviation of 0.1 s for the P-phase and 0.2 s for the S-phase were then added to the synthetic arrival times, calculated for the checkerboard model, to account for picking errors in the real data. We also conducted restoring resolution tests (Zhao et al. 1992) by taking the final 3-D velocity structure as the synthetic model.

Results and discussion

In Figure 3a, two regions of positive amplitude (indicated by red) and one region of negative amplitude (indicated by blue) in the receiver functions can be clearly traced, and lines are drawn as interpretations of P-S converted phases at velocity discontinuities. The negative amplitude region at depths of 30 to 45 km and dipping northward is interpreted as the plate interface. The positive amplitude region at depths of 30 to 35 km is interpreted as the land Moho; the oceanic Moho can be observed between the horizontal distances (Y) of 15 and 65 km. If we change the water level and high-cut frequency, the structure is almost reproduced (Figure 3b,c). However, because the signal-to-noise ratio of the waveform is not particularly high (10 ~ 30 at the frequency of 0.2 to 2.0 Hz), we selected the results for the higher water level of 0.1.

Figure 3

Depth sections of receiver functions beneath line A to B in Figure1c. Gray lines show the boundaries estimated from the receiver function. Dotted line denotes the SSE fault. (a, b) The depth section when the water level is 0.1 (a) and 0.01 (b). (c) The depth section when the water level is 0.1 and the high-cut filter is set to 1.0 Hz.

The plate interface is estimated from the receiver function analysis as dipping northward with an angle ranging from 6° to 14°. The fault locations of the SSEs correspond with the location of the plate interface (Figure 3). The plate interface is in contact with the land Moho at a depth of 30 km, which is the upper part of the SSE faults. The negative amplitudes are obscured from Y = 5 to 20 km. As the plate interface contacts or crosses the overlying land Moho from Y = 5 to 20 km, the S-wave velocity contrast becomes smaller at the plate interface and the negative amplitude of the receiver function weakens.

The area of the low Vp anomaly (-4%), which corresponds to the oceanic crust, is distributed approximately 10 km above the cluster of slab earthquakes at X = -20 km and Y = 10 to 70 km (Figure 4a). However, the area of low Vp anomaly is elongated within and above the cluster of slab earthquakes at the trace of X = 10 km and Y = 0 to 70 km (Figure 4b). The resolution of Vp is good in the area of Y = 0 to 60 km and the depth range of 10 to 50 km (Figure 4e,f). The area of high Vp/Vs (1.78 to 1.83) in the oceanic crust is distributed above the cluster of slab earthquakes (Figure 4c,d). The area of high Vp/Vs is also distributed at depths of 20 to 45 km when Y = 0 to 40 km along the traces of X = -20 km and X = 10 km. The high Vp/Vs region is parallel to the cluster of slab earthquakes. The deep extent of the high Vp/Vs anomaly can be traced to a depth of 60 km. However, the distribution of the high Vp/Vs anomaly, which continues into the deeper region, cannot be determined because of poor resolution (Figure 4g,h).

Figure 4

North-south vertical cross section of seismic velocities. (a, b) Vp at the location of X = -20 km (a) and X = 10 km (b). Gray lines show the boundaries estimated from the receiver function. Black spots denote the relocated hypocenters within a horizontal range of 10 km from the cross section. Dotted line denotes the SSE fault. CMT solutions of the relocated hypocenters are after F-net, NIED. (c, d) Vp/Vs at the locations of X = -20 km (c) and X = 10 km (d). (e, g) The checkerboard resolution test of Vp at the locations of X = -20 km (e) and X = 10 km (g). (i, k) The checkerboard resolution test of Vp/Vs at the locations of X = -20 km (i) and X = 10 km (k). (f, h) The restoring resolution test of Vp at the locations of X = -20 km (f) and X = 10 km (h). (j, l) The restoring resolution test of Vp/Vs at the locations of X = -20 km (j) and X = 10 km (l). Dots in (e), (g), (i), and (k) show the grid points used in this study.

Slab earthquakes dipping at approximately 10° occur at shallow depths within the oceanic crust of the PHS (0 < Y < 40 km) (Figure 4b), whereas slab earthquakes dipping at approximately 25° occur at greater depths within the lower crust near the oceanic Moho (40 < Y < 60 km). The centroid moment tensor (CMT) solutions of the intraslab earthquakes are strike-slip-fault type (down-dip-extension type) (Figure 4a). Slab earthquakes occur on a distinct plane 8 km below the plate interface in the Nankai (Shelly et al. 2006), Tokai (Kato et al. 2010), and Cascadia regions (Hyndman and Wang 1993), where low-frequency earthquakes, tremors, and slow-slip events occur. On average, subducting oceanic crust is characterized by relatively low velocities and high Vp/Vs ratios, which is interpreted as being indicative of fluids released by progressive metamorphic dehydration reactions in the oceanic crust (Hacker et al. 2003). As the occurrence of intraslab earthquakes is controlled by thermo-petrologic conditions (Hacker et al. 2003), the similarity in occurrence of the slab earthquakes studied here suggests that the thermal conditions in the southern Ryukyu arc are similar to those in the Nankai, Tokai, and Cascadia regions. Yamasaki and Seno (2003) divided the morphology of the double seismic zone into hot- and cold-slab types, for which dehydration occurs at depths of approximately 45 km and more than 75 km, respectively. They estimated that dehydration reactions occur mostly before a depth of approximately 60 km in the Taiwan (southern Ryukyu) region and before a depth of approximately 45 km in the hot-slab areas of the Nankai region and Cape Mendocino region south of Cascadia. Similar thermal conditions in the southern Ryukyu region as in the Nankai and Cascadia regions would generate repeated SSEs in the southern Ryukyu Trench.

McCormack et al. (2013) showed that a serpentinized layer, whose thickness is 6 km, overlies the upper slab interface of subducted PHS beneath the Ryukyu arc. However, the serpentinized layer would likely be a region of moderate Vp/Vs (1.75 ~ 1.83) at a depth of 40 km (Christensen 2004). Moreover, the thickness of the serpentinized layer estimated by McCormack et al. (2013) is only 6 km, which is smaller than the node intervals (horizontal 10 km; vertical 7.5 km) used in this study; thus, the resolution is too low to distinguish a 6-km-thick layer within the surrounding mantle.

Chou et al. (2009) investigated the velocity structure of the mantle wedge above the subducted PHS in the Taiwan-Ryukyu junction zone (121.8° E) and showed that approximately 60% serpentinized mantle wedge lies on the PHS at a depth of 50 km. Thus, a serpentinized area is distributed in the mantle wedge beneath the southern Ryukyu arc. However, serpentinized mantle wedge could only be detected near Taiwan (121.8° E); it could not be detected at 122.6° E, 80 km east of Taiwan, because the resolution of the tomography was poor. Because station distribution is limited to islands, the detection of serpentinized mantle beneath ocean areas would be difficult unless station distribution were improved in the back-arc region.

The depths of the crustal earthquakes are restricted to 0 to 25 km (Figure 4a,b). The lower depth limit of the crustal earthquakes is 15 km in the northern region (Y = 60 km) and 25 km in the southern region (Y = 10 km). A previous study estimated the depth of the land Moho, based on an analysis of arrival times of the Pn-phase, to be 25 to 35 km beneath Iriomote Island (Nakamura and Umedu 2009), which is similar to the depth of the land Moho estimated from the receiver function analysis.

Areas of low Vp/Vs are distributed at a depth of 20 km in the X = -20 km and Y = 10 km and X = 10 km and Y = 20 km regions, which are above the area of high Vp/Vs of the subducted oceanic crust. The Vp anomaly is small at X = -20 km and Y = 10 km, whereas high Vp is distributed at X = 10 km and Y = 20 km. These areas could be reproduced by the restoring resolution test, although the anomaly tends to be weak (Figure 4c,e). A high Vp body with an area of low Vp/Vs in the crust, which would be an impermeable layer, was observed in the Tokai region (Kato et al. 2010). Additionally, in Hikurangi in New Zealand, Reyners and Eberhart-Phillips (2009) observed a cap-rock structure of a low Vp/Vs schist layer overlying high Vp/Vs oceanic crust. They interpreted this structure as an impermeable barrier or cap rock that restricts fluid movement and generates fluid overpressure (Audet et al. 2009; Reyners and Eberhart-Phillips 2009; Kato et al. 2010). The fluid beneath the SSE region of the southern Ryukyu arc would be generated under similar conditions. Adjacent and to the north of the low Vp/Vs body, a high Vp/Vs body is distributed. This suggests that the Vp/Vs structure of the lower crust changes at around Y = -20 km. As the resolution of the southern part of the area is poor, relatively large ambiguity exists for the southward extension of the low Vp/Vs area.

These structures are similar to that in the Tokai district. An area of high Vp/Vs (approximately 1.95) in the oceanic crust and an overlying impermeable layer have been detected in the area of long-term SSEs in the Tokai district. Overpressured fluid is assumed in the oceanic crust of the area of long-term SSEs. Similar to the case of Tokai, an impermeable layer overlying the area of high Vp/Vs oceanic crust occurs in southern Ryukyu. However, in this case, the degree of the Vp/Vs anomaly is not large. This suggests that the generated fluid volume is small because dehydration of the oceanic crust is not as active as it is in the Tokai district. Unfortunately, regarding possible alternatives for the discrepancy in the degree of the Vp/Vs anomaly, the fine-scale structure of the slab and mantle wedge cannot be resolved because the seismic stations are too sparse and the node size is too large. A dense seismic survey in the SSE area of the southern Ryukyu would show the fluid distribution in the SSE fault area, which would reveal the relation between the fluid and the generation of SSEs in the Ryukyu subduction zone.

Conclusions

Based on fine-scale variations of seismic velocities and converted teleseismic waves, the tectonic similarities between SSEs in the southern Ryukyu arc region and the hot-slab area where SSEs occur are revealed. Slab earthquakes occur about 10 km below the plate interface where the SSEs occur, which suggests that the slab earthquakes occur in the lower part of the oceanic crust, near the oceanic Moho in the slab. The results also show that an area of low Vp and high Vp/Vs is distributed in the subducted oceanic crust underlying the fault of the SSEs. The impermeable structure of the land crust caps the area of low Vp and high Vp/Vs of the subducted oceanic crust, which would generate high fluid pressure in the oceanic crust. The similarities of this structure with cases in the Tokai district imply that the thermal condition of the southern Ryukyu arc is similar to the hot-slab area where SSEs occur.