Imaging the subducting slabs and mantle upwelling under the Japan Islands

Summary

A high-resolution P-wave tomography of the crust and mantle down to 700km depth beneath the Japan Islands is determined using a large number of high-quality arrival-time data from local earthquakes and teleseismic events simultaneously. The tomography shows that the Philippine Sea slab is subducting aseismically down to 430km depth under southwest Japan, though the seismicity within the slab ends at 180km depth. A low-velocity (low-V) zone in the mantle wedge under Tohoku and Kyushu is found to extend westward from the volcanic front to the backarc under the Japan Sea and East China Sea. Significant low-V anomalies are revealed in the deep portion of the mantle wedge (400–500km depth) above the Pacific slab under southwest Japan, which may reflect hot mantle upwelling associated with fluids from the deep dehydration of the Pacific slab. Low-V anomalies appear at 420–700km depths beneath the Pacific slab under eastern Japan, which may reflect hot mantle upwelling associated with the deep subduction of the Pacific slab and its collapsing down to the lower mantle.

1 Introduction

Four lithospheric plates exist in and around the Japan Islands (Fig. 1). The Pacific Plate is moving toward the northwest and is subducting beneath Hokkaido and northern Honshu from the Kuril and Japan trenches, whereas the Philippine Sea (PHS) Plate is moving northwestwards and is descending beneath southwest (SW) Japan from the Nankai trough. Hokkaido and northern Honshu belong to the Okhotsk (or North American) Plate, whereas SW Japan belongs to the Eurasian (or Amur) Plate (Bird 2003; DeMets et al. 2010). The strong interactions among the four plates cause intense seismic and volcanic activities in and around the Japan Islands. A striking recent example of the intense seismicity in the region is the great 2011 Tohoku-oki earthquake (M_w 9.0) that took place in the interplate megathrust zone under the northeast (NE) Japan forearc because of the subduction of the Pacific Plate beneath the Okhotsk Plate (Fig. 1). There are 110 active arc and backarc volcanoes on the Japan islands, in addition to several active intraplate volcanoes in the NE Asian continental margin, for example, Changbai, Jeju and Ullung, etc. (Siebert & Simkin 2002) (Fig. 1). It is of great importance to investigate the detailed 3-D crustal and mantle structure of this typical subduction zone so as to better understand the seismotectonics, volcanism and subduction dynamics.

So far many researchers have used seismic tomography to investigate the 3-D seismic velocity structure beneath Japan (for details, see recent reviews by Hasegawa et al. 2009; Zhao et al. 2011). However, most of the previous tomographic studies used only the arrival-time data from local earthquakes in the crust and the subducting Pacific and PHS slabs under Japan, which could reveal the structure down to about 200km depth including the crust, upper-mantle wedge and the upper portion of the Pacific slab, but could not determine the deeper structure for the lower portion of the Pacific slab and the mantle below the slab. This problem can be resolved by adding data from teleseismic events to conduct a joint inversion of the local and teleseismic data, but so far only a few studies have taken this approach to study the deep structure under the Japan islands (e.g. Zhao et al. 1994; Abdelwahed & Zhao 2007). One drawback of the few previous studies is that the teleseismic data used were mainly recorded by the old and sparse J-array seismic network (Yomogida 1996) and so the teleseismic data used were not sufficient in both quantity and quality.

In this work, we have made great efforts to collect a large number of high-quality teleseismic data from the seismograms recorded by both the J-array and the High-Sensitivity Seismic Network (Hi-net) that covers the Japan islands densely and uniformly (Fig. 2a). The local and teleseismic data are used jointly in the tomographic inversion, which results in a high-resolution 3-D P-wave velocity model down to 700km depth under Japan. The present results shed new light on the deep structure and dynamics of the Japan subduction zone.

2 Data and Method

In this study we used two sets of data to conduct tomographic inversions. One is arrival times from local earthquakes which occurred in and around the Japan islands (Fig. 2b). The other is relative traveltime residuals from teleseismic events (Fig. 3).

