The nature of the Pacific plate as subduction inputs to the northeastern Japan arc and its implication for subduction zone processes

Abstract

Devastating megathrust earthquakes and slow earthquakes both occur along subducting plate interfaces. These interplate seismic activities are strongly dependent on the nature of the plate interface, such as the shape of the plate interface and the materials and physical conditions along the plate interface. The oceanic plate, which is the input to the subduction zone, is the first order control on the nature of the plate interface. To reveal the nature of the subduction inputs to the northeastern Japan arc, we have conducted large-scale controlled-source seismic surveys of the northwestern part of the oceanic Pacific plate. The obtained seismic data have revealed (1) oceanic plate structural evolution caused by plate bending prior to subduction, suggesting the promotion of the oceanic plate hydration; (2) spatial variation of the oceanic plate structure, such as variations in the thickness of sediment and crust; (3) that the spatial variations are caused by both ancient plate formation processes and more recent volcanic activities; and (4) that spatial variations of the nature of the subduction inputs show good correlation with the along-strike variations in the seismic structure and seismic activities after subduction, including the coseismic slip distribution of the 2011 Tohoku earthquakes and the structural differences between the northern and the southern Japan Trench. These observations indicate that the incoming oceanic plate structure is much more spatially variable than previously thought and also imply that the spatial variation of the subduction inputs is a key controlling factor of the spatial variation of various processes in subduction zones, including interplate seismic activities and evolution of the forearc structure.

1 Introduction

In plate subduction zones, an oceanic plate is subducting beneath a continental plate, and many fault slip behaviors, ranging from devastating megathrust earthquakes to slow earthquakes, occur along the plate interface (e.g., Obara and Kato 2016). The nature of the plate interface, such as its geometry and the materials along the plate interface, is thought to control the occurrence of the various types of interplate earthquakes. Because both megathrust and slow earthquakes occur in the forearc region beneath the ocean, detailed seismic structures beneath the forearc region, especially at around the depth of the plate interface, are key to understanding the controlling factors of a wide variety of fault slip behaviors.

Marine controlled-source seismic surveys are effective tools for revealing the detailed seismic structure beneath the seafloor. In the forearc region of the Japan Trench, controlled-source seismic surveys have been conducted for more than 50 years. Initially, because of limitations in the performance of the seismic survey system and processing techniques, it was difficult to evaluate the plate interface at depths greater than about 10 km, where most large interplate earthquakes occur (e.g., Ludwig et al. 1966; Asano et al. 1981; Suyehiro et al. 1985; Suyehiro and Nishizawa 1994; von Huene et al. 1994). Advances in the seismic survey system and data processing techniques, e.g., use of a more powerful controlled seismic source, longer hydrophone streamer cables, and a denser array of ocean bottom seismometers (OBSs), as well as improved tomographic inversion techniques, have made it possible to model deeper structures than before. As a result, the plate interface has been modeled at a depth of up to 20–30 km along many survey lines perpendicular to the Japan trench (e.g., Tsuru et al. 2000, 2002; Miura et al. 2003; Takahashi et al. 2004; Ito et al. 2004, 2005; Fujie et al. 2006; Mochizuki et al. 2008). These results have enabled discussion of, for example, the relationship between the updip and downdip limits of the rupture zone of large interplate earthquakes and seismic structural features, such as changes in the plate dip angle and the depth of the forearc Moho.

Along-strike variation is a notable feature of interplate seismic activities in the Japan Trench. First, there is a remarkable difference in the distribution of large interplate earthquakes between the northern and the southern parts. Large interplate earthquakes $(M >7)$ have repeatedly occurred in the northern Japan Trench, but few large interplate earthquakes have been observed in the southern Japan Trench (e.g., Kawakatsu and Seno 1983; Yamanaka and Kikuchi 2004). The boundary between the northern and southern parts is located at around $37.5^\circ$N, where the strike of the Japan Trench changes. Tsuru et al. (2002) compiled many seismic reflection profiles perpendicular to the strike of the Japan Trench and found that there was a channel-like sediment unit in between the overriding plate and the subducting plate only in the southern Japan Trench; they proposed that the existence of the channel-like sediment unit along the plate interface works as an interplate stress releaser in the southern Japan Trench. Nakata et al. (2021) confirmed this idea in a numerical simulation.

Even within the northern Japan Trench, where large interplate earthquakes have repeatedly occurred, there are several seismically inactive regions where no large interplate earthquakes have been observed (hereafter, referred to as seismic gaps) (e.g., Yamanaka and Kikuchi 2004; Kawakatsu and Seno 1983). To reveal structural differences between the seismic gaps and rupture zones of large interplate earthquakes, Hayakawa et al. (2002) and Fujie et al. (2013a) conducted along-strike controlled-source seismic surveys. They found that compared to the rupture zones of large interplate earthquakes, the seismic structure of the upper plate in the seismic gaps shows thinner sediment, deeper forearc Moho, and high seismic attenuation, although the causes for these relationships were not determined.

Fujie et al. (2002) and Mochizuki et al. (2005) focused on seismic gaps where not only large earthquakes but also microearthquakes have rarely been observed (Fig. 1a, black dashed curves, hereafter, referred as aseismic areas). They conducted wide-angle seismic surveys across these aseismic areas and found that the plate interface in these areas was highly reflective. On the basis of this observation, they proposed that a thin layer with low seismic velocity related to fluid or hydrated minerals existed along the plate interface in the aseismic areas and that dehydration of the subducting oceanic plate or the serpentinized forearc mantle was a possible source of the fluid supply to the plate interface. These seismic structure studies showed good correlations between the along-strike variations of the seismic structure and seismic activities in the forearc region. However, the factors causing the along-strike variations in the forearc seismic structure remain poorly understood. Note that microseismicity has been very low after the 2011 Tohoku earthquakes and slow earthquakes have been observed in the aseismic areas, suggesting relatively weak plate coupling there (Nishikawa et al. 2019, 2023), although large coseismic slips were observed during the 2011 Tohoku earthquake there (e.g., Iinuma et al. 2012).

Fig. 1

a Seismic survey lines in the northwestern part of the Pacific plate. Red lines show seismic survey lines along which we shot the large-volume tuned airgun arrays of R/V Kairei and R/V Kaimei while towing a 6-km-long multi-channel hydrophone streamer cable to obtain multi-channel seismic (MCS) reflection data. White circles represent Ocean Bottom Seismometers (OBSs) aligned at a spacing of 6 km. The green contours show the slip distribution of the 2011 Tohoku earthquake derived from both terrestrial GPS and seafloor geodetic data (Iinuma et al. 2012); the contour interval is 10 m, and the areas with values $>30$ m are shaded in green. Dashed black curves outline regions of low microseismicity, where a thin low-velocity layer was identified between the overlying and subducting plates (Mochizuki et al. 2005). Yellow dashed lines represent discontinuities of the magnetic anomaly lineations (Nakanishi 2011). The thick black arrow shows the plate convergence direction. b Slope of the seafloor bathymetry of the oceanic plate, showing the horst and graben structure along the Japan Trench and the Kuril Trench. Red dashed lines are magnetic anomaly lineations (Nakanishi 2011). The blue dashed curve shows the Joban seamount chain

To explain the along-strike variations of the interplate earthquakes, Tanioka et al. (1997) studied the nature of the subduction inputs. They found that the seafloor roughness of the incoming oceanic plate correlates well with the distribution of large interplate earthquakes after subduction in the northern Japan Trench; smooth seafloor correlates with large earthquakes at the deeper part of the seismogenic plate interface and rough seafloor correlates with tsunami earthquakes occurring at the shallower plate interface. They proposed that interplate seismic coupling after subduction is affected by the seafloor roughness of the incoming plate because a smooth seafloor topography would be expected to result in a uniform frictional property along the plate interface, whereas a rough seafloor topography should make the frictional property heterogeneous. The seafloor roughness of the incoming plate is mainly dependent on the development of a horst and graben structure, which is formed by plate-bending-related normal faulting (hereafter, bend faulting) in the outer trench area (e.g., Masson 1991). Thus, their observation suggested that the along-strike variations in the bend faulting may be key in controlling along-strike variations in the interplate seismic activities after subduction.

The throw of a horst and graben can be 500 or more meters in the vicinity of the trench axis (e.g., Nakanishi 2011), suggesting that the offset is too large to have slipped all at once and slip is a drawn out process. The fault plane of large bend faulting can extend from the seafloor to a depth of more than 40 km (e.g., Hino et al. 2009; Kanamori 1971; Obana et al. 2012), implying that the oceanic plate is damaged to a depth of more than 40 km by repeated bend faulting. In the eastern Pacific subduction zones, pioneering work has revealed that bend faults penetrate from the seafloor into the topmost mantle (Ranero et al. 2003) and that P-wave velocities ($V_p$) within the oceanic crust and the topmost mantle become lower toward the trench, implying fracturing within the crust and the topmost mantle caused by bend faulting and seawater penetration and hydration through the bend faults (e.g., Ivandic et al. 2008; Grevemeyer et al. 2007).

