Interface geometry of potential mega-thrust earthquakes beneath the westernmost Ryukyu subduction system

The interface geometry at the subduction boundary plays one of the most important roles on evaluation of both seismic and tsunami impacts as the mega-thrust earthquakes occur between two plates. Although the general feature of the subducted slab at the mantle depth is clearly delineated from the seismicity in the Wadati-Benioff zone, the exact interface geometry at the crustal depth is hardly obtained from abundant earthquakes scattering in and around the interface between the two plates. Examination of seismic data recorded at the dense seismic array in the Tatun volcano group of northern Taiwan shows two unambiguous P-waves generated by shallow earthquakes offshore the Hualien area in eastern Taiwan. The detailed analyses of travel-times of both P-waves show that the 1st P-wave was the direct wave propagating within the upper crust while the 2nd P-waves were reflected from the subducted slab dipping northward ~23 degrees. This observation reveals the general interface geometry between the subducted Philippine Sea and the overlaid Eurasian plates at the crustal depth, where is often considered as the locked zone and accumulated a lot of seismic energy for release. Thus, the interface geometry obtained here provides important parameters for estimating the possible rupture plane along the mega-thrust fault between two plates in the westernmost Ryukyu subduction system in the future. Article history: Received 2 August 2017 Revised 21 January 2018 Accepted 22 January 2018


IntRoductIon
It is well known that the largest earthquakes in the Earth (Ms > 9.0), such as the 1960 Chile earthquake, 1964 Alaska earthquake, the 2004 Indonesia earthquake, and the 2011 Japan earthquake, were always taken place along the mega-thrust fault in the subduction system (Plafker 1972;Kanamori and Cipar 1974;Cifuentes 1989;Briggs et al. 2006;Chu et al. 2011;Kurahashi and Irikura 2011).Since the mega-thrust fault is the boundary between the subducted slab and the overlaid plate, strong convergent energy will be released by mega-earthquakes along a low-angle interface between two plates.The mega-thrust might be rupture starting from the trench down to the seismogenic zone at depth around 40 -50 km in general (Hyndman et al. 1997).Thus, the detailed geometry of the interface between two plates, such as the dipping angle and depths, plays on an important role to estimate rupture behaviors and their impact of both ground motions and tsunamis.
However, the interface between two plates at the crustal depth is not easy to be identified from seismicity directly because the seismicity in the shallow part is extremely scattering than that in the deep subduction zone.It is generally clear that seismicity in the Wadati-Benioff zone often provides clear information to delineate the interface of the subducted slab in that all of earthquakes basically occur within the subducted slab.The slab interface is approximately plotted on the upper limitation of seismicity according to thermal structures (Abers et al. 2006).On the other hand, crustal earthquakes occurred not only along the plate boundary but also scattered in both of the overlaid plate and the subducted plate even although major convergence rate was often taken place in the Terr.Atmos. Ocean. Sci., Vol. 29, No. 4, 405-415, August 2018 plate boundary on the surface.Thus, it is hard to find the interface geometry between two plates from seismicity since a lot of scattering earthquakes were detected at the shallow depths.As a result, some other methods might be developed to identify the interface geometry between the subducted slab and the overlaid plate in the subduction zone.
In the Taiwan area, it won't be surprising to see a lot of earthquakes taken place in that there are two subduction systems (Tsai et al. 1981;Lin 2002).Among the abundant seismicity, most of the large earthquakes might be resulted from ruptures along the subduction interface between two plates.For instance, one of the large earthquakes on 31 March 2002 (often called the 331 Earthquake, Table 1) was rupturing along the subduction interface in the northeastern Taiwan area and caused some damages in Taiwan (Chen 2003;Chen et al. 2004).It is obviously that such an interface will potentially rupture again and might generate larger earthquakes in the future.Although nobody has any idea when the mega-thrust earthquakes will occur exactly, we have to carefully prevent the possible impact by knowing all possible source parameters.
In order to improve the understanding of interface geometry between the subducted slab and the overlaid plate in eastern Taiwan, we carefully examine the seismic data recorded by a dense seismic array in the Tatun volcano group for identifying particular later phases.Those later phases might be associated with the major boundary in the subduction zone at the crustal depths.For comparing with the observations, a well-developed traditional ray-tracing algorithm (MacRay, developed by J. Luetgert) was employed to calculate the travelling times recorded at the dense array.Finally, the interface geometry between the PSP and EUP at the crustal depth is obtained beneath the westernmost Ryukyu subduction system.This result provides the important parameters for estimating possible impact of both seismic shaking and tsunami by a megathrust earthquake in the future.