The local earthquakes were recorded by the updated Kiban Seismic Network that includes over 150 seismic stations operated by the Japan Meteorological Agency (JMA), ∼350 stations operated by the Japanese national universities and over 700 Hi-net stations operated by the National Research Institute for Earth Sciences and Disaster Prevention in Japan (Okada et al. 2004). The short-period Hi-net stations are installed throughout Japan (Fig. 2a). Each Hi-net station is equipped with three-component velocity-type seismometers of natural frequency of 1 Hz. Most of the Hi-net seismometers are installed at the bottom of a borehole with a depth of 100–200 m. The data are digitized at each station with a sampling rate of 100 Hz and then transmitted to the data center after being attached with the absolute time information from the global positioning system clock. Among the great number of earthquakes recorded by the Kiban seismic network, we selected a best set of events for our tomographic study. The study area is divided into many cubic blocks with a size of 30 × 30 × 20km³. The local events are selected by a specific scheme of selection according to the maximum number of recording stations and the minimum formal uncertainty of the hypocentral location. The shallow offshore events were removed from the data set to exclude the poorly located events. As a result, 1180 local events were selected (Fig. 2b), which generated ∼253 500 P-wave and ∼101 200 S-wave arrival times. The uncertainty of hypocenter locations is <4km for all the events and <2km for the events beneath the land areas where the seismic network is located. The accuracy of the P-wave arrival times is 0.05–0.10 s and that of the S-wave data is ∼0.10 s, according to the seismic network staffs who collected the data.

We collected ∼45 400 first P-wave arrivals from 360 teleseismic events (Fig. 3), among them ∼33 900 data from 333 events were recorded by the J-array seismic network during 1988 November to 2005 June. The remaining ∼11 500 data from 27 teleseismic events were collected from the high-quality seismograms recorded by the Hi-net during 2001 January to 2010 October. The picking accuracy of the P-wave arrivals is estimated to be ∼0.10 s. The size of the teleseismic events selected ranges from M 6.0 to M 8.0. The event epicentral distances are in the range of ∼25°–95°. The azimuthal coverage of these events is fairly good in all directions except for the Pacific Ocean. Most of the events occurred in the subduction zones in northern Pacific, Fiji-Tonga, eastern and southern Asia, North America and the Eurasian continent (Fig. 3). The hypocentral parameters of the teleseismic events determined precisely by E. R. Engdahl are used in this study (for the event location procedure, see Engdahl 2006).

We used the tomographic method of Zhao et al. (1994, 2009) to invert the relative traveltime residuals of the teleseismic events and local-earthquake arrival times simultaneously. This method includes a forward scheme that can deal with a general velocity model in which complex velocity discontinuities exist and seismic velocity changes in three dimensions (see Zhao et al. 1994 for details). Following the previous tomographic studies in Japan (e.g. Zhao et al. 1994; Abdelwahed & Zhao 2007; Nakamura et al. 2008; Huang et al. 2011; Tong et al. 2012), the Conrad and Moho depth variations and the subducting Pacific slab are taken into account in the modelling space (Figs S2 and S3). Previous studies show that the Pacific slab becomes stagnant in the mantle transition zone under East Asia (e.g. Zhao 2004; Huang & Zhao 2006; Zhao et al. 2009, 2011; Gu et al. 2012), which was considered in the starting velocity model for the 3-D tomographic inversion in this study. Two sets of 3-D grid nodes are arranged in the modelling space: one set is in the crust and upper-mantle wedge, the other set is in the subducting Pacific slab and the subslab mantle (Figs S1 and S2). Hypocentral parameters of local earthquakes and velocity perturbations at the grid nodes from a starting velocity model are taken as unknown parameters. The velocity perturbation at any point in the model is calculated by linearly interpolating the velocity perturbations at the eight grid nodes surrounding that point. An efficient 3-D ray tracing technique (Zhao et al. 1994) is used to calculate traveltimes and ray paths accurately and rapidly. The large and sparse system of observation equations was solved by using the LSQR algorithm (Paige & Saunders 1982) with damping and smoothing regularizations (Zhao 2001, 2004; Zhao et al. 2009).

3 Analysis and Results