If bend faulting actually promotes hydration of the oceanic plate, the amount of water transported by the subducting oceanic plate as a part of the hydrated minerals might be notably larger than previously assumed (Peacock 2001). The dehydration process itself and the water expelled from a hydrated subducting oceanic plate affect various processes in subduction zones, such as dehydration-induced intermediate-depth earthquakes and arc magmatism (e.g., Tatsumi 1989; Meade and Jeanloz 1991; Kirby 1995; Hyndman and Peacock 2003; Rüpke et al. 2004). In addition, water expelled from the subducting plate might ascend along the plate interface and form a layer with low seismic velocity (Fujie et al. 2002; Mochizuki et al. 2005), which could promote forearc mantle serpentinization and affect interplate seismic coupling. Therefore, to understand various processes occurring in the subduction zone, it is important to reveal the nature of the subduction inputs, especially the degree of oceanic plate hydration and its spatial variations.

To reveal the seismic structure of the northwestern part of the oceanic Pacific plate, Shimamura et al. (1983) conducted a controlled-source seismic survey in 1980. Their primary survey line was 1800 km long, and they modeled the 1-D $V_p$ structure of the oceanic plate based on traveltimes from the extremely long profile. This model was a kind of standard, averaged velocity model of the northwestern part of the Pacific plate. No special effort was made thereafter to reveal the detailed seismic velocity structure or its regional variations in the northwestern part of the Pacific plate. Therefore, the seismic structure of the northwestern part of the Pacific plate was regarded to be approximately homogeneous laterally, and the structural impacts caused by bend faulting were not studied for several decades.

In 2009, the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) started a campaign of extensive controlled-source seismic surveys to reveal the nature of the subduction inputs to the northeastern Japan arc. To reveal the normal (intact) oceanic plate structure and its evolution due to plate bending, the survey area has covered the area from the trench to 500 km east of the trench, where plate bending is not expected to affect the oceanic plate structure (Fig. 1a). Parts of this seismic survey campaign were conducted in the collaboration with the GEOMAR-Helmholtz-Centre for Ocean Research Kiel, Earthquake Research Institute (University of Tokyo), and Research Center for Prediction of Earthquakes and Volcanic Eruptions (Tohoku University). The large amount of seismic data obtained through this survey campaign confirmed systematic structural evolution due to bend faulting prior to subduction and showed that the oceanic plate, previously thought to be nearly homogeneous from the sediments to the topmost mantle, exhibits significant spatial variations. The data discussed in the Section 4 of this paper are seismic reflection and refraction profiles acquired in 2017 as part of this campaign, aiming to reveal the differences between the northern and the southern Japan Trench.

In this paper, we first briefly summarize the findings of our previous studies, which focused on primarily the central and northern Japan Trench (lines A2, A3, A4 and R2, Fig. 1a) to reveal (1) the structural evolution of the oceanic plate due to bend faulting, (2) regional variations in bend faulting and its implications, and (3) along-trench variations in the oceanic plate and its potential impacts on interplate seismic coupling. We then present new results in detail from the first wide-angle seismic survey of the southern Japan Trench (line A6), to compare with the previous results and expand the study of incoming plate structure farther south. Finally, we present a map of the thickness of the oceanic crust derived from all of the data and discuss the spatial variations in thickness and new implications for interplate seismic coupling and structural development of the subduction zone.

2 Data and method

2.1 Data acquisition

To reveal the seismic structure of the oceanic plate and its spatial variation in the northwestern part of the Pacific plate, we have conducted controlled-source marine seismic surveys that consist of both multi-channel seismic (MCS) reflection surveys and wide-angle seismic reflection and refraction surveys. As the seismic vessel, we have mainly used R/V Kairei and sometimes R/V Kaimei, both operated by JAMSTEC. At the time of the surveys, R/V Kairei was equipped with a large-volume (7,800 cubic inch) tuned airgun array and a 440-ch, 5500 m long hydrophone streamer cable with a channel interval of 12.5 m. R/V Kaimei was equipped with a 10,800 cubic inch tuned airgun array and an adjustable length hydrophone streamer cable with a channel interval of 3.125 m. During the project, the cable length was between 5500 (1760-ch) and 8700 m (2780-ch). For the MCS data acquisition (Fig. 1a, red lines without white circles), we fired the tuned airgun arrays at a spacing of 50 m.

To obtain the wide-angle seismic reflection and refraction data, we first deployed OBSs at a spacing of 6 km (Fig. 1a, white circles) and then fired tuned airgun arrays at a spacing of 200 m. Each OBS was equipped with a three-component geophone and a hydrophone, except for some ultra-deep OBSs deployed in water depths greater than 6000 m (lines A3, A4, and A6). While shooting for OBSs, we towed the hydrophone streamer cable to obtain MCS data even though the shooting spacing was 200 m. In addition, we separately obtained MCS reflection data with a shot spacing of 50 m along lines A2, A3, R2, and P2, as well as the western part of A4.

2.2 Methods

All of the MCS reflection data were first processed by post-stack time migration after a standard pre-processing flow including band-pass filtering, predictive deconvolution, velocity analysis, normal moveout, and stacking. Some of the MCS data were then reprocessed by using pre-stack time and pre-stack depth migration (e.g., Kodaira et al. 2017). To discuss the spatial variations of the oceanic plate structure using all of the profiles, we use time-migrated seismic sections instead of the depth sections. Lines A2, A3, and A4 and all of the along-strike lines were processed only by post-stack migration, but the other lines were processed by both pre-stack and post-stack time migration. We used pre-stack time-migrated sections preferentially when they existed.

Along all of the OBS survey lines, we developed P-wave velocity models by traveltime inversion using first arrivals and wide-angle Moho reflections picked from the OBS receiver gathers and two-way reflection times from the basement and Moho picked up on time-migrated MCS reflection sections (e.g., Fujie et al. 2013b, 2016). In addition, along lines A2 and A3, where we could stably observe P-to-S converted phases that were converted at the boundary between the sediment and the top of the oceanic crust, we also utilized the traveltimes of P-to-S converted phases to estimate the S-wave velocity structure (Fujie et al. 2013b, 2018).

The inversion results are generally dependent on the starting model. To evaluate its dependency and show the uncertainty of the obtained model, we applied Monte Carlo-type uncertainty analysis (e.g., Fujie et al. 2018, 2020; Korenaga et al. 2000). In addition, we applied checkerboard resolution tests (CRT) to evaluate the spatial resolution.

3 Summary of previous results

3.1 Structural evolution of the incoming oceanic plate by bend faulting

The seafloor depth of an oceanic plate gradually deepens as the plate’s age increases, but the surface of old oceanic plates such as the northwestern part of the Pacific plate is nearly flat (e.g., Stein and Stein 1992; Parsons and Sclater 1977). In order for the oceanic plate to subduct beneath the continental plate, the seafloor depth needs to increase rapidly toward the trench, causing the oceanic plate to bend downward (e.g., Walcott 1970; Watts and Talwani 1974). The plate bending causes an extensional stress field at the surface and a compressional stress field at depth (e.g., Christensen and Ruff 1983; Levitt and Sandwell 1995), leading to normal fault earthquakes (bend faulting) at the surface of the oceanic plate and thrust faulting at depth in the outer trench slope (e.g., Masson 1991; Gamage et al. 2009). The normal fault earthquakes generated at the bend faults, generally called outer rise earthquakes, occasionally are larger than magnitude 8 and produce devastating tsunamis, such as the 1933 Sanriku earthquake (Mw8.4) (e.g., Kanamori 1971; Chapple and Forsyth 1979; Uchida et al. 2016).

In the northwestern part of the Pacific plate, a horst and graben structure running roughly parallel to the trench axis, which is considered to have been formed by repeated bend faulting, is clearly observed on the bathymetric map and seismic reflection profiles (Figs. 1b and 2) in the area between the trench and 80 km seaward (e.g., Nakanishi 2011). In this region, Fujie et al. (2013b, 2018) modeled $V_p$ structure using controlled-source seismic survey data and showed a reduction of $V_p$ within the oceanic crust and the topmost mantle toward the trench. $V_p$ begins to decrease at more than 100 km from the trench, but it is initially confined to the upper oceanic crust, coinciding with the fracturing at the sediment basement recognized on the seismic reflection profile $V_p$ within the lower oceanic crust and the topmost oceanic mantle begins to decrease at $\sim 80$ km from the trench, which correlates well with the development of the horst and graben structure at the seafloor. These observations suggest that, at first, minor fracturing occurs at the top of the oceanic plate more than 100 km away from the trench. Large normal faults that cut across the entire oceanic crust and reach into the topmost oceanic mantle then start to develop $\sim 80$ km from the trench accompanied by the development of the horst and graben structure, resulting in a considerable $V_p$ reduction to a depth of the topmost mantle in the northwestern part of the Pacific plate.

Fig. 2