tectonIc bAckgRound
The island of Taiwan is located at a portion of convergent zone between the Eurasian plate (EUP) and Philippine Sea plates (PSP) (Fig. 1).The PSP subducts beneath the EUP in the northeastern Taiwan area while the EUP underthrusts the PSP in the southern Taiwan area.The suture zone between two plates is identified in the island of Taiwan from Hualien (~24.0°N) to Taitung (~22.6°N)along the Longitudinal Valley, which is often considered as the plate boundary between two plates.Such a suture zone in the island of Taiwan is clearly connected to the westermost Ryukyu trench offshore Hualien.The subduction slab beneath the NE  Taiwan area can be clearly delineated by the deep earthquakes in the Wadati-Benioff zone (Fig. 1b).A clear interface in the upper mantle is unambiguously marked by a thick-red line at the depths greater than 40 km, but it is hard to see the interface geometry from the abundant earthquakes scattering in and around the plate boundary.The dipping angles of the subduction slab at the shallow depths have a lot of uncertainties, approximately ranging from 10 -40 degrees.

seIsMIc dAtA
For monitoring volcanic activities and improving the understanding of the volcanic characteristic at the Tatun volcano group (TVG) in the northern Taiwan area (Lin et al. 2005a, b;Lin 2016Lin , 2017a, b;, b;Lin and Pu 2016), we have deployed a dense seismic array since 2003 (Fig. 2).This seismic array has been gradually updated in the past decade.The seismic array is basically composed of 4 groups including YM, YC, YL, and KL for their own particular purposes and deployed at different times.At first, we started to deploy the YM group around the central part of Yang-Min-Shan National Park for detecting a lot of volcanic earthquakes and tremors in the TVG.Then the YC group is roughly a circle outside the YM group for increasing the ability of detecting seismicity.The YL group is basically alienated along the NW-SE direction for recording abundant earthquakes in eastern Taiwan, such as the Hualien and Ilan areas (Fig. 3).In other words, the source and receiver will be roughly located at a two dimensional profile for examining the deep structures across the northern Taiwan area.The KL group is deployed along the northern coast near the Keelung city for detecting possible seismicity associated with the submarine volcano offshore the Keelung harbor.In total, the seismic stations in the TVG are composed of 40 broadband seismic stations within a small aperture now.Each of seismic stations is equipped with a 3-component sensor (Garulp 6TD) and recording seismic data with a sampling rate of 100 Hz.
It is very interesting to see that two P-waves (P1 and P2) generated by a felt earthquake (M 3.8) at the eastern Taiwan area were well recorded at many seismic stations in the TVG.This earthquake was located by the Central Weather Bureau offshore the Hualien area (121.97°E,24.52°N) with the focal depth of 13.1 km (Table 1).A plot of vertical seismograms with the epicentral distances between 74 and 90 km shows that two P-waves separated by ~1 sec were roughly identified at the entire seismic array, even although some seismograms were overlap on each other (Fig. 4).The detailed arrivals of the 2 nd P-waves are marked in Fig. 5.This feature can be checked in details at either different epicentral distances or seismic groups.First, examination of the epicentral distances ranging from 78 -81 km along a roughly two-dimensional profile shows the 2 nd P-waves were almost arriving ~1 sec later (Fig. 6).Similar features were recorded at other stations at groups of YC and KL at the epicentral distances from 74 and 90 km (Fig. 7) as well as YL group (Fig. 8).All of those observations consistently show that two P-waves separated by ~1 sec were identified at most of seismic stations within the seismic network.
A similar feature of two P-waves (P1 and P2) was also recorded from another earthquake (Event 2 on Table 1).The arrival-time difference between P1 and P2 were identical of ~1 sec at the epicentral distances between 81 and 94 km (Fig. 9).Some arrivals of the 2 nd P-waves are clearly shown in Fig. 10.Although both earthquakes were separated by ~15 km in the horizontal distance (Fig. 3), the consistent features of two P-waves recorded at the dense seismic array in northern Taiwan strongly indicate that some major structures exist beneath the eastern Taiwan area.The 2 nd P-waves were hardly generated around the receiver-side structures since the delays between the direct P-waves and the 2 nd Pwaves were almost identical.Any lateral structures nearby the receiver-sites would generate the significant delay time variations recorded at the whole seismic array.Also, the observation of the 2 nd P-waves with seismic amplitudes greater than the direct waves rule out the possibility of any later phases generated by the layer structures beneath the seismic network in the TVG area.In other words, some multiple paths were produced by the subsurface structures between the earthquakes and the dense seismic array.