Applying the tomographic method to our data set, we conducted many tomographic inversions by changing the grid interval, the starting velocity model and the damping and smoothing parameters. We used only the S-wave data with short epicentral distances for locating the local earthquakes together with the P-wave data. The hypocentral parameters of the local events, in particular, those of deep earthquakes outside the seismic network (Fig. 2b), can better be constrained when the S-wave data are used. Figs 4–6 show the map views of our optimal 3-D P-wave velocity model obtained, which is selected from many models we have determined. In this tomographic model, the grid interval is 0.33° in the lateral direction and 5–30km in depth (see Fig. S1). Note that 1°= 111.2km in the north–south direction and 78.6–96.3km in the east–west direction for this study region. Fig. S4 shows the ray path coverage in different depth ranges under the study region. Figs 7–10 show 36 vertical cross-sections of the tomographic image from the surface down to 700km depth under different portions of the Japan islands.

One drawback of the LSQR algorithm (Paige & Saunders 1982) is that no formal resolution matrix is calculated and so all resolution analyses must be made based on the synthetic data tests. In this work we have conducted extensive resolution tests to examine the reliability and resolution scale of the tomographic images obtained. One such test is the checkerboard resolution test that is conducted with the following procedure. We first assign alternative positive and negative velocity anomalies to the 3-D grid nodes in the model and then construct a synthetic data set by calculating theoretical traveltimes of the rays passing through the checkerboard model. The numbers of seismic stations, events and ray paths in the synthetic data set are the same as those in the real data set. To simulate the picking errors contained in the observed data, we added random errors with a standard deviation of 0.1 s to the synthetic traveltime data. The synthetic data are then inverted by applying the same tomographic method to produce a recovered image of the input checkerboard model. The resolution is considered to be good in areas where the checkerboard pattern is recovered well. Results of three checkerboard resolution tests with grid intervals of 0.33°, 0.5° and 1.0° are shown in Figs S5–S7. The checkerboard pattern is recovered down to 200, 500 and 700km depths for the tests with grid intervals of 0.33°, 0.5° and 1.0°, respectively. Hence we think that the spatial resolution of the tomographic model is 30–35km at depths of 0–200km, ∼50km at depths of 200–500km and ∼100km at depths of 500–700km. In addition to these checkerboard resolution tests, other synthetic tests were also conducted, see the following section for details.

The overall pattern of the tomographic results obtained by this study (Figs 4–10) is generally consistent with that of the previous local tomography studies for the shallow parts (0–200km depth) and with that of Zhao et al. (1994) and Abdelwahed & Zhao (2007) for the deeper areas. But the images are much improved because of the much better data set used in this work, as well as some improvements in technical details. For example, in Abdelwahed & Zhao (2007), the images of the subducting Pacific slab look zigzag, which is an artefact due to the large intervals (35–40km) between the grid nodes adopted. This problem was resolved in this study by imposing a constraint on the geometry of the Pacific-slab upper boundary as that shown in Fig. S3 and so the Pacific slab looks smooth in the present tomographic images (Figs 7–10).

We also conducted a joint inversion of the local and teleseismic data using a homogeneous 1-D starting velocity model that does not include the subducting Pacific slab. The results show that the overall features of the tomographic image (Figs S8–S14) are similar to those shown in Figs 4–10 except for the subducting Pacific slab. The Pacific slab is generally imaged as high-velocity (high-V) zones, but the high-V zones are intermittent and scattered in many places and they are not well consistent with the Wadati–Benioff deep seismic zone (Figs S8–S14). This problem was caused by the imperfect ray path coverage in and around the Pacific slab because the seismic stations are located only on the narrow Japan islands which are 200–300km wide in the trench-normal directions (Figs 1 and 2a). According to many previous studies (see the reviews by Hasegawa et al. 2009 and Zhao et al. 2011), the Pacific slab has a thickness of 80–90km and a P-wave velocity 4–6 per cent higher than that of the surrounding mantle, hence the Pacific slab is the most significant velocity anomaly in the mantle under the Japan islands. However, as shown in Figs S8–S14, the Pacific slab itself cannot be fully recovered by the tomographic inversion with a homogeneous 1-D starting velocity model, because of the imperfect distributions of the seismic stations, local and teleseismic events and the ray paths (Figs 1–3). Introducing the Pacific slab into the starting velocity model (Fig. S2) makes it possible to calculate the theoretical traveltimes and ray paths more accurately and the tomographic inversion can be conducted from a more realistic starting velocity model, which leads to a better final 3-D velocity model (Figs 4–10). That is why so far most of the tomographic studies of the Japan subduction zone have adopted a starting velocity model including the subducting Pacific slab (Fig. S2) and so they have obtained reasonable tomographic results for understanding the subduction dynamics under the region (e.g. Zhao et al. 1994; Hasegawa et al. 2005, 2009; Abdelwahed & Zhao 2007; Nakamura et al. 2008; Huang et al. 2011; Tong et al. 2012).