Time-migrated seismic reflection profiles in the outer trench area of lines A2, A3, and A4 showing horst and graben structure developed by faulting related to plate bending (Fujie et al. 2013b, 2018, 2020)

A $V_p$ reduction in the outer trench area has been observed in many subduction zones around the world, for example, in the Chile Trench (e.g., Contreras-Reyes et al. 2007, 2008b, 2008), the Middle American Trench (e.g., Grevemeyer et al. 2007; Ivandic et al. 2008, 2010; Van Avendonk et al. 2011), the Alaska Trench (e.g. Shillington et al. 2015), and the Tonga Trench (e.g., Contreras-Reyes et al. 2011). In most of these studies, two mechanisms have been proposed to explain the reduction of $V_p$: (1) an increase of porosity owing to fracturing caused by bend faulting and (2) seawater penetration along the faults and promotion of hydration (e.g., Van Avendonk et al. 2011). To distinguish these two mechanisms, S-wave velocity ($V_s$) is critically important because $V_s$ is much more sensitive to the presence of fluid and the lithology than $V_p$ (e.g., Christensen 1996; Takei 2002). However, $V_s$ modeling is not straightforward by marine controlled-source seismic surveys because the controlled source in the water cannot emit S-waves. To model $V_s$, it is necessary to utilize P-to-S converted phases, but P-to-S conversion is generally not stably observed and difficult to utilize. Therefore, most of the studies cited above using controlled-source seismic survey data modeled only $V_p$ structure and could not unequivocally distinguish between the above two possible explanations.

In the Kuril Trench and the Japan Trench, Fujie et al. (2013b, 2018) succeeded in determining $V_s$ within the upper oceanic crust by using controlled-source seismic survey data (Fig. 3b,e). The determination of $V_s$ was successful because (1) all OBSs were equipped with a gimbal-mounted three-component geophone that recorded high-quality horizontal components of the seismic signals, and (2) the basement (top of the igneous oceanic crust) worked as an efficient and stable P-to-S conversion interface, and no other efficient conversion interface existed. The obtained $V_p$ and $V_s$ models suggest that the $V_p/V_s$ ratio becomes large toward the trench axis, accompanied by the development of the bend faults, implying that the water content increases toward the trench axis at least up to the depth of the upper oceanic crust. This observation supports the idea that bend faulting promotes water penetration into the oceanic plate and possible hydration of the oceanic plate.

Fig. 3

Seismic velocity models along lines A2 and A3 modeled by traveltime inversion (Fujie et al. 2018). All of the seismic velocity models are the average of successful inversion results from 400 randomly generated starting models. a–c $V_p$ and $V_p$ perturbation within the mantle, and the $V_p/V_s$ ratio within the crust along line A2 (Kuril Trench). d–f $V_p$ and $V_p$ perturbation within the mantle, and the $V_p/V_s$ ratio within the crust along line A3 (Japan Trench). The inverted triangles show OBS locations. DbM is depth below the Moho, and Dbb is depth below the basement

The degree of mantle hydration is critically important because the mantle peridotite has the potential to contain much more water than crustal rocks by hydration (e.g., Peacock 2001). Unfortunately, Fujie et al. (2013b, 2018) could not constrain mantle $V_s$ because P-to-S converted phases, which propagate as S-waves in the mantle, were not stably observed because of weak signals. One way to constrain $V_s$ in the topmost oceanic mantle near the trench would be to utilize passive seismic data. Obana et al. (2018, 2019) presented tomographic $V_s$ models determined by using natural earthquakes that had occurred in the trench and outer trench area. Their results did not show significant changes in the $V_p/V_s$ ratio within the oceanic mantle even in the well-developed horst and graben area, where significant reduction of mantle $V_p$ was observed, implying that water penetration and oceanic plate hydration are mostly limited to the oceanic crust because hydration of mantle peridotite is expected to cause increase in $V_p/V_s$ ratio (Korenaga 2017; Christensen 2004). However, owing to limitations of the adopted 3-D tomographic inversion method, Obana et al. (2018, 2019) assumed an isotropic seismic structure, even though a strong mantle anisotropy (nearly 10%) was confirmed by controlled-source seismic survey data in this region (e.g., Kodaira et al. 2014). Although the influence of fracturing and hydration due to bend faulting on the seismic anisotropy is poorly understood, to discuss the degree of mantle hydration caused by bend faulting, we need to build more reliable seismic velocity models that can explain both the controlled-source and passive-source seismic data, including the seismic anisotropy.

3.2 Regional variations in bend faulting

As shown in Fig. 3, Fujie et al. (2018) observed systematic structural changes due to bend faulting, a reduction of $V_p$ and $V_s$, and an increase in the $V_p/V_s$ ratio toward the trench axis in both the Japan Trench (line A3) and the Kuril Trench (line A2). However, the degree of seismic structural changes was more remarkable in the Japan Trench than in the Kuril Trench, suggesting that the amount of water penetrating into the oceanic plate is larger in the Japan Trench.

The frequency of dehydration-induced intermediate-depth earthquakes occurring within the subducting oceanic plate after subduction (Hyndman and Peacock 2003; Rüpke et al. 2004) is often considered to be a good criterion for the degree of oceanic plate hydration (e.g., Shillington et al. 2015; Ranero et al. 2005). The frequency of intermediate-depth earthquakes is much higher in the Japan Trench subduction zone than that in the Kuril Trench subduction zone (Kita et al. 2006), suggesting that the degree of hydration of the subducting oceanic plate is higher in the Japan Trench, which is consistent with the seismic structure changes prior to subduction. Because the same oceanic plate of similar age is subducting into both the Japan Trench and the Kuril Trench, differences in the degree of hydration after subduction are considered to be caused by processes just prior to or after subduction. Therefore, Fujie et al. (2018) proposed that the degree of oceanic plate hydration is primarily determined by bend faulting just prior to subduction and it, in turn, determines the occurrence of dehydration-induced earthquake activity deep in the subduction zone.

Differences in bend faulting can be attributed to the difference in the type of bend faults. The strike of the bend faults tends to develop parallel to the trench axis because the axis of plate bending is basically parallel to the trench axis. If there is an existing weak zone subparallel to the trench axis, the plate-bending stress concentrates on the weak zone and the bend fault will likely develop along that zone. When the oceanic plate is formed at the oceanic spreading ridge, normal faults called abyssal-hill faults are created parallel to the ridge axis (e.g., Macdonald et al. 1996). Ancient abyssal-hill faults are considered to be weak zones even after a long time, and they are reactivated as bend faults in the outer trench area when the angle between the ancient spreading ridge axis and the current trench axis is small. Billen et al. (2007) studied bend faults around the world and reported that when the angle between the ancient spreading ridge axis and the current trench axis is smaller than $25^\circ$, abyssal-hill faults will be reactivated as bend faults, but new faults are formed when the angle is larger than $25^\circ$. In the northwestern part of the Pacific plate, the strike of the ancient spreading ridge axis suggested by magnetic anomaly lineations is subparallel to the Kuril Trench (with a trench-ridge angle of about $10^\circ$), but it is oblique to the Japan Trench (Fig. 1b). Thus, the bend faults in the Kuril Trench are reactivated abyssal-hill faults, and those in the Japan Trench are newly formed ones (e.g., Nakanishi 2011; Nakanishi et al. 1992a).

Because the overall flexure of the plate is about the same in both subduction zones, the extensional stress related to plate bending is also considered to be about the same. However, the spacing and the throw of the bend faults are significantly different between these trenches. In the Kuril Trench, the spacing of the bend faults is a few kilometers and is determined by the spacing of the ancient abyssal-hill faults (Macdonald et al. 1996), with a maximum throw of about 400 m. The spacing of the bend faults in the Japan Trench, however, is two to three times larger, and the throw is also much larger than in the Kuril Trench. Because the formation of new faults is considered to be more difficult than the reactivation of the existing weak zones, a much larger extensional stress is concentrated on each newly formed fault in the Japan Trench, resulting in much more developed bend faults that cut across the entire crust and reach into the oceanic mantle in the Japan Trench than that in the Kuril Trench (Fig. 2), which probably causes more mantle hydration in the Japan Trench.

In summary, Fujie et al. (2018) proposed that the difference in structural changes and hydration due to bend faulting between the Japan and Kuril Trenches can be explained by differences in the bend faults, which are determined by the angle between the current trench axis and the ancient spreading ridge axis.

Recently, Boneh et al. (2019) studied the relationship between the throw of bend faults prior to subduction and the activity of intermediate-depth earthquakes in subduction zones around the world and suggested that the hydration of the oceanic plate is generally controlled by plate bending and not by some inherent feature, which is consistent with the conclusion of Fujie et al. (2018).

3.3 Along-trench variations in the oceanic plate and their potential impacts on interplate seismic coupling

3.3.1 Ancient scar on the oceanic plate