RAY-tRAce ModeLIng
In order to check the possible ray-paths for explaining both P-waves (P1 and P2) observed above, we have employed a well-developed software, called MacRay (Luetgert 1992), to calculate the traveling times and ray-paths by assuming a simplified 2-D velocity model.The P-wave velocity is given by 5.0 km s -1 on the surface, and it gradually increases to 6.5 km s -1 at the bottom of the crust and jumps to 8.0 km s -1 at the uppermost mantle.A variety of forward modeling calculations have been conducted to get the best model for fitting the observations of both P-wave arrivals recorded at the seismic array in the TVG.At first, both P-waves (P1 and P2) recorded at the epicentral distances between 75 and 90 km were too short to be the Pwave refracted from the Moho-discontinuity.Also it was almost no way to find any horizontal reflection boundary for fitting the observations because the delay times between two P-waves recorded at the distances from 75 -90 km were almost identical (~1.0 sec).It is obvious that the travel-time difference between the direct and reflect waves increases with the epicenter distances if we consider a horizontal reflection boundary (Fig. 11a).Alternatively, we have considered a dipping boundary to fit the later P-waves with almost identical delays (~1 sec) recorded at the entire seismic array.The final result suggests that the 1 st P-wave was the direct waves while the 2 nd P-wave (PsP) was reflected from the boundary dipping ~23 degrees northward with some        1).uncertainties of a few of degrees (Fig. 11b).The 2 nd P-waves with larger amplitude might be generated by the variation of the source radiation.In order to confirm this possibility, we have determined the focal mechanism from the 1 st -motions (Fig. 3) and then projected along the depth-profile later.The results show the 2 nd P-waves reflected from the subduction slab had a larger amplitude generated at the source.On the other hand, the first direct P-waves with smaller amplitudes were come around the nodal plane.
As we mentioned early, similar feature of two P-waves separated by ~1 sec was also found from another earthquake (M = 3.0) offshore the Hualien area.The 2 nd earthquake was located nearby the 1 st earthquake with almost the same depth (Table 1 and Fig. 3).The epicentral distances from the 1 st earthquake was ranging from 75 -90 km, while that from the 2 nd one was ranging from 79 -94 km.No matter what epicentral distances were, two P-waves separated by ~1 sec were consistently observed at most seismic stations generated from two earthquakes (Figs. 5 and 10).