4 Discussion

4.1 Aseismic subduction of the PHS slab

So far many researchers have used local-earthquake tomography to image the subducting PHS slab down to about 200km depth under central and SW Japan (e.g. Hirahara 1981; Ishida & Hasemi 1988; Sadeghi et al. 2000; Salah & Zhao 2003; Wang et al. 2006; Matsubara et al. 2008; Nakajima et al. 2009). The PHS slab images in the shallow mantle revealed by this study (Figs 9 and 10) are generally consistent with those by the previous local-tomography studies.

A new feature revealed by this study is that a high-V zone appears below the seismic portion of the PHS slab under SW Japan (Fig. 10). The seismicity within the PHS slab ends at 60–80km depth under western Honshu and at ∼180-km depth under Kyushu, whereas the high-V zone is visible down to ∼370km depth under the Japan Sea off western Honshu and down to ∼430km depth under western Kyushu (Fig. 10), which may reflect the aseismic subduction of the PHS slab. The aseismic PHS slab is also imaged clearly by the tomographic inversion with the 1-D homogeneous starting model (Fig. S14). To confirm this feature, we conducted detailed synthetic tests (Figs S15–S17). The procedure of the synthetic tests is the same as that of the checkerboard resolution tests; the only difference between them is in the input model (e.g. Fig. S15). The test results show that the high-V zone under SW Japan can be well recovered and so the aseismic PHS slab is a reliable feature (Figs S16 and S17c,d). More synthetic tests and tomographic inversions with different initial models were conducted by Yanada (2011).

Fig. 11 shows the geometry of the aseismic PHS slab under SW Japan estimated from the present tomographic images. The aseismic PHS slab looks discontinuous and intermittent (Fig. 11). There are two possibilities for that. One is that the whole aseismic PHS slab was not fully imaged because the teleseismic ray path coverage is not very dense and uniform under the marginal seas due to the lack of seismic stations in the seas. The other possibility is that the PHS slab itself is segmented and broken and parts of the slab have detached and collapsed down to the deeper mantle. The teleseismic data recorded by the seismic stations on the Korea Peninsula would help to clarify this issue.

The PHS Plate is one of the marginal sea complexes in the western Pacific and it started to subduct northwestwards about 40 Ma ago when the Pacific Plate changed its direction of motion from north-northwest to west-northwest (Karig 1971; Seno & Maruyama 1984). Along the Nankai trough off SW Japan, the PHS Plate is composed of several blocks with ages increasing from the east to west, which are the Izu–Bonin arc and backarc (0–2 Ma), Shikoku basin (15–30 Ma), Kyushu–Palau ridge and Amami Plateau (40–49 Ma) (Hilde & Lee 1984; Hall et al. 1995). The subduction rate of the PHS Plate also increases from 28 mm yr^–1 at the Sagami trough to 50 mm yr^–1 off southern Kyushu (Seno et al. 1993; Bird 2003). It has been inferred that the subduction of the PHS Plate occurred intermittently through the late Cenozoic, depending on the rates and modes of destruction and formation of plates in the PHS (Uyeda & Miyashiro 1974; Hirahara 1981).

Considering the aseismic PHS slab from the present tomographic result, the total length of the subducting PHS slab under SW Japan is estimated to be 800–900km from the Nankai trough to the bottom of the aseismic slab (Figs 10 and 11). Suppose this part of the PHS slab started its subduction 40–15 Ma ago, then its average subduction rate is 20–60 mm yr^–1 from our tomographic result. This value is roughly consistent with the estimates (28–50 mm yr^–1) by the above-mentioned studies using other approaches. The convergence rate may have varied in different periods of the PHS Plate subduction (Uyeda & Miyashiro 1974; Hirahara 1981).