To reveal the structural variations of the incoming oceanic plate along the Japan Trench, Fujie et al. (2016) conducted a wide-angle seismic reflection and refraction survey along line R2, located 120 km east of the northern Japan Trench (Fig. 1). The obtained $V_p$ and $V_s$ models showed that the central part of the survey line, at around $39^\circ$N, showed different structural features than the northern and southern parts. In the central part, the basement is shallow and the Moho is deep, meaning that the oceanic crust is thick in this area. In addition, $V_p$ and $V_s$ are low and $V_p/V_s$ is high as compared with the northern and southern parts (Fig. 4). Based on a magnetic anomaly lineation, Nakanishi (2011) found a plate age discontinuity at around $39^\circ$N that crosses perpendicular to line R2 (Fig. 1, yellow dashed line) and pointed out that this age discontinuity is a pseudofault associated with the ridge propagation generated near the spreading ridge when the plate was formed.

Fig. 4

a Pre-stack depth migration (PSDM) section along the northern part of line R2 (Fujie et al. 2016). PSDM was applied to the MCS reflection data obtained by hydrophone streamer cable, but the velocity model for PSDM was developed by traveltime inversion using OBS data. The Moho is relatively deeper and the basement is relatively shallower near the area of the pseudofault than in the adjacent area. The red arrow on the top shows the pseudofault, and the black shows the intersection with the line A3. b $V_p/V_s$ ratio along the northern part of line R2 derived by the tomographic inversion of OBS data (Fujie et al. 2016); $V_p/V_s$ is shown to be relatively high near the pseudofault. c A map around the survey line R2. The red line shows the survey line corresponding to a and b. Yellow dashed lines represent discontinuities of the magnetic anomaly lineations, and the western one is interpreted as a pseudofault (Nakanishi 2011)

Because the oceanic plate around the pseudofault is expected to be pervasively fractured (e.g., Hey 1977), the reduction of $V_p$ and $V_s$ and the increase in $V_p/V_s$ in the central part can be interpreted as the result of fracturing and water penetration related to the pseudofault. Therefore, the amount of water transported by the oceanic plate is considered to be larger in the central part of line R2 than that in the northern and the southern parts.

The landward extent of the pseudofault is difficult to constrain because magnetic anomaly lineations become obscured and untraceable after subduction. However, if we extrapolate the pseudofault linearly toward the subduction zone, it would reach the seismic gaps where even microearthquakes are rarely observed (Fig. 1, black dashed curves). As described above, Fujie et al. (2002) and Mochizuki et al. (2005) proposed that there is a thin layer with low $V_p$ along the plate interface in these seismic gaps and speculated that the formation of the thin low-velocity layer is related to dehydration processes occurring deep in the subduction zone. Based on the good agreement of the location between the pseudofault and the seismic gap, Fujie et al. (2016) suggested that the pseudofault, an ancient scar related to the ridge propagation generated at around the mid-ocean ridge, determines the amount of water transported into the subduction zone and affects interplate seismic coupling in this subduction zone.

As described above, Tanioka et al. (1997) pointed out that the seafloor roughness of the incoming plate correlated well with the distribution of large interplate earthquakes. The rough seafloor area pointed out by Tanioka et al. (1997) corresponds exactly to the central part of line R2. This good agreement suggests that the pseudofault might be responsible for the along-strike variations in the development of the horst and graben structure prior to subduction as well as for the distribution of large interplate earthquakes after subduction.

3.3.2 Recent volcanic activity on the oceanic plate

The surface topography and materials of the incoming oceanic plate are key factors determining the nature of the plate interface and controlling interplate slip behaviors. For example, surface topographic irregularities of the incoming oceanic plate, such as seamounts, cause spatial variations in the stress distribution after subduction and can affect the distribution of both ordinary and slow earthquakes (e.g., Kodaira et al. 2000; Mochizuki et al. 2008; Sun et al. 2020; Shiraishi et al. 2020). The thickness of the incoming sediment is considered to be one possible control on the maximum size of interplate earthquakes in each subduction zone (e.g., Ruff 1989; Heuret et al. 2012; Scholl et al. 2015). In addition, surface materials of the incoming plate are expected to directly affect the frictional properties along the plate interface after subduction and to play important roles in determining interplate seismic coupling.

One notable feature of the M9.0 2011 Tohoku earthquake was the large coseismic slip ($>50$ m) in the vicinity of the trench axis (e.g., Ide et al. 2011; Iinuma et al. 2012; Fujiwara et al. 2011; Kodaira et al. 2012). Immediately after the 2011 Tohoku earthquake, the Japan Trench Fast Drilling Project (JFAST) was carried out on Integrated Ocean Drilling Program Expedition 343 and 343T and recovered core samples from the fault plane (decollement) that hosted the large coseismic slip (e.g., Chester et al. 2013). The key material that enabled the large coseismic slip was smectite, a very weak clay mineral, which is abundant around the depth of the slip interface and has extremely low frictional properties (Ujiie et al. 2013). The scientific drilling results on the northwestern part of the oceanic Pacific plate show that the origin of the smectite at the plate interface is smectite-rich pelagic clay within the incoming oceanic sediments (Kameda et al. 2015), which are uniformly distributed on the northwestern part of the Pacific plate (Moore et al. 2015). In summary, the incoming sediments, containing smectite-rich pelagic clay, enabled the large coseismic slip of the 2011 Tohoku earthquake.

The sediment on the northwestern part of the Pacific plate has typically been considered to be 300–400 m thick and nearly uniform (Divins 2003), consisting of (from bottom to top) a chert layer, a thin smectite-rich pelagic clay layer, and a thick hemi-pelagic sediment layer (e.g., Moore et al. 2015; Shipboard Scientific Party 1980). However, newly obtained seismic data gathered during the seismic survey campaign, shown in Fig. 1, and another seismic survey campaign that used JAMSTEC’s portable MCS system in the vicinity of the trench axis (Nakamura et al. 2014, 2020; Qin et al. 2022) revealed remarkable spatial variations in the sediment thickness and clarified the existence of areas of extremely thin sediments, less than a few tens of meters (Fig. 5).

Fig. 5

a Seismic reflection profiles along lines A4 and R2. The inverted triangles show OBS positions. On the northwestern part of the Pacific plate, sediment thickness is generally 300–400 s in two-way time, corresponding to 300–400 m assuming the average $V_p$ within the sediments is 2.0 km/s, but we found areas of thin sediments such as the area from 260–380 km along this line. b Magnified seismic reflection section (left) and receiver function (RF) of controlled-source OBS data (right) in the area of thick sediments. One of characteristics of the sedimentary layer in this region is the existence of the chert unit at the bottom of the sediments, immediately above the basaltic crust (Moore et al. 2015; Shipboard Scientific Party 1980). The two P-to-S conversion interfaces imaged at about 2 s on the OBS RF section are interpreted as the top of the chert and the top of the basaltic crust, respectively. c The same figures as b but in the thin-sediment area. The characteristic appearance of the chert unit could not be recognized in the seismic reflection section, but multiple P-to-S conversion interfaces were observed in the OBS RF sections, suggesting that the strong reflection of the MCS reflection section is the apparent basement and multiple seismic structure boundaries exist immediately beneath the apparent basement (Fujie et al. 2020)

The distribution of the smectite-rich pelagic clay layer is important in relation to the large coseismic slip of the 2011 Tohoku earthquake. However, because this layer is difficult to recognize on the seismic reflection profiles, Fujie et al. (2020) mapped areas where they could clearly observe the characteristic appearance of the chert layer as an alternative to mapping the distribution of the smectite-rich pelagic clay layer (Fig. 6b). The distribution of the chert layer shows good agreement with spatial variations in the sediment thickness, suggesting that the smectite-rich pelagic clay layer is also missing in the thin-sediment areas. One of the thin-sediment areas that is missing chert occurs near the Japan Trench at $\sim 39^\circ$N (northwest thin-sediment and chert-missing area of Fig. 6). According to coseismic rupture distribution models derived from seismic and geodetic data (e.g., Lay 2018; Iinuma et al. 2012), the large coseismic slip of the 2011 Tohoku earthquake did not propagate beyond $39^\circ$N, suggesting that the absence of the smectite-rich pelagic clay layer might relate to the prevention of the large coseismic slip from propagating further northward.

Fig. 6

a Thickness of sediment in two-way traveltimes determined by MCS surveys (Fujiのはe et al. 2020). Red circles show petit-spot cluster sites where many small and young volcanoes are observed (Hirano et al. 2006). Petit-spot sites A and C are located at the areas of thin-sediment cover. b Distribution of chert unit along survey lines (red) and coseismic slip distribution of the 2011 Tohoku earthquake determined by geodetic study (black contour lines and gray shaded area) (Iinuma et al. 2012). Yellow circles represent the area of petit-spot clusters (e.g., Hirano et al. 2006; Hirano 2011), which corresponds to the area of missing chert (Fujie et al. 2020). Yellow squares represent tremors determined by S-net seafloor observatories (Nishikawa et al. 2019). c Schematic diagram of the sedimentary structure proposed by Fujie et al. (2020). The thick-sediment area is based on the MCS reflection profiles and drilling results (Moore et al. 2015; Shipboard Scientific Party 1980). The thin-sediment area is interpreted based on OBS RF and the existence of the nearby petit-spot clusters. The key point of the thin-sediment area model is that the bottom of the original sediment is heavily disturbed and metamorphosed by recent magmatic intrusions