dIscussIon
A dipping boundary obtained from the forward modeling above reveals the interface geometry of the subducted slab beneath the westernmost Ryukyu subduction system (Fig. 11).The later P-waves were reflected from the top of the subducted PSP at the depths of 25 -30 km beneath the location of ~24.5°N.The dipping angle of ~23 degrees is found at the depths of 25 -30 km beneath ~24.5°N.The southward projection of such a dipping angles beneath 24.5°N on the surface is around the Hualien area between 23.5 and 24.0°N, which is roughly consistent with the western end of the Ryukyu trench offshore the Hualien area of eastern Taiwan.It is also very interesting to see that the location of the 331 earthquake in 2002 is approximately overlaid on the interface with the dipping angle of ~23 degrees at the crustal depth (Fig. 12).Such a result is generally consistent with one of the previous studies (Wu et al. 2009) that indicated the dipping angle of the subducted slab was around 20 degrees based on the seismicity.Thus, those consistent features indicate that the interface geometry obtained in this study is an acceptable result, even although some of uncertainties might be existed.
Dramatic variation of the dipping angles changing from ~23 degrees at the crustal depth to 60 -70 degrees in the uppermost mantle provides valuable information to estimate the largest possible earthquakes beneath the northeastern Taiwan.It is very clear that the subduction angle at depths below ~40 km is about 60 -70 degrees based on the seismicity in the subduction zone (Fig. 12).Such a high dipping angle is hardly taken place on the mega-thrust earthquakes according to both earthquake physics and observations.Instead, the dipping angle of ~23 degrees is more often observed and acceptable to generate the mega-thrust earthquakes in the subduction zone.Thus, the interface geometry dramatically changed from the low to high angles might be considered as the deepest part of the seismogenic zone for generating mega-thrust earthquakes in the future.Some preliminary simulation of the potential impact by such a mega-thrust earthquake can be modeled for evaluating seismic and tsunami threats in the island of Taiwan.The more reliable results might be improved by examining seismic data generated in and around the Hualien area in the future.

concLusIon
Examination of seismic data recorded at the dense seismic array in the Tatun volcano group shows two unambiguous P-waves were respectively generated by two shallow earthquakes offshore the Hualien area.Careful modeling of travel-time arrivals and ray-paths along a two dimensional profile indicates that the first P-wave was the directed wave propagated within the upper crust while the second P-wave was reflected from the subducted PSP slab beneath the westernmost Ryukyu subduction system.The general pattern of the subduction interface can be delineated here.The subduction of the PSP starts around 24.0°N with a dipping angle of ~23 degrees toward the north.The subduction angle dramatically increase to 60 -70 degrees at the depth of ~100 km beneath the Ilan area (25.0°N).The interface geometry dramatically changed from the low to high angles might be considered as the deepest part of the seismogenic zone for generating mega-thrust earthquakes.In other words, some large earthquakes, such as 331 earthquakes in 2002, might be taken place along the dipping angle of ~23 degrees in the future.Potential impact by such a mega-thrust earthquake beneath the westernmost Ryukyu subduction system can be modeled for evaluating seismic and tsunami threats in the island of Taiwan.
Fig. 1.(a) Simplified tectonics in the Taiwan area and (b) Seismicity profile beneath the NE Taiwan area, marked in green at (a).The red-lines marks the possible subduction boundary between two plates.

Fig. 3 .
Fig. 3. Locations of YL seismic stations and two felt earthquakes (Events 1 and 2) occurred in December 2015.The focal mechanisms determined by the 1 st -motions are also shown here.

Fig. 2 .
Fig. 2. Locations of seismic array (triangles) in the Tatun volcano group.The insert map shows the general tectonics in the Taiwan area.The Tatun volcano group is marked a black square and the plate boundaries of the Manila and Ryukyu Trenches is shown by the thick lines with triangles.

Fig. 4 .Fig. 5 .
Fig. 4. Vertical Seismograms recorded at the TVG and plotted at the epicentral distances for roughly showing two P-waves.The detailed plots are separated into 3 groups in Fig. 5.

Fig. 7 .
Fig. 7. Vertical seismograms recorded at the YC and KL seismic stations and plotted at the epicentral distances.

Fig. 8 .
Fig. 8. Vertical seismograms recorded at the YL group and plotted at the epicentral distances.

Fig. 10 .
Fig. 10.Vertical seismogram recorded at YL stations.The arrivals of the 2 nd P-waves are marked by grey boxes.

Fig. 12 .
Fig. 12. Seismicity profile beneath northeastern Taiwan marked in green on the insert map.The plate boundary between the subducted slab and the overlaid plate is shown in a thick-red line.The locations of both Event 1 used here and the 331 earthquake in 2002 are marked in the white circles.The focal mechanism of Event 1 is also plotted on the depth profile for showing seismic amplitude variations between the direct (Pg) and reflect waves (PsP).

Table 1
. Earthquake parameters (Provided by Central Weather Bureau).