4.2 Low-velocity zones in the mantle wedge

Many previous tomographic studies have revealed a low-velocity (low-V) zone in the central part of the mantle wedge under the land area of NE Japan and suggested that the low-V zone represents the hot upwelling flow associated with the slab dehydration and corner flow in the mantle wedge and it forms the source of arc magmatism and volcanism (e.g. Zhao et al. 1994; Iwamori & Zhao 2000; Hasegawa et al. 2005, 2009; Huang et al. 2011; Tong et al. 2012). The present results show that the low-V zone in the mantle wedge extends westwards to the eastern margin of the Japan Sea (Fig. 8), not just confined beneath Honshu Island as suggested by the previous tomographic studies. The results of the synthetic test (Fig. S17b) indicate that this is a reliable feature. However, it is still unclear whether the mantle–wedge low-V zone extends further westwards under the entire Japan Sea or not.

Recent regional and global seismic studies show that the Pacific slab has subducted down to about ∼575km depth under the Japan Sea and then the slab becomes stagnant in the mantle transition zone under eastern China and western Japan (e.g. Zhao 2004; Huang & Zhao 2006; Zhao et al. 2011; Gu et al. 2012; Fig. 11). The western edge of the stagnant slab has reached Beijing, about 2000km west of the Japan Trench. It was suggested that a big mantle wedge (BMW) has formed above the stagnant slab (Zhao 2004) and the BMW exhibits low velocities under the Japan Sea and eastern China (Zhao et al. 2009, 2011; Zheng et al. 2011). Deep dehydration of the stagnant slab may take place and corner flow may exist in the BMW, which induced the upwelling of hot and wet asthenospheric materials and lead to the formation of the intraplate volcanoes (such as Changbai, Jeju and Ullung in Fig. 1) and lithospheric thinning in eastern China. The global and regional tomographic models still have a lower resolution, so it is unclear whether the BMW low-V zones are continuous from the volcanic front in Honshu to eastern China or they exist intermittently beneath Honshu, Japan Sea and eastern China. The spatial extent and detailed structure of the BMW low-V zones are an important issue related to the deep slab dehydration, mantle-wedge corner flow and origin of backarc and intraplate magmatism. To clarify this issue, a portable ocean-bottom-seismometer (OBS) network should be installed in the Japan Sea to determine a high-resolution upper-mantle structure under the Japan Sea and East Asia.

Beneath Kyushu, a significant low-V zone is visible in the crust and upper-mantle wedge above the PHS slab (Figs 10e–h). The low-V zone extends down to 300–400km depth dipping toward the west beneath the East China Sea, which is considered to be associated with the shallow and deep dehydration of the PHS slab and a convective circulation process in the mantle wedge, similar to the processes occurring in the mantle wedge above the Pacific slab (Fig. 8). These processes may have led to the formation and opening of the Okinawa Trough whose northern edge has reached the Unzen volcano in western Kyushu (e.g. Nakamura et al. 2003; Zhao et al. 2011; Fig. 1).

Significant low-V anomalies are imaged between the subducting PHS slab and the Pacific slab at depths of 100–500km under western Japan (Figs 9e–i and 10a–i). The low-V zones are connected with the Pacific slab, which may reflect fluids from the deep dehydration of the Pacific slab, as well as the convective circulation process in the mantle wedge. The Pacific slab is old (100–130 Ma) and has a large subduction rate (7–10 cm yr^–1) at the Japan Trench and so the slab has a lower temperature, which allows deep slab dehydration at 400–500km depth, similar to that occurring in the subducting Pacific slab under the Tonga arc and the Lau backarc spreading centre (Conder & Wiens 2006). Water transported to the deep upper mantle and the mantle transition zone is stored in some hydrous minerals such as phase A, phase E, superhydrous phase B, phase D and nominally anhydrous minerals such as wadsleyite and ringwoodite (Inoue et al. 2004; Komabayashi et al. 2004; Ohtani et al. 2004; Ohtani & Zhao 2009). Dehydration reactions proceed by decomposition of the hydrous minerals due to the temperature increase of the deep slab. The decrease of the maximum water solubility in wadsleyite and ringwoodite may also cause the deep slab dehydration (Inoue et al. 2004; Komabayashi et al. 2004; Ohtani et al. 2004).