To investigate the detailed seismic structure beneath the thin sediments, Fujie et al. (2020) applied a receiver function (RF) analysis to OBS-airgun data, which is a technique used to detect the P-to-S conversion interface (e.g., Vinnik 1977). To highlight the sediment-crust boundary, they chose an offset range of $9 \sim 25$ km, which corresponded to the crustal refractions. In the thick-sediment (normal) area, two P-to-S conversion interfaces were imaged at around 2.0 s, which were interpreted as the top of the chert and the top of the basaltic crust (Fig. 5). In contrast, in the thin-sediment area, multiple P-to-S conversion interfaces were detected between 0 and 2.0 s, suggesting that multiple structural interfaces exist between the seafloor and the top of the crust.

On the northwestern part of the Pacific plate, many young (1–10 Ma), small monogenetic volcanoes, called petit-spots, have been found in clusters (e.g., Hirano et al. 2006; Hirano 2011). The thin-sediment and chert-missing areas (Fig. 6b, yellow circles) correspond to the areas where the petit-spot clusters were identified by rock sampling and seafloor topography (e.g., Hirano et al. 2006; Hirano 2011). Because the alkalic basalt of the petit-spot volcanoes has extremely low viscosity, sills and dikes are considered to have frequently intruded into the sediments (Hirano et al. 2006). The multiple P-to-S conversion interfaces detected by OBS RF can be explained by these magmatic intrusions. Thus, Fujie et al. (2020) proposed that the magmatic intrusions and thermal metamorphism associated with petit-spot volcanism disturbed the incoming sediments, causing the transformation of smectite into illite (Fig. 6c), and that the subduction of such disturbed sediments prevented the northward propagation of the large coseismic slip of the 2011 Tohoku earthquake. Recently, Nishikawa et al. (2019) investigated various earthquake activities by using S-net seafloor observatories along the Japan Trench and found that various types of seismic activities, such as microearthquakes, tremors and very low frequency earthquakes (VLFEs), occurred near where the chert-missing area is subducting. In particular, tremors were concentrated in the extremely thin-sediment area (i.e., the northern half of the chert-missing area). The disturbed sediment acted as a barrier to the propagation of the large coseismic slip near the trench during the 2011 Tohoku earthquake, but it might also host various types of small earthquakes, especially tremors, in other circumstances.

4 Incoming plate of the southern Japan Trench

4.1 Wide-angle seismic survey in the outer trench area of the southern Japan Trench

In 2017 and 2018, in cooperation with JAMSTEC and the Earthquake Research Institute (University of Tokyo), we conducted a wide-angle seismic reflection and refraction survey along line A6, which runs normal to the Japan Trench and intersects it at $\sim 37^\circ$N (Figs. 1 and 7). This was our first wide-angle seismic survey on the incoming plate of the southern Japan Trench. Because of limitations in ship time and number of OBSs, we conducted the seismic survey along line A6 in two separate cruises. In 2017, we deployed OBSs between the western end and site S67 and collected wide-angle seismic reflection and refraction data onboard R/V Kairei. In 2018, we deployed OBSs from site S59 to the eastern end and collected seismic data onboard R/V Kaimei. To connect the two separate seismic survey datasets smoothly, we overlapped the area between two surveys and redeployed 10 OBSs in 2018 at the same sites as in 2017. We then reshot the line with the tuned airgun array so that we could obtain seismic data up to an offset distance of at least 100 km. Based on these overlapping OBS data, we confirmed the consistency of the traveltime data between two separate surveys. Therefore, we used both datasets to model $V_p$ structure by traveltime inversion.

Fig. 7

Map of line A6 and examples of OBS record sections. Note that the amplitude of first arrivals becomes weak at offsets larger than $\sim 20$ km at many eastern sites, such as S80 and S90, which is generally called the shadow zone and often observed in this part of the Pacific plate (e.g., Fujie et al. 2016, 2018). Such shadow zones suggest a very small or negative vertical velocity gradient within the oceanic crust (i.e., velocity reversal within the oceanic crust) (e.g., Zelt 1999). See the supplementary materials for more OBS record sections

4.2 $V_p$ structure modeling

The most notable bathymetric feature of line A6 is the Iwaki seamount, which is the largest seamount of the Joban seamount chain (Nakanishi 2011), located about 100 km east from the trench axis (Fig. 7). Although the seismic energy attenuated across the seamount, we observed high-quality P-wave refractions that had an offset distance of up to $\sim 100$ km at most sites (Figs. 7 and S1). In this study, we focused on the incoming oceanic plate and adopted the same $V_p$ modeling approach as before (Fujie et al. 2018) in order to directly compare models along lines A2, A3, A4, and A6; we divided seismic structure models into four layers (seawater, sediment, crust, and mantle) separated by three interfaces (seafloor, basement, and Moho) and determined $V_p$ by using first arrivals (refractions) and Moho reflections picked on OBS record sections, and two-way traveltimes from the basement and Moho picked on the time-migrated MCS reflection profiles.

To evaluate the dependency of the resultant model on the starting model, we applied a Monte Carlo-type uncertainty analysis (Korenaga et al. 2000; Sambridge and Mosegaard 2002; Fujie et al. 2018) as follows. First, we randomly constructed a wide range of starting models. The random parameters for the starting models were, $V_p$ at the top of the oceanic crust, the vertical $V_p$ gradient within the oceanic crust, thickness of the oceanic crust, $V_p$ at the top of the oceanic mantle, and the vertical $V_p$ gradient within the mantle. Note that sediment thickness and sediment $V_p$ are common to all of the starting models. The average and standard deviation of the starting $V_p$ models are shown in Fig. 8a and b. We then applied traveltime inversion from each randomly generated starting model and collected 400 final models that fulfilled the condition $\chi ^2 < 1$, where $\chi ^2$ is the normalized sum of the root mean square of traveltime misfits divided by the uncertainty of each traveltime pick. After that, we calculated the average and standard deviation of the 400 final $V_p$ models (Fig. 8c, d). In this study, to obtain a representative $V_p$ model, we applied one more traveltime inversion using the average $V_p$ model as the starting model. Figure 8e shows the representative $V_p$ model. To evaluate the spatial resolution of the representative $V_p$ model, we applied a CRT with a checkerboard pattern size of 30 $\times$ 3 km and a perturbation of 5%. The checkerboard patterns were well recovered within the upper crust and the topmost mantle (Fig. 8g).

Fig. 8

The results of $V_p$ modeling of traveltime inversion along line A6. a Average of the 400 starting $V_p$ models. The inverted triangles on the seafloor show the OBS positions. The yellow line is the average depth of the Moho, and the error bars along the Moho are the standard deviation of the Moho depth. b Standard deviation of the 400 starting $V_p$ models. c Average of the 400 final $V_p$ models. d Standard deviation of the 400 final $V_p$ models. e The representative $V_p$ model determined by traveltime inversion using the average inversion model in c as the starting model. The area with no ray-path coverage is shaded in gray. f Ray-paths of the wide-angle Moho reflection (PmP) on the representative $V_p$ model of e. g Recovery of the checkerboard resolution test (CRT). In this CRT, a 5% velocity perturbation was added to the representative $V_p$ model of e with a pattern size of 30 km in the horizontal direction and 3 km in the vertical direction. h The same $V_p$ model as e flattened by the depth of the basement (the top of the oceanic crust). The yellow line is the relative Moho depth from the basement, showing the variation in thickness of the oceanic crust along the profile. The red line is plotted as a reference depth of 6.5 km. Note that this figure is vertically exaggerated (2$\times$) relative to the other panels

4.3 Obtained model and its characteristics

The obtained oceanic plate structure model is basically similar to those in survey lines A3 (Fujie et al. 2018) and A4 (Fujie et al. 2020); the crustal thickness is about 6.5 km, except for in the seamount area, and the mantle $V_p$ becomes lower toward the trench. In the area 250–350 km from the trench axis, $V_p$ reversal with increasing depth is observed within the oceanic crust (Fig. 8c, e), which is similar to what is often observed in this part of the Pacific plate (e.g., Fujie et al. 2016, 2018) and supported by the “shadow zone” of the first arrivals observed on the OBS record section in this area (Figs. 7 and S1, S4).