The low-V zones at depths of 100–500km between the Pacific and PHS slabs should have a higher temperature due to the corner flow in the mantle wedge above the Pacific slab. Thus, the PHS slab would be heated from below and so its brittleness is lost at a shallower depth and so intermediate-depth earthquakes do not occur within the slab below 80-km depth under Chugoku and below 180-km depth under Kyushu. The heating by the hot upwelling flow from below may also reduce the stiffness of the PHS slab so that the PHS slab suddenly bends where the intraslab seismicity stops under western Honshu (Figs 10a–d).

4.3 Hot mantle upwelling under the Pacific slab

In the mantle below the Pacific slab under NE and central Japan, low-V anomalies are visible down to 700km depth, the bottom of the present tomographic model (Figs 8 and 9a–f). The result of the synthetic test indicates that this is a reliable feature (Figs S17b and S17d).

Being driven by the gravitational force, the subducting Pacific slab continues to sink down to the boundary between the upper and lower mantle, the 660km discontinuity, where the sudden increase in viscosity will not allow direct slab penetration into the lower mantle (e.g. Turcotte & Schubert 2001). Thus the Pacific slab becomes stagnant in the mantle transition zone and will finally collapse down to the lower mantle as a result of very large gravitational instability from phase transitions (Maruyama 1994; Zhao 2004; Maruyama et al. 2007). This would cause turbulences and thermal instability in the mantle transition zone and the lower mantle, and so mantle upwelling occurs as a consequence of the thermal instability at a thermal boundary layer (Albers & Christensen 1996). The thermal plume model (Morgan 1971) proposes that the thermal disturbance at the core-mantle boundary is responsible for the generation of mantle plumes. The low-V anomalies in the mantle below the Pacific slab imaged by this study (Figs 8 and 9a–e) may be a consequence of the thermal instability caused by the collapsing of the Pacific slab materials down to the lower mantle (Zhao 2001, 2004; Abdelwahed & Zhao 2007; Zhao et al. 2011).

The deep-subduction related hot mantle upwelling has been also revealed in other regions. For example, recent geophysical and geochemical studies show that a young mantle plume exists under the Hainan hotspot in southern China and the Hainan plume is caused by the thermal instability associated with the Indian plate's deep subduction in the west and the Pacific and PHS slabs' deep subduction in the east down to the lower mantle (Zhao 2004; Lei et al. 2009; Huang et al. 2010; Wang et al. 2012).

In this work and also the previous studies, only P-wave tomography is determined for the deeper part of the Japan subduction zone. In future studies, a large number of S-wave data should be collected from the teleseismic seismograms recorded by the Hi-net stations, so that S-wave tomography and Poisson's ratio images under the Japan islands can also be obtained. In addition, portable OBS stations should be deployed in the Pacific Ocean and the Japan Sea around the Japan islands to monitor the suboceanic seismicity and to record the teleseismic events, thus the structure and dynamics in the forearc and backarc regions will also be clarified. For a better interpretation of the tomographic images, information from other branches of Earth sciences, such as mineral physics and computer simulations, is needed, which may provide us new insights into the deep structure and dynamics of the Japan subduction zone.

5 Conclusions

We made great efforts to collect a large number of high-quality teleseismic data recorded by the dense seismic networks on the Japan islands. As a result, a high-resolution 3-D P-wave velocity model down to 700km depth under the Japan islands is determined by inverting the local and teleseismic data simultaneously. The novel tomographic model reveals the following new features.

• 1

  The PHS Plate has subducted aseismically down to 430km depth under SW Japan.

• 2

  The low-V zone in the mantle wedge under NE Japan and Kyushu exists not only under the land area but also extends westward to the backarc region under the Japan Sea and East China Sea.

• 3

  Significant low-V anomalies are revealed in the deep portion of the mantle wedge (400–500km depth) above the Pacific slab under central and western Japan, which may reflect the hot and wet mantle upwelling caused by mantle–wedge corner flow and fluids from deep dehydration of the Pacific slab.

• 4

  At depths of 420–700km beneath the Pacific slab under eastern Japan, low-V anomalies are visible, which may reflect upwelling of hot mantle materials from the lower mantle, induced by the deep subduction of the Pacific slab and its collapsing down to the core–mantle boundary.

Acknowledgments

The authors wish to thank the staffs of National Research Institute for Earth Sciences and Disaster Prevention and the Hi-net data centre for maintaining the seismic network and providing the high-quality data used in this study. Figures in this paper were produced using the GMT software (Wessel & Smith 1995). This work was partially supported by research grants (Kiban-A 17204037; Kiban-S 11050123) to D. Zhao from Japan Society for the Promotion of Science, as well as from the Global-COE program of Earth and Planetary Sciences of Tohoku University. We are very grateful to Dr. J. Nakajima who kindly provided the digital files for the shallow parts of the PHS slab in Fig. 11 and to Dr. Zhouchuan Huang for his help in data processing and computer graphics. Prof. Jean Virieux (the Editor) and two anonymous referees provided thoughtful review comments that have improved the manuscript.

Supporting Information

Additional Supporting Information may be found in the online version of this article:

Figure S1. Grid configuration adopted in this study. The horizontal grid interval is 0.33°. The vertical grid interval is 5-20 km for depths <160 km and 30 km for depths of 160-700 km.

Figure S2. The seismic velocity model adopted for the tomographic inversions. The subducting Pacific slab is considered in the starting velocity model. The depth variations of the Conrad and Moho discontinuities (Zhao et al. 1994) are also taken into account.

Figure S3. Depth contours to the upper boundary of the subducting Pacific slab (Hasegawa et al. 2009). The dashed lines show the major plate boundaries. Red triangles denote the active and quaternary volcanoes. The orange contour lines show the rupture areas of the large earthquakes in the megathrust zone under the northern Japan forearc.

Figure S4. Distribution of ray paths in different depth ranges under the study region.

Figure S5. Results of a checkerboard resolution test. The horizontal grid interval is 0.33°. The open and solid circles denote low and high velocities, respectively. The velocity perturbation (in per cent) scale is shown at the bottom.

Figure S6. The same as Fig. S5 but the horizontal grid spacing is 0.5°.

Figure S7. The same as Fig. S5 but the horizontal grid spacing is 1.0°.

Figure S8. The same as Fig. 4 but the P-wave tomography is obtained by adopting a homogeneous 1-D starting velocity model that does not include the subducting Pacific slab. Note that for the arrival-time data of local P waves from the intermediate-depth and deep earthquakes that travel longer within the subducting Pacific slab, their traveltime residuals amount to -2 to -3 s or even greater, which can be explained well when the Pacific slab is considered in the starting velocity model (Figs S2 and S3). Thus even for the same data set, the traveltime residuals are different in the homogeneous 1-D starting velocity model and in the inhomogeneous starting model that contains the subducting Pacific slab (Figs S2 and S3). To stabilize the inversion process, we have used only the local and teleseismic data that have traveltime residuals less than 2.0 s. Thus the actual number of data used becomes different in the two inversions with different starting models. Some data with larger residuals (-2 to -3 s or even larger) due to their rays travelling longer within the Pacific slab were not used in the inversion with the homogeneous 1-D starting model. These differences in the data number and traveltime residuals may have also caused the differences in the tomographic images resulted from the two inversions (Figs 4-10 and S8-S14).

Figure S9. The same as Fig. S8 but for the images at 190-430 km depths.

Figure S10. The same as Fig. S8 but for the images at 460-700 km depths.

Figure S11. The same as Fig. 7 but the P-wave tomography is obtained by adopting a homogeneous 1-D starting velocity model that does not include the subducting Pacific slab.

Figure S12. The same as Fig. S11 but for vertical cross-sections beneath Tohoku as shown on the inset map.

Figure S13. The same as Fig. S11 but for different vertical cross-sections as shown on the inset map.

Figure S14. The same as Fig. S11 but for vertical cross-sections beneath western Japan as shown on the inset map.

Figure S15. The input model of a synthetic test conducted for confirming the aseismic subduction of the PHS Plate under western Japan. The locations of the vertical cross-sections are shown on the inset map. The velocity perturbation (in per cent) scale is shown beside the map.