The Iwaki seamount shows distinct features. Its height is about 4 km above the seafloor and its diameter of $\sim 60$ km. Because of the poor ray coverage of the Moho reflection (Fig. 8f), we could not adequately constrain the Moho below the seamount. However, if the $V_p$ contour lines, such as 7.0 and 7.5 km/s, can be considered as a proxy of the Moho there, the Moho below the seamount would be relatively deeper than in the surrounding area, which is consistent with other seamounts on the northwestern Pacific plate, such as the Daiichi–Kashima seamount of the Joban seamount chain (Nishizawa et al. 2009) and seamounts of the Marcus–Wake seamount chain (Kaneda et al. 2010). A cone-shaped intrusive high $V_p$ core was clearly imaged at the center of the Iwaki seamount body. Similarly high $V_p$ intrusive cores have been reported in many seamounts (e.g., Kaneda et al. 2010; Contreras-Reyes et al. 2010; Gallart et al. 1999). On the other hand, $V_p$ at a depth of about 10 km of the Iwaki seamount (about 6.4–6.5 km/s) is lower than the $V_p$ at the same depth of the surrounding normal oceanic crust (6.7–6.8 km/s). A similar relatively low $V_p$ layer ($\sim 6.5$ km/s) was modeled at the mid-depth of the Daiichi–Kashima seamount (Nishizawa et al. 2009), suggesting that relatively low $V_p$ at around the mid-depth of the crust might be a common feature of the Joban seamount chain. In summary, the seismic velocity structure of the Iwaki seamount is basically similar to other seamounts found on the Pacific plate.

In the northern Japan Trench, the impacts of bend faults become remarkable at about 80 km from the trench axis (Fujie et al. 2018, 2020). Along line A6, unfortunately, the western foot of the Iwaki seamount is located about 80 km from the trench axis, and distinguishing $V_p$ reduction caused by bend faulting from the inherent low $V_p$ at mid-depth of the seamount body is not a straightforward task. However, in the lower crust and topmost mantle in the vicinity of the trench (0–30 km from the trench axis), $V_p$ values are significantly lower than those 50 km from the trench. Therefore, we infer that, at least in the vicinity of the trench, $V_p$ values within the crust and topmost mantle become lower because of bend faulting. Although the Moho depth in the vicinity of the trench axis could not be well constrained owing to the lack of Moho reflection picks (Fig. 8f), it is natural to assume the crustal thickness is roughly constant on the normal oceanic plate without seamounts like this area, suggesting that the thickness of the oceanic crust in the vicinity of the trench is about 6.5 km, nearly the same as that estimated by the PmP reflections at around 40 km from the trench. Based on this assumption, $V_p$ at the topmost mantle became lower toward the trench and estimated to be $\sim 7.5$ km/s or less in the vicinity of the trench on line A6 (Fig. 8h), which is comparable to values found along other survey lines in the northern part of the Japan Trench (Fig. 3e). Therefore, we infer that $V_p$ reduction due to bend faulting is similar between the northern and the southern Japan Trench, although the altered area might be narrower along line A6 because of the existence of the large seamount.

Although we described the crustal thickness to be roughly 6.5 km above, the crust is relatively thinner in the eastern area of line A6 than in the western area (Fig. 8h); roughly speaking, the crust is thicker than 6.5 km in the western area but thinner in the eastern part. The same trend has been confirmed along lines A2, A3, and A4. In the next section, we discuss the thickness variations of the oceanic crust in this region based on both OBS $V_p$ structure and time-migrated MCS reflection profiles, focusing on a few hundred-km scale structural variations characterized by the thickness of the oceanic crust.

5 Thickness of the oceanic crust and its implications

5.1 Spatial variations in crustal thickness

To better understand spatial variations in the thickness of the oceanic crust, we first investigated the $V_p$ model derived from the wide-angle seismic reflection and refraction data along lines A2, A3 (Fujie et al. 2018), A4 (Fujie et al. 2020), and A6. To quantitatively compare the crustal thickness of all these lines quantitatively, we applied a Monte Carlo type uncertainty analysis to estimate the uncertainty in the crustal thickness determined by the traveltime inversion; for each line, we (1) applied traveltime inversions using randomly generated starting models and collected 400 successful $V_p$ models, (2) built a representative model by applying traveltime inversion using the average of the 400 successful $V_p$ models as the starting model, and (3) flattened the representative model with the basement depth. This is the same approach used to generate Fig. 8h along line A6. We adopted the results of the Monte Carlo-type uncertainty analysis of lines A2 and A3 from Fujie et al. (2018) and those of line A4 from Fujie et al. (2020).

Figure 9 shows the flattened $V_p$ models along A2, A3, A4, and A6, showing crustal thickness along each line. Roughly speaking, the crustal thickness is $\sim 6.5$ km (red lines in Fig. 9), except in the seamount areas. However, the crust becomes thinner toward the east along all lines. In addition, we found that the thickness of the crust shows a good correlation with the seafloor depth, except in the outer trench slope area; the thick crust areas mostly correspond to shallow seafloor depth areas and vice versa (Fig. 10), implying a relation to an isostatic equilibrium.

Fig. 9

$V_p$ models flattened by the basement (top of the crust) along lines A2, A3, A4, and A6, showing the along-profile variations in crustal thickness. The horizontal axis is the distance from the trench axis, and the vertical axis is the relative depth from the top of the crust. The A6 panel is the same as Fig. 8h. The yellow lines indicate the Moho depth, showing the crustal thickness along each line determined by the traveltime inversion. The red lines are plotted at a reference depth of 6.5 km to show variations in the thickness of the oceanic crust along each line and among all of the lines

Fig. 10

Seafloor bathymetry along lines A2, A3, A4, and A6. Except for the outer trench slope area (0–80 km) and large seamount areas, the seafloor depth shows the following trend; (1) in the southeast-side (350–500 km), water depth is roughly the predicted seafloor depth by the plate age-depth curve (e.g., Crosby et al. 2006; Müller et al. 2008), and (2) in the central part (100–250 km), the seafloor depth of the northern two lines (A2 and A6) is shallower than the predicted ones, but that of the southernmost line A6 is roughly equal to the predicted ones. The line A4 shows intermediate characteristics. This observation suggests that the thick crust area (central part of the line A2 and A3, Fig. 9) corresponds to the shallow water depth area

To evaluate the spatial variations in the thickness of the oceanic crust in more detail, we mapped the thickness of the oceanic crust using all the MCS reflection profiles shown in Fig. 1. As described previously, some of the MCS data were not processed in depth, so we used time-migrated profiles to map the crustal thickness. The processing flow was similar to that used to flatten the $V_p$ models: (1) Pick the basement (top of the crust), (2) flatten the time-migrated section using the basement picks, and (3) pick the Moho reflections on the flattened profiles (Fig. 11). Moho picks obtained in this manner are equivalent to the thickness of the oceanic crust in two-way time.

Fig. 11

a Time-migrated MCS section along line A6. The red marks indicate the picks of the basement (top of the crust) reflections. The reflection from the Moho is located at about 2 s after the basement reflections. b Time-migrated MCS section after flattening by the basement picks. Note that we could not apply flattening in the area with no basement picks. c Crustal thickness. Blue marks show Moho picks on the flattened time-migrated section. The dashed red line plotted at 2.0 s is a reference line to help illustrate variations in the crustal thickness along the profile

We applied this processing flow to all MCS reflection data (Figs. 12 and S2) and plotted the thickness of the oceanic crust in two-way time on the seafloor bathymetry map (Fig. 13). As expected, we confirmed a good correlation between the seafloor depth and the thickness of the oceanic crust. The thick crust, 2.1–2.3 s in two-way time, was observed on the outer trench area of the Kuril Trench where the seafloor depth is $\sim 5200$ m. In contrast, the thinnest crust, $\sim 1.5$ s in two-way time, was located in the eastern part of lines A4 and A6, where the seafloor is deepest in our survey area, $\sim 5800$ m, except for the trench area.

Fig. 12

Example of time-migrated MCS sections after flattening by the basement (top of the oceanic crust) reflection picks. The processing flow is the same as in Fig. 10, and the line positions are shown in Fig. 1. The blue line segments are the interpreted Moho, showing the thickness of the oceanic crust. The dashed red line plotted at 2.0 s is a reference line for variations in the crustal thickness along the profile. See the supplementary materials for more examples

Fig. 13

Map of the crustal thickness in two-way time derived from all the MCS lines shown in Fig. 1. The processing flow was the same as in Fig. 10. Note that the crustal thickness could be mapped only the area where we could identify both basement and Moho reflections. The seafloor topography was plotted in grayscale, and water depth contour lines were plotted as the thin black line every 200 m. The thick contour line corresponds to a 5400 m water depth, showing the approximate extent of the Hokkaido Rise (Kobayashi et al. 1998)

5.2 The Hokkaido Rise

The area where the seafloor depth is shallow and the oceanic crust is the thickest is a part of the Hokkaido Rise (e.g., Kobayashi et al. 1998) (Fig. 14). The Hokkaido Rise is a prominent topographic high located in the outer trench area of the Kuril Trench, and the peak around our survey area is as shallow as 5100 m, which is nearly 1000 m shallower than the oceanic basin located 500 km to the south (Figs. 1 and 13). Because the Hokkaido Rise is located just outside of the trench axis and elongated along the trench for more than 1200 km (Fig. 14), many previous studies have tried to explain this prominent topographic high by the formation of the outer rise (i.e., by the flexure of the elastic oceanic lithosphere) by assuming several parameters, such as the elastic thickness of the lithosphere and flexure rigidity (e.g., Liu 1980; McAdoo et al. 1978; Walcott 1970; Hanks 1971; Levitt and Sandwell 1995; Turcotte et al. 1978; Bry and White 2007; Hunter and Watts 2016). However, the estimated parameters (e.g., elastic thickness) vary significantly among these studies (e.g., Levitt and Sandwell 1995; Zhang et al. 2020). In addition, seafloor topographic highs like the Hokkaido Rise are not formed in adjacent areas northeast and southwest of the Hokkaido Rise even though plate flexure is expected to be similar (Fig. 14). These results and observations suggest that the prominent topographic high of the Hokkaido Rise cannot be explained by a simple flexure of the elastic lithosphere. Thus, the formation process of the Hokkaido Rise has remained enigmatic.

Fig. 14

Wide-area bathymetry on the northwestern part of the Pacific plate. Contour lines of 5400 and 5500 m water depth are plotted as solid lines to show the approximate extent of the Hokkaido Rise (Kobayashi et al. 1998). The dashed black lines are contour lines of the plate age (Seton et al. 2020; Gee and Kent 2007). The red line located in the forearc region of the Japan Trench is a segment boundary of the forearc structure defined by the residual topography and residual gravity anomalies (Bassett et al. 2016). The bathymetry data are from GEBCO Compilation Group (2021)

As described above, the good correlation between seafloor depth and crustal thickness implies that at least a part of the topographic high of the Hokkaido Rise can be explained by isostasy. Assuming local isostasy, the difference in the seafloor depth ($\Delta$) between the thick crust area and the thin crust area can be calculated as follows:

$$\begin{aligned} \Delta &= D_{\textrm{thin}} - D_{\textrm{thick}}\\ &= (H_{\textrm{thick}} - H_{\textrm{thin}}) \frac{(\rho _m - \rho _c) }{(\rho _m - \rho _w) } \end{aligned}$$

(1)

where $D_{\textrm{thick}}$ and $D_{\textrm{thin}}$ are the seafloor depth, and $H_{\textrm{thick}}$ and $H_{\textrm{thin}}$ are the crust thickness of each area, respectively, and $\rho _m$, $\rho _c$, and $\rho _w$ are the densities of the mantle, crust, and water, respectively. In addition, assuming the mean crustal velocity (v), we can rewrite equation 1 as:

$$\begin{aligned} \Delta \simeq \frac{v}{2} (T_{\textrm{thick}} - T_{\textrm{thin}}) \frac{(\rho _m - \rho _c) }{(\rho _m - \rho _w) } \end{aligned}$$

(2)

where $T_{\textrm{thick}}$ and $T_{\textrm{thin}}$ are the two-way traveltime in the thick and thin crust areas, respectively. Based on the $V_p$ models shown in Fig. 9, we adopted a mean crustal velocity of $v = 6.5$ km/s. For densities, we used $\rho _w = 1.03$, $\rho _c= 2.86$ and $\rho _m = 3.30$ after Hoggard et al. (2017).

Using Eq. (2) and these parameters, we compared two areas: (1) the thick crust area at the peak of the Hokkaido Rise where $T_{\textrm{thick}}$ is $\sim 2.2$ s and the seafloor depth is $\sim 5200$ m, and (2) the thin crust area at the eastern part of lines A4 and A6 where $T_{\textrm{thin}}$ is $\sim 1.5$ s and the seafloor depth is $\sim 5800$ m. Substituting these two-way times into Eq. (2), $\Delta$ was estimated to be 440 m. Although this is a rough estimation, it suggests that the isostatic equilibrium can explain much more than half of the difference in the seafloor depth of $\sim 600$ m between two areas. In both the thick and thin crust areas, the sedimentary layers, including the chert unit at the bottom of the sediments, are not disturbed (Fig. 11), suggesting that the oceanic plate has not experienced significant tectonic processes after its formation. Therefore, we consider that the difference in the oceanic crust thickness between the Hokkaido Rise and the surrounding area originates from the plate formation process at the mid-ocean ridge, meaning that the Hokkaido Rise has been relatively shallower than the surrounding area since its formation.

The Hokkaido Rise elongates to more than 1200 km along the Kuril Trench, and our study area was located at the southwestern part, where the height of the Hokkaido Rise is the lowest. The seafloor depth on the Hokkaido Rise becomes shallower toward the northeast, suggesting that the thickness of the oceanic crust becomes thicker in that direction by assuming isostatic equilibrium. Interestingly, the horizontal extent of the Hokkaido Rise is roughly comparable to that of the Shatsky Rise, one of the Large Igneous Provinces (LIPs), which was formed at the trace of the junction of three spreading ridge axes. It has a different tectonic history and a much thicker crust ($\sim 30$ km) at its peak (e.g., Sager et al. 1999, 2016; Korenaga and Sager 2012). However, our observation suggests that the thickness of the oceanic crust, which has conventionally been considered to be nearly uniform other than near LIPs and large seamounts, actually has large-scale spatial variations. These variations can be as large as the horizontal extent of the Shatsky Rise and are much larger than the seamount-scale or even the scale of a seamount chain.

5.3 The deep ocean basin near the Nosappu Fracture Zone

The thinnest oceanic crust we found was located at the eastern parts of lines A4 and A6. The crustal thickness in these areas was less than 1.5 s in two-way time, corresponding to a thickness of less than 4.9 km by assuming a mean $V_p$ of 6.5 km/s. The seafloor depth in this area reaches more than 6000 m, which is deeper than the depth predicted by the plate age-depth curve based on the thermal cooling model (e.g., Müller et al. 2008; Crosby et al. 2006).

Compilation studies of oceanic crustal structure around the world have shown that the thickness of the oceanic crust is basically proportional to both the plate age and the spreading rate at the mid-ocean ridge (Grevemeyer et al. 2018; Van Avendonk et al. 2017; Christeson et al. 2019). Considering the plate age ($\sim 135$ Ma) and spreading rate (140 mm/year) of this deep ocean basin (Seton et al. 2020), the expected thickness is $\sim 6.5$ km, much thicker than the observed crustal thickness of less than 5 km.

One of characteristics of this deep ocean basin is its weak magnetic anomalies (e.g., Nakanishi et al. 1989, 1992b; Nakanishi 1993; Choe and Dyment 2020). The geomagnetic field was recorded in the oceanic crust when the oceanic plate formed at the spreading ridge, and the oceanic crust retains this geomagnetic information until the crust is heavily altered (e.g., Vine and Matthews 1963). Because the amplitudes of the magnetic anomalies are expected to be dependent on the thickness of the oceanic crust (Vine and Matthews 1963), the thin oceanic crust in this deep ocean basin might be responsible for the weak and obscured magnetic anomaly lineations.

In addition, the topmost mantle showed strangely high $V_p$ values ($\sim 8.5$ km/s) along lines A4 and A6 in this deep ocean basin, which were considerably higher than those in the central part of these lines (7.8–8.0 km/s; Figs. 8h and 9). The differences in mantle $V_p$ along these lines were also readily visible on the record section as changes in apparent velocity (Figs. 7 and S1, Supplementary Fig. S2 of Fujie et al. (2020)).

Generally, the topmost oceanic mantle formed at a fast-spreading ridge exhibits strong $V_p$ anisotropy due to the preferred orientation of the olivine caused by the fast mantle flow at the spreading ridge, and the orientation of the high $V_p$ is perpendicular to the strike of the ridge (e.g., Hess 1964). In the northwestern part of the Pacific plate, Kodaira et al. (2014) confirmed $V_p$ anisotropy of nearly 10% in the topmost mantle near the intersection of lines A2 and P1 (Fig. 1); they found a $V_p$ of $\sim 8.5$ km/s along the spreading direction (line A2) and $\sim 7.8$ km/s along the paleo-ridge axis (line P1). Considering the angle between the survey lines and the magnetic anomaly lineation, and assuming sinusoidal patterns of seismic anisotropy (e.g., Shinohara et al. 2008), $V_p$ values near the central parts of lines A4 and A6 were consistent with the expected values, but $V_p$ values near eastern parts of these lines were notably higher.

This deep ocean basin is located on the southwestern side of the Nosappu Fracture Zone (NFZ). Based on the intensity of the magnetic anomalies around the NFZ, Nakanishi (1993) suggested that the oceanic plate formation processes at the mid-ocean ridge might have been unstable around the NFZ. The strangely high $V_p$ values along the eastern part of lines A4 and A6 imply the possibility that the $V_p$ anisotropy is oblique to the paleo-spreading axis and might be explained by the unstable and unusual plate formation processes near the NFZ. The age of the oceanic plate of this deep ocean basin is 134–140 Ma, which corresponds to the age when one of the prominent massifs of the Shatsky Rise, the Ori Massif, was formed at the nearby triple junction (Sager et al. 2016; Seton et al. 2020), although the relation between the thin crust here and LIP formation is not clear. To better understand the formation process of this deep basin and its potential relation to the nearby LIPs, we need to collect more geophysical data to constrain the actual thickness of the crust and the seismic anisotropy around the NFZ.