Figure S16. The same as Fig. S15 but the inversion results of the synthetic test. The results indicate that aseismic subduction of the PHS Plate under western Japan is a reliable feature.

Figure S17. The input model (left) and the inversion result (right) of a synthetic test conducted for confirming the main features revealed by this study. The shape of the subducting Pacific slab was fixed in the inversion process. The vertical cross-sections are along the profiles in Hokkaido (a), NE Honshu (b), Shikoku and Chugoku (c) and Kyushu and SW Japan (d), as shown on the inset map.

Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

Supporting Information

Additional Supporting Information may be found in the online version of this article:

Figure S1. Grid configuration adopted in this study. The horizontal grid interval is 0.33°. The vertical grid interval is 5-20 km for depths <160 km and 30 km for depths of 160-700 km.

Figure S2. The seismic velocity model adopted for the tomographic inversions. The subducting Pacific slab is considered in the starting velocity model. The depth variations of the Conrad and Moho discontinuities (Zhao et al. 1994) are also taken into account.

Figure S3. Depth contours to the upper boundary of the subducting Pacific slab (Hasegawa et al. 2009). The dashed lines show the major plate boundaries. Red triangles denote the active and quaternary volcanoes. The orange contour lines show the rupture areas of the large earthquakes in the megathrust zone under the northern Japan forearc.

Figure S4. Distribution of ray paths in different depth ranges under the study region.

Figure S5. Results of a checkerboard resolution test. The horizontal grid interval is 0.33°. The open and solid circles denote low and high velocities, respectively. The velocity perturbation (in per cent) scale is shown at the bottom.

Figure S6. The same as Fig. S5 but the horizontal grid spacing is 0.5°.

Figure S7. The same as Fig. S5 but the horizontal grid spacing is 1.0°.

Figure S8. The same as Fig. 4 but the P-wave tomography is obtained by adopting a homogeneous 1-D starting velocity model that does not include the subducting Pacific slab. Note that for the arrival-time data of local P waves from the intermediate-depth and deep earthquakes that travel longer within the subducting Pacific slab, their traveltime residuals amount to -2 to -3 s or even greater, which can be explained well when the Pacific slab is considered in the starting velocity model (Figs S2 and S3). Thus even for the same data set, the traveltime residuals are different in the homogeneous 1-D starting velocity model and in the inhomogeneous starting model that contains the subducting Pacific slab (Figs S2 and S3). To stabilize the inversion process, we have used only the local and teleseismic data that have traveltime residuals less than 2.0 s. Thus the actual number of data used becomes different in the two inversions with different starting models. Some data with larger residuals (-2 to -3 s or even larger) due to their rays travelling longer within the Pacific slab were not used in the inversion with the homogeneous 1-D starting model. These differences in the data number and traveltime residuals may have also caused the differences in the tomographic images resulted from the two inversions (Figs 4-10 and S8-S14).

Figure S9. The same as Fig. S8 but for the images at 190-430 km depths.

Figure S10. The same as Fig. S8 but for the images at 460-700 km depths.

Figure S11. The same as Fig. 7 but the P-wave tomography is obtained by adopting a homogeneous 1-D starting velocity model that does not include the subducting Pacific slab.

Figure S12. The same as Fig. S11 but for vertical cross-sections beneath Tohoku as shown on the inset map.

Figure S13. The same as Fig. S11 but for different vertical cross-sections as shown on the inset map.

Figure S14. The same as Fig. S11 but for vertical cross-sections beneath western Japan as shown on the inset map.

Figure S15. The input model of a synthetic test conducted for confirming the aseismic subduction of the PHS Plate under western Japan. The locations of the vertical cross-sections are shown on the inset map. The velocity perturbation (in per cent) scale is shown beside the map.

Figure S16. The same as Fig. S15 but the inversion results of the synthetic test. The results indicate that aseismic subduction of the PHS Plate under western Japan is a reliable feature.

Figure S17. The input model (left) and the inversion result (right) of a synthetic test conducted for confirming the main features revealed by this study. The shape of the subducting Pacific slab was fixed in the inversion process. The vertical cross-sections are along the profiles in Hokkaido (a), NE Honshu (b), Shikoku and Chugoku (c) and Kyushu and SW Japan (d), as shown on the inset map.

Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

References