5.4 Subduction of the Hokkaido Rise

As shown in Fig. 14, the Hokkaido Rise is located just outside of the Kuril Trench and the northern Japan Trench. When evaluating the link between the Hokkaido Rise and the subduction zone, one key question is whether the Hokkaido Rise has already started to subduct beneath the forearc. Figure 13 shows the oceanic crust is thick near the axis of the Kuril Trench, suggesting that the Hokkaido Rise has already entered into the trench. In the northern Japan Trench, however, data coverage is relatively poor, and it is difficult to identify the crustal thickness beneath the well-developed horst and graben structure. It is therefore more difficult to answer the question based on our mapping results alone, so we focused on the bathymetric data because the Hokkaido Rise shows a relatively shallower seafloor depth than the surrounding areas. As Fig. 15 shows, the seafloor depth along the trench axis shows systematic variations: the typical maximum depth is $\sim 7,100$ m in the Kuril Trench, $\sim 7,500$ m in the northern Japan Trench, and $\sim 7,800$ m in the southern Japan Trench. These maximum seafloor depths are roughly 2, 000 m deeper than those at the top of the outer rise in each area. Considering the thickness of the trench sediments varies by at most a few hundred meters as shown in Figs. 2 and 11a, the seafloor depth variations among these areas are considered to be a good proxy of the depth variations of the top of the oceanic crust. Therefore, we conclude that the seafloor topographic high of the Hokkaido Rise has already reached at least the trench axes of the Kuril Trench and the northern Japan Trench.

Fig. 15

Seafloor bathymetry focusing on the trench axis. The maximum seafloor depths are $\sim 7100$ m in the Kuril Trench, $\sim 7500$ m in the northern Japan Trench, and $\sim 7800$ m in the southern Japan Trench

The Hokkaido Rise has a thicker crust than the surrounding areas, so it is considered to be more buoyant. Thus, subduction of the Hokkaido Rise would be expected to affect the gravity and seafloor topography of the forearc region (e.g., Vogt 1973). Bassett et al. (2016) calculated residual topography and residual gravity anomalies in the Japan Trench and found a clear boundary that divides the Japan Trench into northern and southern parts (the red line in the left middle part of Fig. 14). There were positive residual gravity anomalies to the north of the line and negative ones to the south. Generally, the residual gravity shows positive anomalies if the subducting oceanic plate has a thick and buoyant oceanic crust, such as in a seamount area. Indeed, Bassett et al. (2016) showed that positive gravity anomalies where the Joban seamount chains are considered to be subducting. Although we cannot rule out the possibility that the variations in the residual topography and residual gravity anomalies are mainly caused by the upper plate structure based on the seismic structure before subduction alone, the positive anomalies to the north of the boundary are consistent with the subduction of the Hokkaido Rise. In addition, the boundary of the residual gravity is consistent with the trenchward extension of the southern edge of the Hokkaido Rise (Fig. 14). Therefore, we propose that the Hokkaido Rise has already entered into the Japan Trench north of $37.5^\circ$N, and the forearc structural boundary of the residual topography and the residual gravity is related to the subduction of the Hokkaido Rise.

5.5 Implications for the subduction process

Interplate seismic activities show clear differences between the northern and southern Japan Trench. As Bassett et al. (2016) pointed out, all of the observed large interplate earthquakes (M$>7$) in the Japan Trench subduction zone to date, including the 2011 Tohoku earthquake, have been located on the northern side of the boundary of the residual topography and residual gravity anomalies (the red line in Fig. 14), suggesting differences in interplate seismic coupling between the northern and southern Japan Trench. Factors controlling the differences in the interplate seismic coupling are still under debate, but the presence or absence of the Hokkaido Rise may contribute to the differences in the interplate seismic coupling because subduction of buoyant oceanic crust such as that of the Hokkaido Rise is generally expected to increase normal stress at the plate interface (e.g., Kelleher and McCann 1976; Scholz and Small 1997)

There are remarkable differences between the northern and the southern Japan Trench in the seismic structure after subduction. In the northern Japan Trench, a prism-shaped low-velocity wedge (frontal prism composed of accretionary sediments) is developed at the toe of the inner trench slope (Kodaira et al. 2017; Tsuru et al. 2002). In contrast, in the southern Japan Trench, there is no prism-shaped zone; rather, a channel-like sediment unit exists along the plate interface. The presence of this type of sediment unit along the plate interface is considered to be the key factor weakening the interplate seismic coupling (Tsuru et al. 2002; Nakata et al. 2021). Factors forming these structural features are poorly understood, but the good agreement of the boundary between the northern and southern Japan Trench and the southern edge of the Hokkaido Rise (Fig. 14) implies that the subduction of the shallow and buoyant Hokkaido Rise influences the formation of these structural features. To test this hypothesis, future works, including 3-D numerical simulation of the subduction of large-scale topographic irregularities based on practical geophysical and petrological constraints, are necessary.

Morell (2016) pointed out, based on various geophysical observations in the southern Central America subduction zone, that there are strong links between the upper plate deformation and the nature of the incoming plate, such as its surface topography and the crustal thickness. The impact of subduction of the larger size of bathymetric features such as the Hokkaido Rise has a much larger and longer-lasting impact than smaller features on the scale of seamounts. In fact, the structural impacts caused by a large seamount subduction are on a scale of a few tens of kilometers (e.g., Kodaira et al. 2000; Mochizuki et al. 2008; Bangs et al. 2006; Collot et al. 2017), whereas structural differences between the northern and the southern Japan Trench are a few hundred kilometers in scale. Sediment thickness on the oceanic plate also shows significant spatial variations (Fig. 6a), but the characteristic scale is significantly smaller than that of the large-scale structural difference between the northern and the southern Japan Trench (Fujie et al. 2020). Therefore, we suggest that the subduction of the buoyant Hokkaido Rise along the Japan Trench is one plausible explanation for the large-scale structural differences between the northern and the southern Japan Trench.

To investigate whether the Hokkaido Rise has already subducted, one effective and certain approach is to investigate the thickness of the oceanic crust after subduction. Wide-angle seismic reflection and refraction data might have more potential to constrain the thickness of the subducting oceanic crust than MCS data. Some previous wide-angle seismic studies have reported a slightly thicker oceanic crust, $7.5 \sim 9$ km, in the Kuril Trench (Nakanishi et al. 2004) and in the northern Japan Trench (Fujie et al. 2006; Ito et al. 2004; Fujie et al. 2002). However, it is difficult to evaluate the credibility or uncertainty of the estimates of the thickness of the subducting oceanic crust in these earlier studies because of limitations in the amount and quality of seismic data as well as in the modeling technique used. To confirm the subduction of the Hokkaido Rise, we need to model the thickness of the subducting oceanic crust with a small enough uncertainty that we can distinguish between thicknesses of 6 and 7 km. For this purpose, it is desirable to acquire high-quality wide-angle seismic data that would enable us to apply modern imaging techniques, such as full-waveform inversion (e.g., Górszczyk et al. 2017).

6 Conclusions

To investigate the nature of the subduction inputs to the northeastern Japan arc, we have conducted large-scale controlled-source seismic surveys on the northwestern part of the oceanic Pacific plate for more than 10 years. The large amount of obtained seismic data have revealed that the incoming oceanic plate shows significant spatial variations caused by various factors. Those factors include:

• Structural evolution just prior to subduction due to plate bending, which shows reduction in $V_p$ and an increase in the $V_p/V_s$ ratio within the crust and the topmost mantle toward the trench and is accompanied by the development of the horst and graben structure on the seafloor.

• Differences in bend faulting and structural evolution between the Kuril Trench and the Japan Trench.

• Structural heterogeneity associated with a pseudofault, an ancient scar of the oceanic plate generated near the mid-ocean ridge.

• Disturbance of the sediment structure due to recent monogenetic volcanic activities related to the petit-spot and its potential impacts on the coseismic slip propagation of the 2011 Tohoku earthquake.

• Thickness variations of the oceanic crust, suggesting that the Hokkaido Rise, a remarkable topographic high located in the outer trench area of the Kuril Trench and the northern Japan Trench, is explained by isostatic equilibrium to a certain extent.

• Subduction of the buoyant thick-crust Hokkaido Rise, which might be related to differences in seismic velocity structure and activities, including megathrust earthquakes between the northern and southern Japan Trench.

In conclusion, the oceanic Pacific plate as a subduction input to the northeastern Japan arc shows a great degree of structural variations on various scales. The structural variations of the subduction inputs show good pattern agreement with the seismic structure and the distribution of the seismic activities after subduction. To reveal the link between the subduction inputs and processes in the subduction zone, it is desirable to conduct comprehensive future works, including structural studies across the trench based on various geophysical data, numerical simulations, petrological observations, and experiments.