The 2016 Mw 5.9 earthquake off the southeastern coast of Mie Prefecture as an indicator of preparatory processes of the next Nankai Trough megathrust earthquake

The 2016 Mw 5.9 earthquake off the southeastern coast of Mie Prefecture as an indicator of... Megathrust earthquakes have occurred repeatedly at intervals of 100 to 150 years along the Nankai Trough, situated in the southwest of Japan. Given that it has been 70 years since the last event, the occurrence of the next devastating earthquake is anticipated in the near future. On April 1, 2016, a moderate earthquake (M 5.9, M 6.5) occurred off the southeastern coast of Mie Prefecture in the source region of the 1944 Tonankai JMA earthquake (M 8.2). In this study, we investigated the influence of the 2016 earthquake on future megathrust earthquakes. We first determined the hypocenter distributions using a precise velocity structure obtained from seismic surveys in the source region. Using data obtained from the DONET ocean-bottom observation network, we found that this earthquake occurred along the plate boundary fault, which is also believed to have slipped during past megathrust earthquakes. We then performed a preliminary numerical simulation to reproduce the occurrence of a moderate earthquake in the middle of a megathrust earthquake cycle. We used a hierarchical asperity model, in which smaller asperities causing moderate earthquakes are embedded in a hyper-asperity that serves as the source region of megathrust earthquakes. The simulation shows that moderate earthquakes, caused by ruptures of a smaller asperity, occur as a result of shrinkage of strongly coupled areas in the hyper-asperity. This result is consistent with the observation that the hypocenter of the 2016 earthquake was located at the edge of a strongly coupled region along the plate boundary. The simulation also reproduced postseismic slip (afterslip/slow slip) along the plate boundary updip of the hyper-asperity, which is consistent with the observations of slow slip events and very-low-frequency earthquakes after the 2016 earthquake. Thus, the occurrence of the moderate earthquake offshore southeastern Mie Prefecture in the middle of the megathrust earthquake cycle implies that the shrinkage of the strongly coupled area along the plate boundary is occurring as a preparatory process for the next megathrust earthquake in the region. Keywords: Hypocenter determination, Crustal structure, Ocean floor observations, Accretionary prism, Earthquake cycle simulation, Hierarchical asperity model * Correspondence: mnakano@jamstec.go.jp R&D Center for Earthquake and Tsunami, Japan Agency for Marine-Earth Science and Technology, Yokohama, Japan Full list of author information is available at the end of the article © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 2 of 17 Introduction seismic activities (Kaneda et al. 2015; Kawaguchi et al. At convergent plate boundaries, large earthquakes occur 2015). These observations revealed that, in recent years, repeatedly with intervals of one to several hundred years. earthquake activities along the Nankai Trough have The Nankai Trough, where the Philippine Sea plate mostly occurred in the oceanic crust or upper mantle of (PSP) subducts beneath the Amur plate (Bird 2003; the subducting PSP (Obana et al. 2004, 2005; Obana and DeMets et al. 2010), is one such region where mega- Kodaira 2009; Mochizuki et al. 2010; Akuhara et al. thrust earthquakes and tsunamis have occurred at inter- 2013; Nakano et al. 2013, 2014, 2015; Akuhara and vals of 100 to 150 years (e.g., Ando 1975; Ishibashi Mochizuki 2014). The most recent notable events during 2004). The 1944 Tonankai (M 8.2; Ichinose et al. 2003) this time have been the 2004 M 7.1, 7.4, and 6.5 w JMA and 1946 Nankai (M 8.4; Murotani et al. 2015) earth- earthquakes off the Kii Peninsula, which occurred in the quakes were two such devastating events, which rup- subducting plate (e.g., Sakai et al. 2005; Seno 2005). tured the eastern and western segments of the trough, Significant improvements in supercomputers have fa- respectively (Fig. 1), causing severe damage in western cilitated intensive numerical simulations of megathrust and central Japan. Future earthquakes with matching in- earthquakes (e.g., Hori 2006; Kodaira et al. 2006; Hori et tensities will almost certainly result in severe damage, al. 2009; Hyodo and Hori 2013). These studies have given the strong ground motion and tsunamis likely to mainly focused on reproducing the slip distributions, be associated with them. Recent studies have indicated rupture segmentations, and recurrence intervals of his- that such megathrust earthquakes would also cause torical earthquakes, based on realistic fault models and massive damage due to long-period ground motion in subduction zone structures obtained from seismic sur- distant regions, including metropolitan Tokyo (e.g., Fur- veys. In contrast, simulations based on simplified fault umura et al. 2008). models can provide us with insights into earthquake In order to clarify the deformation process along the rupture patterns, cycles, and related phenomena such as convergent margin and the potential for such disasters, slow slip events (SSEs) (e.g., Matsuzawa et al. 2010; interdisciplinary studies, including seismic surveys of Noda et al. 2013, 2014; Ariyoshi et al. 2014). Hori and subsurface structures, seismic monitoring, geodetic ob- Miyazaki (2011) simulated occurrences of M7–8 earth- servations, and numerical simulations, have been con- quakes during M9 earthquake cycles to reproduce the ducted. A number of seismic surveys have investigated occurrence of the 1978 M 7.5 earthquake that occurred geophysical and geological structures in detail (e.g., offshore Miyagi before the 2011 Tohoku-oki earthquake Mochizuki et al. 1998; Nakanishi et al. 1998, 2008; Park (M 9.0). They used the hierarchical asperity model, in and Kodaira 2012; Park et al. 2002, 2003, 2010; Kamei et which a smaller asperity causing M7−8 earthquakes was al. 2012). One of the most remarkable structures in the embedded in a hyper-asperity that caused the M9 earth- accretionary prism along the Nankai Trough is the quake; namely, the fundamental rupture mode in the megasplay fault or the out-of-sequence thrust branching subduction zone. Numerical simulations can also indi- from the decollement, imaged as strong reflections in cate the temporal evolution of coupling ratios and stress multi-channel seismic survey (MCS) profiles (Park et al. states along fault interfaces. These parameters are, how- 2002; Moore et al. 2007). Recent studies, based on de- ever, still very difficult to estimate from observations but tailed seismic surveys, have interpreted the megasplay are considered to be crucial for monitoring the prepara- fault as the main plate boundary fault (Bangs et al. 2009; tory processes of megathrust earthquakes (e.g., Kodaira Kamei et al. 2012, 2013; Tobin et al. 2014; Tsuji et al. et al. 2006). Thus, combinations of numerical simula- 2014). The megasplay fault branches in the shallower tions with seismic and geodetic observations could pro- part; one section reaches the ocean bottom and the vide information about the current state of stress and other connects to the decollement beneath the toe of strain along plate boundary faults. the accretionary prism (Tsuji et al. 2014). Both branches On April 1, 2016, an earthquake occurred off the are considered to have slipped during past megathrust southeastern coast of Mie Prefecture (M = 5.9 as per earthquakes (Sakaguchi et al. 2011; Tsuji et al. 2014). the U.S. Geological Survey) (“off-Mie earthquake” here- Extensive observations of seismic activities along the after). It occurred updip of the hypocenter of the 1944 Nankai Trough have been conducted using temporary Tonankai earthquake (Fig. 1a). The hypocenter was be- deployments of recoverable ocean-bottom seismometers neath the DONET stations, and real-time ocean-bottom (OBS) (e.g., Obana et al. 2004; Obana and Kodaira 2009; observations were obtained. The epicenter was within Mochizuki et al. 2010). The recent development of a the slip area of the 1944 Tonankai earthquake (Ichinose permanent seafloor observation network called the et al. 2003) and close to the eastern margin of the 1946 Dense Oceanfloor Network System for Earthquakes and Nankai earthquake slip area (Murotani et al. 2015). Dur- Tsunamis (DONET) (Fig. 1 and Additional file 1: Figure ing the 70-year period since the occurrence of the last S1) has facilitated continuous monitoring of offshore megathrust earthquake, earthquakes with magnitudes Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 3 of 17 Fig. 1 (See legend on next page.) Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 4 of 17 (See figure on previous page.) Fig. 1 Past seismic activities in the study area. a Regional location map (inset) and seismic activity along the Nankai Trough. Orange star represents the epicenter location of the 2016 off-Mie earthquake; results from the Japan Meteorological Agency (JMA), U.S. Geological Survey (USGS), and the Global CMT Project (GCMT) are shown together with mechanism symbols. Epicenters of earthquakes larger than M6 before and after 1976 are represented by diamonds and circles, respectively, from the JMA unified earthquake catalogue. Blue and red contours represent the slip distributions of the 1944 Tonankai and 1946 Nankai earthquakes, respectively (Ichinose et al. 2003;Murotani et al. 2015). Gray triangles represent locations of DONET stations. The ellipsoid labeled with “A” designates the region of interplate seismic activity reported by Akuhara and Mochizuki (2014). Orange arrows represent the subduction direction of the Philippine Sea plate. Plate name abbreviations are as follows: AM, Amur; OK, Okhotsk; PA, Pacific; PS, Philippine Sea. b Occurrence date of earthquakes plotted against their source longitude. Symbols are the same as in a larger than 6 have been very rare along the eastern Nan- that the sources were distributed very close to the plate kai Trough (Fig. 1b). Exceptions were aftershocks of the interface. However, strong horizontal heterogeneities of 1944 and 1946 earthquakes that continued for about seismic wave velocities due to plate subduction could 10–15 years and the sequence of the 2004 intraplate affect the precision of determining the hypocenter loca- earthquakes off the Kii Peninsula. Hence, large interplate tions. Takemura et al. (2016) showed that an accurate earthquakes have not been known to occur until the 3D velocity model reflecting the plate subduction is middle of a megathrust earthquake cycle. Wallace et al. necessary to explain the spatial distributions of P-wave (2016) consider that the 2016 off-Mie earthquake most first motion polarities observed for this earthquake. likely occurred on the plate interface, based on inte- Therefore, a velocity model reflecting the horizontal grated analyses of hypocenter distributions, seafloor heterogeneities is crucial for precise determination of crustal deformations, and tsunami modeling using sea- hypocenters. floor geodetic and seismological data from DONET. In this study, we determined the hypocenters of the After this earthquake, SSEs along the plate boundary, 2016 off-Mie earthquake and its aftershocks using a vel- very-low-frequency earthquakes (VLFEs), and tectonic ocity model that includes horizontal heterogeneities tremors have been observed updip of the hypocentral re- along the subduction zone. We used a velocity structure gion (Fig. 2) (Annoura et al. 2017; Araki et al. 2017; obtained from a wide-angle OBS survey along a line Kaneko et al. 2018; Nakano et al. 2018). Recent studies passing through the source region shown by Wallace et have shown that signals from these slow events radiated al. (2016). We compared the hypocenter distribution from a common source (Kaneko et al. 2018; Nakano et with the MCS reflection profiles to confirm whether this al. 2018). The development of a comprehensive model of earthquake occurred on the plate boundary. Since this earthquake sources that can explain the generation of earthquake was located along the plate boundary fault, these regular and slow earthquakes would help us to we conducted a preliminary numerical simulation of understand the processes of stress accumulation and/or plate boundary earthquakes to understand the implica- release along the plate boundary. It is important to de- tions of this earthquake for the Nankai Trough mega- termine the precise source location of the 2016 off-Mie thrust earthquake cycle. Using the hierarchical asperity earthquake to evaluate its influence on the generation of model, we found that the 2016 off-Mie earthquake oc- the next megathrust earthquake and its preparatory curred due to shrinkage of the strongly coupled area processes. along the plate interface, which is likely to be a prepara- The centroid moment tensor (CMT) solution for the tory process for the next megathrust earthquake. Inter- mainshock shows a thrust-type mechanism consistent estingly, this model can also explain the absence of large with plate boundary slip (Fig. 1a). However, the source earthquakes until the middle of the megathrust earth- centroid location, which is based on regional or teleseis- quake cycle. mic data, is uncertain; depending on the institution that determined the location (Japan Meteorological Agency, Results and discussion U.S. Geological Survey, and the GlobalCMT project), Hypocenter distribution of the 2016 off-Mie earthquake there is a variation of ~ 20 km in the horizontal direc- and aftershocks tion and 7 km (between 12 and 19 km) in depth. Such a We determined the hypocenters of the mainshock and discrepancy may have occurred because of strong struc- aftershocks of the 2016 off-Mie earthquake using a local tural heterogeneities that exist in the dip direction of velocity model obtained in the source region (Fig. 2). plate subduction. Precise determination of the source lo- We used a horizontally heterogeneous velocity structure cation, especially for the depth, is crucial to clarify based on that obtained along a survey line that passes whether this earthquake occurred along the plate bound- through the hypocentral region (see the section ary. Wallace et al. (2016) determined the hypocenter dis- “Methods/experimental”). The hypocenters were distrib- tribution using a 1D velocity structure, which showed uted several kilometers shallower than those based on Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 5 of 17 Fig. 2 (See legend on next page.) Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 6 of 17 (See figure on previous page.) Fig. 2 Hypocenters of the off-Mie earthquake and aftershocks. a Hypocenter distribution of the 2016 off-Mie earthquake and aftershocks (red circles). Pink focal mechanism symbols indicate shallow VLFEs, and magenta dots are locations of tectonic tremors in April 2016 (Kaneko et al. 2018; Nakano et al. 2018). White triangles represent DONET stations used for our hypocenter relocation. Green and blue lines represent location of MCS survey lines shown in Fig. 5. The thin black line represents the location of MCS survey line (COP2), and the thick black segment represents the location of a topographic high (Tsuru et al. 2005). The purple dashed line indicates the boundary of the topographic high in the basin. Red and blue dashed contours represent regions of low- and high-velocity anomalies, respectively, in the overriding plate (Yamamoto et al. 2017). A high-velocity body off the Kii Peninsula is represented by an orange oval (Kodaira et al. 2006). b Velocity structure used for our hypocenter determination obtained along MCS survey line KI03 (Wallace et al. 2016) our routine hypocenter determinations using a 1D vel- but the depth difference from the plate boundary fault ocity structure (Fig. 3). The mainshock depth was 9.2 ± was smaller in our result. The hypocenter location corre- 0.6 km below sea level, about 4 km shallower than that sponds to the shallower edge of the strongly coupled re- obtained by the routine analysis, and about 2 km shal- gion along the plate boundary, as noted by Yokota et al. lower than that obtained by Wallace et al. (2016). The (2016) (Additional file 2: Figure S2). differences can be attributed to the horizontal heteroge- neities of the velocity structure considered in our ana- Earthquake cycle simulation based on the hierarchical lysis. Aftershocks were distributed in a very limited area, asperity model about 10 km downdip of the mainshock, at depths of be- In order to investigate the generating mechanism of the tween 10 and 14 km, with source location errors of off-Mie earthquake in the middle of the megathrust 1.9 km in both horizontal and depth directions. There is earthquake cycle along the Nankai Trough, we con- a seismic gap distinct from the mainshock hypocenter, ducted a preliminary earthquake cycle simulation based which agrees with the results obtained by Wallace et al. on the hierarchical asperity model (Hori and Miyazaki (2016). The aftershock activity continued for 2 weeks, 2010, 2011) that has been used previously to simulate and a burst of reactivation occurred on April 19 (Fig. 4), multiscale earthquakes. In this model, smaller asperities, coinciding with the end of the SSE following the main- causing moderate-sized earthquakes, are embedded in a shock (Araki et al. 2017). larger asperity (“hyper-asperity”), causing the fundamen- The aftershock distribution can be further separated tal rupture mode of the fault system. To simulate the into two clusters: northeastern and southwestern generation of earthquakes on a megathrust, we solved (Fig. 3a). Events in the northeastern cluster were distrib- the quasi-dynamic equation of shear stress along fault uted almost parallel to the strike of the PSP, with a segments (Rice 1993) using the rate- and state- near-identical source depth, while those in the south- dependent friction (RSF) laws to represent frictional western cluster showed a rather scattered distribution in properties across the fault (see the section titled the dip direction. Event sizes were systematically larger “Methods/experimental”). The occurrence of an earth- in the northeastern cluster, and decay times of these ac- quake on the plate boundary is modeled as a perturb- tivities were clearly different. Activity in the northeastern ation from the steady-state relative plate motion. Shear cluster lasted for only a day, whereas, in the southwest- stress builds up due to secular motion of the subducting ern cluster, it continued for much longer (Fig. 4). plate which is finally released as earthquakes when the Figure 5 compares the earthquake depths and reflec- accumulated shear stress exceeds the fault strength gov- tors imaged from MCS surveys along lines passing erned by RSF. through the source region (see the section “Methods/ex- In our simulation, we attempted to reproduce the fol- perimental”). Clear reflections indicative of the plate lowing characteristics observed for the off-Mie earth- boundary are observed in the MCS profiles. The main- quake. First, the earthquake occurred after 70 years of a shock was located at a slightly shallower depth than the quiet period characterized by the absence of M6 or lar- upper reflector, possibly corresponding to the megasplay ger earthquakes after the previous M8 earthquake. Sec- fault. The aftershocks were distributed about 10 km ond, a slow slip and/or afterslip episode was observed landward from the mainshock, at depths about 2 km immediately after this earthquake (Wallace et al. 2016; deeper than the reflector (Fig. 3). This result agrees with Araki et al. 2017). In order to reproduce such compli- that of Wallace et al. (2016), although the separation cated slip behavior during an earthquake cycle, we need depth was comparable to the location errors. to properly define the frictional parameters along the Considering the location errors associated with deter- megathrust, as well as the sizes and locations of asper- mining the hypocenters and reflector depths (±1.1 km), ities. Suitable parameters need to be determined empir- we consider that this earthquake has slipped along the ically and the computational cost for conventional 3D plate boundary fault, as Wallace et al. (2016) reported, simulations is high; hence, we conducted 2D simulations Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 7 of 17 Fig. 3 (See legend on next page.) Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 8 of 17 (See figure on previous page.) Fig. 3 Comparison of routine and relocated hypocenter locations. Yellow and red circles represent routine and relocated hypocenters, respectively. Black and red crosses at the top left of each panel represent the uncertainties of hypocenter locations for the mainshock and aftershocks, respectively. a Map view. Gray triangles represent locations of DONET stations. b Cross-section along line E–F. Black and blue lines represent strong reflections picked from the MCS profile shown in Fig. 5a (temporal evolution of slip on a fault, represented by line updip of the hyper-asperity, and continues for more than segments) for the parameter search. In this study, we 20 days (Fig. 7b). After the occurrence of the M6 event, present one of the simulation results that reproduced the strongly coupled area in the hyper-asperity further the characteristics described above well. shrinks because of the secular aseismic slip in the sur- In the numerical simulation, the frictional parameters rounding area, and the shear stress of the deeper asper- of asperities are set as unstable slip conditions, while re- ity further increases (Fig. 8). When the stress at the gions outside of them are set as stable. Since the hypo- deeper asperity reaches the threshold, a rupture initiated center of the off-Mie earthquake was located at the at the asperity propagates along the entire hyper-asperity shallower edge of the strongly coupled region along the and an M8 earthquake occurs. Nankai Trough (Additional file 2: Figure S2), we defined Using earthquake cycle simulation based on the hier- a smaller asperity, causing an M6-class earthquake, at archical asperity model, we successfully reproduced the the shallower edge of the hyper-asperity. Figure 6 shows occurrence of a moderate earthquake in the middle of the fault geometry and frictional parameters, in which the megathrust earthquake cycle and the SSE (or after- we assume two smaller asperities (dark blue segments) slip) that follows the earthquake. The rupture of the embedded in the hyper-asperity (both dark blue and shallower asperity is caused by shrinkage of the strongly light-blue segments). Based on a simulation using this coupled area along the hyper-asperity, which results frictional parameter distribution, the rupture of a shal- from the secular plate subduction. This finding implies lower asperity placed on a 10-km-long line segment that the occurrence of the moderate 2016 off-Mie earth- along the plate boundary fault, corresponding to a 10– quake in the middle of the earthquake cycle represents 11 km depth range (“shallower asperity” hereafter), shrinkage of a strongly coupled area, which could be a caused the M6 earthquake. Whereas, the rupture of a preparatory process of the next megathrust earthquake deeper asperity, a 30-km-long segment along the plate along the Nankai Trough. We note here that the cycle of boundary fault corresponding to a 13–17 km depth the M8 earthquake and the timing and frequency of the range (“deeper asperity” hereafter), propagates along the M6 earthquake depend on the frictional parameters of entire hyper-asperity to cause an M8 earthquake, the the fault, their distributions, and the slip velocity at the fundamental rupture mode of this fault system. plate interface, all of which were assumed in the Figure 7a shows the spatiotemporal evolution of the simulation. slip rate along the fault. At the beginning of the M8 Aftershocks of the off-Mie earthquake occurred only megathrust earthquake cycle (marked by “M8” at the downdip of the mainshock, with a clear separation of bottom of Fig. 7a), the slip velocity on the fault is very about 10 km (similar to Wallace et al. 2016). The ab- small along the entire asperity, that is, the fault is sence of an updip aftershock is not surprising, because strongly coupled. Outside the hyper-asperity, a secular this region is considered to have conditionally stable aseismic slip occurs during the entire earthquake cycle frictional properties (Wallace et al. 2016). Instead, slow because of the stable sliding frictional properties. The earthquakes as SSEs, VLFEs, and tectonic tremors oc- strongly coupled area on the hyper-asperity gradually curred there, sharing a common updip source region shrinks from the deeper and shallower edges due to the (Araki et al. 2017; Kaneko et al. 2018; Nakano et al. secular aseismic slip, which causes a gradual increase in 2018) (Fig. 2). Both regular and slow earthquakes were shear stress on the shallower and deeper asperities almost entirely absent between the downdip aftershocks (Fig. 8). From the first half to the middle stage of the and updip slow earthquakes, although the occurrence of M8 megathrust earthquake cycle, the rate of increase of an SSE or an afterslip close to the hypocenter region stress is higher in the shallower asperity because of its was implied after the mainshock (Wallace et al. 2016; location closer to the edge of the hyper-asperity. In our Nakano et al. 2018). Our numerical simulation success- simulation, it takes about 120 years for sufficient stress fully reproduced the slow slip updip of the asperity, to accumulate to rupture the shallower asperity (M6 where stable sliding frictional parameters were defined. earthquake) in a 200-year megathrust earthquake cycle. Although the reason for the existence of a gap between The M6 earthquake is followed by postseismic slip the mainshock hypocenter and aftershocks is still elu- (afterslip/slow slip) that occurs mainly in a region of sive, Tsuji et al. (2017) considered that aftershocks oc- velocity-strengthening frictional parameters, that is, an curred close to the location where the plate boundary Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 9 of 17 0 0 7 14 21 Date of April 2016 Fig. 4 Cumulative number and magnitude of mainshock and aftershocks plotted against time. Circles with vertical bars represent the magnitude of the earthquake. The black line represents the cumulative number of aftershocks. Blue and red lines correspond to those occurring in the northeastern and southwestern clusters separated by a dashed line, shown in Fig. 3a decollement soles out into the top of oceanic crust layer in geological structures would cause differences in 2 (TOC). Such heterogeneous structures may cause frictional properties along aftershock faults. The after- stress accumulations, which could explain the very lim- shock region corresponds to the boundary between older ited aftershock distribution. Wallace et al. (2016) sug- (~ 14 Ma) and younger (~ 6 Ma) prisms (Tsuji et al. gested that this gap coincided with the mainshock fault 2017), which would lead to a difference in the degree of that did not trigger aftershocks because stress was heterogeneities in plate boundary structure. already released, or stress was released due to aseismic Previous studies showed that structural heterogeneities slow slip, or both. The migration of VLFE sources updip, exist in overriding and subducting plates around the away from the mainshock hypocenter, supports the latter source region of the off-Mie earthquake (Fig. 2). One of model (Nakano et al. 2018). the characteristic structures off the Kii Peninsula is a The aftershock activity was separated into two clusters, plutonic rock body with high seismic velocity and high which showed different decay times, namely, the p value density, which is coincident with the segment boundary of the Omori–Utsu aftershock decay law (Figs. 3 and 4). of megathrust earthquakes along the Nankai Trough −p The Omori–Utsu law is given as n(t)= K(t + c) , where (Kodaira et al. 2006). A topographic high in the TOC n(t) is the number of events per unit time at t since the can be recognized to the west of the aftershock area in mainshock; K, c, and p are constants (Utsu 1961). Al- the MCS reflection survey obtained along the strike dir- though source properties corresponding to the p value ection of the subducting plate (line COP2 in Tsuru et al. are not well studied, Kisslinger and Jones (1991) ob- (2005)). A boundary of the ocean floor topography (the served that high heat flow causes a high p value. They purple dashed line in Fig. 2a) is also coincident with the considered that the p value represents the stress relax- epicenter distribution of the aftershocks. This region ation time, and the aftershock decay rate is controlled by corresponds to the boundary of low- and high- velocity the temperature at the source depth. In our case, how- anomalies in the overriding plate (Yamamoto et al. ever, the distance of the aftershock clusters was only sev- 2017). The source of the off-Mie earthquake is close to eral kilometers and the temperature would not differ the boundary of the backstop and sediments (Tsuji et al. significantly. Mikumo and Miyatake (1979) conducted 2015) and also corresponds to the location where the de- numerical simulations and showed that the p value be- collement soles out into the TOC (Tsuji et al. 2017). comes smaller as the distribution of fault strength be- Tsuji et al. (2017) also pointed out that the aftershock comes more heterogeneous. Small-scale heterogeneities region corresponds to the boundary between the older Magnitude Cumulative Number Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 10 of 17 a b c d Fig. 5 Hypocenters compared with plate interfaces on sections from a MCS survey. Hypocenter distribution projected on the time-migrated MCS reflection profiles along lines A–B(a) and C–D(b) shown in Fig. 2a. Arrows represent strong reflections indicative of a plate boundary (PB) and younger prisms. These heterogeneous structures subducting or overriding plates, and there is very low ac- concentrate stress and limit the area of the aftershock tivity in the rupture area of the megathrust earthquakes activity. (Obana et al. 2004; Obana and Kodaira 2009; Nakano et In our earthquake cycle simulations, we successfully al. 2015). However, Akuhara and Mochizuki (2014) ob- reproduced the occurrence of a moderate earthquake in served continuous activity of an earthquake cluster with the middle of the megathrust earthquake cycle. The magnitudes smaller than 4 along the plate boundary simulation implies that the M6-class off-Mie earthquake west of the Kii Peninsula (labeled “A” in Fig. 1a). Their was caused by shrinkage of the strongly coupled area hypocenters roughly correspond to the downdip end of along the plate boundary. Present seismic activities in the rupture area of the 1946 Nankai earthquake (Muro- the eastern Nankai Trough are mostly within the tani et al. 2015). These observations imply that the Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 11 of 17 ab Fig. 6 Model setting of 2D earthquake cycle simulation. a Configuration of model plate interface and the schematic illustration of stable or unstable regions. Colors on the model plate interface reflect the assumed distribution of frictional parameters, (a – b)σ : red and blue colors denote stable and eff unstable regions, respectively. The shallow gray segment denotes the steady sliding zone. b Depth profiles of frictional parameters used to simulate an earthquake cycle scenario asperity of the earthquake cluster is continuously loaded distributions on a 3D model would improve our under- because of weak interplate coupling, probably due to the standing of earthquake cycles along plate boundaries. shrinkage of the strongly coupled region. In fact, the slip deficit rate in this region is less than that in the source Conclusions region of the off-Mie earthquake (Yokota et al. 2016). Using a local velocity model, we found that the 2016 Monitoring such seismic activities based on precise hy- Mw5.9 earthquake off the southeastern coast of Mie Pre- pocenters is important for the understanding of the state fecture occurred on the plate boundary along the Nankai of interplate coupling along the plate boundary fault, as Trough. We performed a numerical simulation based on well as geodetic data inversion. To achieve this, a local a hierarchical asperity model, which successfully repro- 3D velocity structure is crucial for accurate hypocenter duced the occurrence of a moderate earthquake in the locations (Akuhara et al. 2013; Yamamoto et al. 2017). middle of a megathrust earthquake cycle. The results of If other asperities exist, their ruptures could occur in the numerical simulation imply that the occurrence of the middle of the megathrust earthquake cycle due to the off-Mie earthquake represents shrinkage of the inter- the decrease in the interplate coupling ratio (Hori and seismically coupled area along the plate boundary and Miyazaki 2011). The timing and magnitude of the earth- serves as an indication of the preparatory processes of quake depends on both the location and size of the as- the next megathrust earthquake. perities. Our simulations can also model the magnitude, timing, frequency of the moderate earthquake, and the Methods/experimental interval of the fundamental rupture mode causing the Hypocenter determinations using a 2.5-D velocity model megathrust earthquake by defining the fault frictional The strong heterogeneity in the seismic velocity struc- parameters appropriately. In our 2D simulations, the ture in the subduction zone prevents precise hypocenter moderate M6 earthquake occurs 120 years after the M8 determinations. We used the 2-D velocity structure ob- earthquake, with a cycle of about 200 years. Both values tained from a tomographic inversion of the first-arrival are much longer than the observed values: the cycle of refraction travel times of a wide-angle OBS survey along megathrust earthquakes in the Nankai Trough is consid- line KI03, which passes through the hypocenter region ered to be 100 to 150 years (Ando 1975; Ishibashi 2004), (Fig. 2) (Wallace et al. 2016). Using this 2D structure, we and the off-Mie earthquake occurred 70 years after the created a 3-D velocity volume, assuming that the struc- previous 1944 Tonankai earthquake. In 3D simulations, ture is identical along the strike of the subducting plate; however, these periods would become shorter with the then, we computed theoretical travel times at each same frictional parameter values because the shrinkage DONET station. Since the PSP strikes about 40° N of the coupled area on the asperity also occurs in the around the hypocenter region (Citak et al. 2012) strike direction of the fault. Studies based on more real- (Additional file 1: Figure S1), the 2D structure is ori- istic fault geometry and frictional parameter ented in this direction (Additional file 3: Figure S3). We Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 12 of 17 ab Fig. 7 Simulation result. a Spatiotemporal evolution of the slip rate V during the resultant earthquake cycles. Colors denote the slip velocity. Black region corresponds to strong coupling on the plate interface. Arrows denote occurrences of earthquakes with denoted magnitude. The vertical color bar in the middle panel denotes the depth distribution of (a – b)σ . b Depth distributions of postseismic slip and slip rate just after the eff occurrence of a M6 earthquake indicated by solid arrows in a. Cumulative slip and slip rate profiles are shown in the left and right side panels, respectively. Colors of lines indicate days after the M6 earthquake. The broken gray line represents the coseismic slip distribution of the M6 earthquake call this the 2.5-D model for convenience because of the time at each station by averaging the residuals between similarity to the computations of wave propagation in the theoretical and observed travel times and redeter- three dimensions using a 2-D velocity model (e.g., Take- mined the hypocenter locations. This procedure was re- naka and Kennet 1996). We did not use the model of peated three times to obtain the final solution. The Kamei et al. (2012) because it did not fully cover the resulting site corrections are within ± 0.09 s, except at study area. KME19, which needed a correction of − 0.17 s. The source For the hypocenter determinations, we used the seafloor location errors are 0.7 km in the horizontal direction and seismic recordings from a broadband seismometer in- 0.6 km in the vertical direction for the mainshock and stalled at each DONET station. We manually picked 1.9 km in both directions for the aftershocks. The P-wave onsets from the vertical component seismograms. root-mean-square (RMS) of travel time residual is 0.02 s. To avoid the effect of possible heterogeneities in the ac- We examined the accuracy of the hypocenter cretionary prism in the strike direction of the PSP, we only depths by changing the dataset or assumed strike of used data from nine stations located within 30 km of the the PSP (see the Additional file 4: Supplementary mainshock source (white triangles in Fig. 2a). Moreover, note). Aftershock depths were several kilometers we did not use S-wave arrival times because of the uncer- shallower when we included distant stations tainty in S-wave velocity, especially in the soft sediments (Additional file 5:FigureS4).Thisdifference couldbe covering the accretionary prism. We used NonLinLoc caused by horizontal heterogeneities of the structure software (Lomax et al. 2000) for the hypocenter determi- along strike (e.g., Park et al. 2010). The increasing nations. This software adopts the eikonal finite-difference RMS travel time residual, including distant stations, scheme of Podvin and Lecomte (1991) to compute implies that theoretical travel times are less accurate first-arrival times at each station in a 3D volume for hypo- at distant stations because of inaccurate velocity center determinations using heterogeneous velocity struc- structures. Stable aftershock depths were obtained tures. We estimated site corrections for the P-wave travel when using S-wave arrival times or changing the Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 13 of 17 along the same line as the velocity structure, the profile in Fig. 5b was obtained along an intersecting line. Hence, the TWTs on this profile do not exactly match those computed from the hypocenter depths due to possible along strike heterogeneity of the subduction zone. We picked TWTs of strong reflections indicative of the plate boundary from the MCS reflection profile (Fig. 5a), converted them to depth values by inverting the process described above, and plotted them on the depth section (Fig. 3b). The uncertainties in picking the reflections are about 0.15 s, corresponding to 1.1 km at reflector depths. The errors in hypocenter depths correspond to a TWT error of 0.31 s for the mainshock and 0.67 s for the aftershocks. Two-dimensional forward simulation of M8-class earthquake cycles along the Nankai Trough Fig. 8 Temporal changes in shear stress on the plate interface. The Here, we attempt to interpret the conceptual relation same time interval as Fig. 7a at the four depths are shown. The blue between M8-class interplate earthquakes (such as the and red lines represent shear stress profiles at the initiation points of 1944 Tonankai earthquake) and the occurrence of a M6 (depth of 11 km) and M8 earthquakes (depth of 17 km), respectively. smaller earthquake (such as the 2016 off-Mie At depths deeper than the hypocenter of the M8 earthquake (black line), the stress accumulation rate is higher than that at the hypocenter earthquake) inside the source region through seismic because of secular loading from stable sliding occurring at a cycle modeling, which assumes the equilibrium of region deeper than the hyper-asperity. At depths shallower than quasi-dynamic shear stresses on the fault and the hyper-asperity (green lines), shear stress is almost constant frictional resistance governed by a rate- and because of stable sliding, except during the postseismic periods state-dependent friction (RSF) law. of M6 and M8 earthquakes. Dashed lines represent the reference stress level of each line, with corresponding colors Governing equations for quasi-dynamic seismic cycle strike of the PSP. Moreover, the mainshock depth was simulation stably located for most cases. Consequently, we con- In quasi-dynamic seismic cycle modeling, the occurrence of sider the obtained hypocenter locations to be stable. interplate earthquakes on plate boundaries has often been We also computed hypocenters by our routine analysis modeled as a perturbation from the steady-state relative based on a 1D velocity structure, using P- and S-wave plate motion by ignoring the contribution of steady (for- readings from all DONET stations (Fig. 3). The main- ward) slip on the whole plate interface. Thus, the interseis- shock hypocenter was located several kilometers east of mic process is expressed as the accumulation of slip deficit the 3D result, at a depth of 13.1 km. Aftershocks were (or backslip) in the seismogenic zone. The shear stress distributed several kilometers west of the 3D result, at a builds up with time as the slip deficit increases, and the ac- depth of about 15 km. We note that this computation is cumulated slip deficit is finally released as earthquakes different from that of Wallace et al. (2016), which only when the accumulated shear stress exceeds the fault resist- uses P-wave readings from stations connected to nodes ance governed by RSF. Below, we briefly describe equations D and E. governing the earthquake cycle simulation used here. At first, by discretizing the target plate interface Comparison of hypocenter depths with reflectors into N small sub-faults, we define physical variables obtained from MCS profiles associated with seismic cycles for each sub-fault: We compared the hypocenter depths with distinct re- shear stress τ,state variable θ representing contact i i flectors in the accretionary prism obtained from MCS state, and slip velocity V . Here, subscript i specifies reflection surveys (Fig. 5). We converted the hypocenter the particular sub-fault (1 ≤ i ≤ N). When the material depths to two-way travel times (TWTs) by integrating containing the fault (or a set of sub-faults) is a the slowness from sea level to the hypocenter depth linear-elastic one, temporal changes in shear stress picked at the epicenter, based upon the velocity struc- can be approximated through the whole seismic cycle ture used for the hypocenter determinations (Fig. 2b). using the quasi-dynamic expression given by Rice While the MCS profile shown in Fig. 5a was obtained (1993) as follows, Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 14 of 17 N rate of friction for a particular slip amount. Because dτ ðÞ t G ¼ − K V −V ðÞ t − V ðÞ t ; these parameters characterize the frictional property of ij pl j i dt 2V j¼1 the fault, and accordingly control the spatial and tem- poral behavior of fault slip, we must pay special atten- ð1Þ tion to the distribution of the frictional parameters on where K represents the static slip response function be- ij the model fault in seismic cycle simulations. tween sub-faults, and V is the relative plate velocity of pl The frictional stress on the plate interface is expressed the target plate interface, assumed to be constant re- as the product of the frictional coefficient μ and effective gardless of sub-faults. V and G represent shear-wave normal stress σ . Changes in fault frictional stress eff velocity and rigidity of the material containing the fault, (d(μσ )/dt) are assumed to balance changes in shear eff respectively. The first term on the right-hand side in Eq. stress (Eq. (1)). We then obtain simultaneous differential (1) denotes the contribution of slip deficit rate distribu- equations with respect to V and θ for 0 ≤ i ≤ N. Under i i tion to the change in shear stress. It should be noted the hypothetical distribution of frictional parameters a , that since only perturbation from the steady-state sub- b ,and L for 0 ≤ i ≤ N, we can simulate the spatiotempo- i i duction is treated, shear stress of sub-fault i is designed ral slip behavior on the target plate interface by integrat- so as not to be affected by sub-fault j if the j-th sub-fault ing the simultaneous differential equations from a is slipping with the same slip velocity as relative plate particular initial condition (and boundary condition if motion, V . Further, because the static slip response pl required). function K is used in the first term, the radiation damp- ij ing term, representing the effect of wave propagation, is Configuration of model plate interface and its discretization added as the second term to better approximate the In order to reproduce the complicated slip behavior dur- stress state through the entire seismic cycle, including ing an earthquake cycle, we need to set appropriate fric- the coseismic stage. tional parameters along the megathrust. To obtain a As previously mentioned, we assume that fault friction suitable set of parameters, we require many computa- governed by RSF resists fault motion during an earth- tions with trial parameter settings. To reduce the com- quake cycle. The friction coefficient μ following RSF is putational cost for the seismic cycle simulation generally defined as, including multiple earthquakes with a larger gap in mag- nitudes compared with the previous studies (Hori and V ðÞ t V θ ðÞ t i  i Miyazaki, 2010, 2011), we conducted 2D simulations in- μ ðÞ t ¼ μ þ a In þ b In i i V L stead of conventional 3D ones. In the 2D simulation, we ð2Þ modeled the plate boundary fault with a set of line seg- ments, as the vertical cross section perpendicular to the The first term μ is a reference friction coefficient at Nankai Trough axis off the Kii Peninsula. The configur- the steady-state sliding with the reference slip rate V , ation of the upper interface of the subducting PSP was while the second and third terms of the right-hand side simply assumed to be a smooth curve with an average on Eq. (2) represent contributions of slip rate (i.e., V) dip angle of about 12° (Fig. 6a). For this vertical cross and state variable (i.e., θ) to the friction coefficient, re- section, we discretized the curved plate interface into spectively. As the state variable θ is an indicator of the 4400 line segments of nearly equal lengths of about contact state at the fault interface, the third term repre- 50 m. We confirmed that this length of line segment is sents the dependency of the frictional coefficient on sufficient to accurately model earthquake cycles with the “contact” state. Several versions have been proposed for distribution of frictional parameters shown later. Each the evolution of the state variable θ to describe state line segment corresponds to a sub-fault; the slip re- variable-dependency on μ. We adopt a widely used ver- sponse function K for this type of sub-fault was calcu- ij sion, known as an “aging law”, given by, lated by the analytic expression given by Rani and Singh (1992), which determines the stress change under the dθ ðÞ t V ðÞ t θ ðÞ t i i i ¼ 1 − ð3Þ plane strain field due to the unit dip slip of each dt L sub-fault in an elastic half space consisting of a Poisson In this study, we simply set V = V . Parameters a, b, solid with G = 30 GPa and V = 3273 m/s. We assumed ∗ pl s and L in Eqs. (2) and (3) are called frictional parameters. V = 5 cm/year for the plate subduction velocity. pl As shown in Eq. (2), a and b are parameters concerning rate- and state-dependencies of friction, respectively. L Distribution of frictional parameters represents a characteristic slip distance necessary for the Colors on the model plate interface in Fig. 6a reflect the friction μ to decrease significantly. In other words, the assumed distribution of frictional parameters, (a − b)σ : eff lower the value of L assumed, the larger the decreasing Red and blue colors denote stable ((a − b)σ > 0) and eff Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 15 of 17 unstable ((a − b)σ < 0) regions, respectively. The first assumed that the shallowest interface (depth > eff blue-colored region roughly corresponds to the depth of 10 km) exhibits stable sliding with (a − b)σ > 0; how- eff the source region with a large slip in simulated Nankai ever, the stable sliding condition did not place a suffi- Trough megathrust earthquakes and is thus termed a cient load on the shallower asperity to trigger an M6 “hyper-asperity.” In the hyper-asperity, we identify two earthquake. Therefore, we added the steady sliding zone segments with L values smaller than those of the other at the top of the plate interface as a boundary condition region (dark blue). The shallower one, identified at a (gray segment in Fig. 6a, depth from 5 to 6 km), where depth of 10–11 km with a length along the plate bound- the plate interface was assumed to be slipping with a ary fault of 10 km (“shallower asperity” hereafter), causes constant slip velocity of plate subduction V . The appli- pl M6 earthquakes. The deeper one is defined at a depth of cation of this boundary condition is realized by exclud- 13–17 km, with a length of 30 km along the plate ing the upper-most sub-faults from the numerical boundary fault (“deeper asperity” hereafter). The deeper computation, because our model only considers the per- asperity is set to model the high slip deficit rate patch turbation from steady sliding on the plate interface. deduced from the inversion analysis of the slip deficit After many trial simulations with the above boundary rate distribution reported by Yokota et al. (2016) condition, we found that the frictional distribution (Additional file 2: Figure S2). The slip deficit rate during shown in Fig. 6b provides the scenario most consistent the interseismic period may reflect a degree of interplate with the above two constraints. In this study, we only coupling. If so, parts with high slip deficit rate may re- present the scenario obtained for the frictional parame- lease larger coseismic slip than parts with relatively low ters shown in Fig. 6b, which has implications for the slip deficit rate. Such heterogeneity in the slip deficit rate Nankai Trough megathrust earthquake cycles. (and possible resultant variation in the coseismic slip distribution) can be modeled by the heterogeneous dis- Additional files tribution of L in the seismic cycle simulation with RSF. Under the constant (a − b)σ condition, small L values Additional file 1: Figure S1. Names of DONET1 stations and nodes. eff Gray contours represent depth distribution of the top of oceanic crust generate a higher slip deficit rate and a larger coseismic layer 2 (Citak et al. 2012). (PDF 1118 kb) slip than large L values. Therefore, we assume a smaller Additional file 2: Figure S2. Comparison of hypocenter distribution L value at the deeper asperity than for the rest of the with the slip deficit rate (SDR). Red circles represent hypocenters hyper-asperity (light-blue region). Similarly, the shal- obtained in this study. SDR is after Yokota et al. (2016). (PDF 184 kb) lower asperity is also set to model the source of the Additional file 3: Figure S3. A schematic image of 2.5-D velocity model used for hypocenter determination. (PDF 165 kb) off-Mie earthquake with a smaller L value. Additional file 4: Supplementary note. (DOCX 35 kb) As proposed by Hori and Miyazaki (2010, 2011), it is Additional file 5: Figure S4. Dependence of hypocenter locations on known that the existence of a smaller asperity, defined (a) stations used, (b) use of S-wave first arrivals, and (c) strike direction of by a region with a much smaller L value relative to its PS plate. Values in brackets represent the RMS travel time residual. Plots surroundings, within the hyper-asperity tends to cause are along the line A–B shown in Fig. 2a. (PDF 82 kb) hierarchical rupture behavior. During the cycle of earth- quakes rupturing the entire hyper-asperity, smaller Abbreviations earthquakes rupturing only the smaller asperity can CMT: Centroid moment tensor; DONET: Dense Oceanfloor Network System for Earthquakes and Tsunamis; MCS: Multi-channel seismic survey; occur without propagating to the surroundings. Accord- OBS: Ocean-bottom seismometers; PSP: Philippine Sea plate; RMS: Root- ingly, both shallower and deeper asperities with smaller mean-square; RSF: Rate- and state-dependent friction; SSE: Slow slip event; L values in our model can produce the complex seismic TOC: Top of the oceanic crust layer 2; TWT: Two-way travel time; VLFE: Very- low-frequency earthquake cycle behavior of a hyper-asperity, depending on the as- sumed gap in L. Here, we simulate earthquake cycle sce- Acknowledgements narios primarily by modifying the gap in L values in the We appreciate comments from Dr. L. Wallace and an anonymous reviewer, hyper-asperity so that the following two behaviors are which greatly improved the manuscript. We also appreciate Dr. Y. Yamamoto for sharing data for plotting Fig. 2. The Japan Meteorological Agency unified satisfied: (i) M8-class earthquake that ruptures the earthquake catalogue was used to plot Fig. 1. All figures were drawn using hyper-asperity only occurs when the deeper asperity be- Generic Mapping Tools (Wessel and Smith 1998). gins to rupture and (ii) the shallower asperity only rup- tures a part of the hyper-asperity and its postseismic slip Authors’ contributions MN conducted the hypocenter determinations and proposed the model. MH (afterslip/slow slip) mainly occurs at the shallow part of performed the numerical simulations of the hierarchical asperity model. MN the shallower asperity. and MH wrote the paper. AN and MY processed the structural data from To satisfy the latter constraint (ii), the shallow side of MCS reflection surveys. TH and TT organized the geophysical interpretations. SKa and KS participated in the data processing for hypocenter determinations. the plate interface must gradually reduce the slip deficit SKo, NT, and YK organized the data acquisitions and the project. All authors during an M8 earthquake cycle so that shear stress on equally contributed to scientific discussions. All authors read and approved the the shallower asperity accumulates. To achieve this, we final manuscript. Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 16 of 17 Competing interests Hori T, Miyazaki S, Mitsui N (2009) A model of earthquake-generation cycle with The authors declare that they have no competing interests. scale-dependent frictional property—preliminary results and research plan for a project of evaluation for coming Tokai, Tonankai, and Nankai earthquakes. J Disaster Res 4:111–117. Publisher’sNote Hyodo M, Hori T (2013) Re-examination of possible great interplate earthquake Springer Nature remains neutral with regard to jurisdictional claims in published scenarios in the Nankai Trough, southwest Japan, based on recent findings maps and institutional affiliations. and numerical simulations. Tectonophys 600:175–186. https://doi.org/10. 1016/j.tecto.2013.02.038. Author details Ichinose GA, Thio HK, Somerville PG, Sato T, Ishii T (2003) Rupture process of the R&D Center for Earthquake and Tsunami, Japan Agency for Marine-Earth 1944 Tonankai earthquake (Ms 8.1) from the inversion of teleseismic and Science and Technology, Yokohama, Japan. Earthquake and Volcano regional seismograms. J Geophys Res 108:2497. https://doi.org/10.1029/ Research Unit, National Research Institute for Earth Science and Disaster 2003JB002393. Resilience, Tsukuba, Japan. Institute of Education, Research and Regional Ishibashi K (2004) Status of historical seismology in Japan. Ann Geophys 47: Cooperation for Crisis Management Shikoku, Kagawa University, Takamatsu, 339–368. Japan. Kamei R, Pratt RG, Tsuji T (2012) Waveform tomography imaging of a megasplay fault system in the seismogenic Nankai subduction zone. Earth Planet Sci Received: 28 December 2017 Accepted: 14 May 2018 Lett 317-318:343–353. https://doi.org/10.1016/j.epsl.2011.10.042. Kamei R, Pratt RG, Tsuji T (2013) On acoustic waveform tomography of wide- angle OBS data—strategies for pre-conditioning and inversion. Geophys J Int References 194:1250–1280. https://doi.org/10.1093/gji/ggt165. Akuhara T, Mochizuki K (2014) Application of cluster analysis based on waveform Kaneda Y, Kawaguchi K, Araki E, Matsumoto H, Nakamura T, Kamiya S, Ariyoshi K, cross-correlation coefficients to data recorded by ocean-bottom Hori T, Baba T, Takahashi N (2015) Development and application of an seismometers: results from off the Kii Peninsula. Earth Planets Space 66:80. advanced ocean floor network system for megathrust earthquakes and https://doi.org/10.1186/1880-5981-66-80. tsunamis. In: Favali P et al (eds) Seafloor observatories. Springer, Heidelberg, Akuhara T, Mochizuki K, Nakahigashi K, Yamada T, Shinohara M, Sakai S, pp 643–662. https://doi.org/10.1007/978-3-642-11374-1_252. Kanazawa T, Uehira K, Shimizu H (2013) Segmentation of the Vp/Vs ratio and Kaneko L, Ide S, Nakano M (2018) Slow earthquakes in the microseism frequency low-frequency earthquake distribution around the fault boundary of the band (0.1–1.0 Hz) off Kii Peninsula, Japan. Geophys Res Lett 45. https://doi. Tonankai and Nankai earthquakes. Geophys Res Lett 40:1306–1310. https:// org/10.1002/2017GL076773. doi.org/10.1002/grl.50223. Kawaguchi K, Kaneko S, Nishida T, Komine T (2015) Construction of the DONET Ando M (1975) Source mechanisms and tectonic significance of historical real-time seafloor observatory for earthquakes and tsunami monitoring. In: earthquakes along the Nankai Trough, Japan. Tectonophysics 27:119–140. Favali P et al (eds) Seafloor observatories. Springer, Heidelberg, pp 211–228. Annoura S, Hashimoto T, Kamaya N, Katsumata A (2017) Shallow episodic tremor https://doi.org/10.1007/978-3-642-11374-1_10. near the Nankai Trough axis off southeast Mie prefecture, Japan. Geophys Kisslinger C, Jones LM (1991) Properties of aftershock sequences in southern Res Lett 44:3564–3571. https://doi.org/10.1002/2017GL073006. California. J Geophys Res 96:11,947–11,958. Araki E, Saffer DM, Kopf AJ, Wallace LM, Kimura T, Machida Y, Ide S, Davis E, IODP Kodaira S, Hori T, Ito A, Miura S, Fujie G, Park J-O, Baba T, Sakaguchi H, Kaneda Y Expedition 365 shipboard scientists (2017) Recurring and triggered slow-slip (2006) A cause of rupture segmentation and synchronization in the Nankai events near the trench at the Nankai Trough subduction megathrust. trough revealed by seismic imaging and numerical simulation. J Geophys Science 356:1157–1160. https://doi.org/10.1126/science.aan3120. Res 111:B09301. https://doi.org/10.1029/2005JB004030. Ariyoshi K, Nakata R, MatsuzawaT,HinoR,HoriT, HasegawaA,Kaneda Y Lomax A, Virieux J, Volant P, Berge C (2000) Probabilistic earthquake location in (2014) The detectability of shallowslowearthquakes by theDense 3D and layered models: introduction of a Metropolis–Gibbs method and Oceanfloor Network system for Earthquakes and Tsunamis (DONET) in comparison with linear locations. In: Thurber CH, Rabinowitz N (eds) Tonankai district, Japan. Mar Geophys Res 35:295–310. https://doi.org/10. Advances in seismic event location, Kluwer, Amsterdam, The Netherlands, pp 1007/s11001-013-9192-6. 101–134. Bangs NLB, Moore GF, Gulick SPS, Pangborn EM, Tobin HJ, Kuramoto S, Taira A Matsuzawa T, Hirose H, Shibazaki B, Obara K (2010) Modeling short- and long- (2009) Broad, weak regions of the Nankai Megathrust and implications for term slow slip events in the seismic cycles of large subduction earthquakes. J shallow coseismic slip. Earth Planet Sci Lett 284:44–49. https://doi.org/10. Geophys Res 115:B12301. https://doi.org/10.1029/2010JB007566. 1016/j.epsl.2009.04.026. Mikumo T, Miyatake T (1979) Earthquake sequences on a frictional fault model Bird P (2003) An updated digital model of plate boundaries. Geochem Geophys with non-uniform strengths and relaxation times. Geophys J R Astr Soc 59: Geosyst 4:1027. https://doi.org/10.1029/2001GC000252. 497–522. Citak SO, Nakamura T, Nakanishi A, Yamamoto Y, Ohori M, Baba T, Kaneda Y Mochizuki K, Fujie G, Sato T, Kasahara J, Hino R, Shinohara M, Suyehiro K (1998) (2012) An updated model of three-dimensional seismic structure in the Heterogeneous crustal structure across a seismic block boundary along the source area of the Tokai–Tonankai–Nankai earthquake. In abstract of AOGS– Nankai Trough. Geophys Res Lett 25:2301–2304. AGU (WPGM) Joint Assembly, Singapore, 13–17 august 2012, abstract no. Mochizuki K, Nakahigashi K, Kuwano A, Yamada T, Shinohara M, Sakai S, OS06-A015. Kanazawa T, Uehira K, Shimizu H (2010) Seismic characteristics around the DeMets C, Gordon RG, Argus DF (2010) Geologically current plate motions. fault segment boundary of historical great earthquakes along the Nankai Geophys J Int 181:1–80. https://doi.org/10.1111/j.1365-246X.2009.04491.x. Trough revealed by repeated long-term OBS observations. Geophys Res Lett Furumura T, Hayakawa T, Nakamura M, Koketsu K, Baba T (2008) Development of 37:L09304. https://doi.org/10.1029/2010GL042935. long-period ground motions from the Nankai Trough, Japan, earthquakes: Moore GF, Bangs NL, Taira A, Kuramoto S, Pangborn E, Tobin HJ (2007) Three- observations and computer simulation of the 1944 Tonankai (M 8.1) and dimensional splay fault geometry and implications for tsunami generation. the 2004 SE Off-Kii Peninsula (M 7.4) earthquakes. Pure Appl Geophys 165: Science 318:1128–1131. https://doi.org/10.1126/science.1147195. 585–607. https://doi.org/10.1007/s00024-008-0318-8. Murotani S, Shimazaki K, Koketsu K (2015) Rupture process of the 1946 Nankai Hori T (2006) Mechanisms of separation of rupture area and variation in time earthquake estimated using seismic waveforms and geodetic data. J interval and size of great earthquakes along the Nankai Trough, southwest Geophys Res 120:5677–5692. https://doi.org/10.1002/2014JB011676. Japan. J Earth Sim 5:8–19. Nakanishi A, Kodaira S, Miura S, Ito A, Sato T, Park J-O, Kido Y, Kaneda Y (2008) Hori T, Miyazaki S (2010) Hierarchical asperity model for multiscale characteristic Detailed structural image around splay-fault branching in the Nankai earthquakes: a numerical study for the off Kamaishi earthquake sequence in subduction seismogenic zone: results from a high-density ocean bottom the NE Japan subduction zone. Geophys Res Lett 37:L10304. https://doi.org/ seismic survey. J Geophys Res 113:B03105. https://doi.org/10.1029/ 10.1029/2010GL042669. 2007JB004974. Hori T, Miyazaki S (2011) A possible mechanism of M9 earthquake generation cycles in the area of repeating M7∼8 earthquakes surrounded Nakanishi A, Shiobara H, Hino R, Kodaira S, Kanazawa T, Shimamura H (1998) by aseismic sliding. Earth Planets Space 63:773–777. https://doi.org/10. Detailed subduction structure across the eastern Nankai Trough obtained 5047/eps.2011.06.022. from ocean bottom seismographic profiles. J Geophys Res 103:27,151–27,168. Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 17 of 17 Nakano M, Hori T, Araki E, Kodaira S, Ide S (2018) Shallow very-low-frequency geology. Elsevier BV, pp 599–640. https://doi.org/10.1016/B978-0-444-62617- earthquakes accompany slow slip events in the Nankai subduction zone. Nat 2.00020-7. Commun 9:984. https://doi.org/10.1038/s41467-018-03431-5. Tsuji T, Ashi J, Strasser M, Kimura G (2015) Identification of the static Nakano M, Nakamura T, Kamiya S, Kaneda Y (2014) Seismic activity beneath the backstop and its influence on the evolution of the accretionary prism in Nankai trough revealed by DONET ocean-bottom observations. Mar Geophys the Nankai Trough. Earth Planet Sci Lett 431:15–25. https://doi.org/10. Res 35:271–284. https://doi.org/10.1007/s11001-013-9195-3. 1016/j.epsl.2015.09.011. Tsuji T, Kamei R, Pratt RG (2014) Pore pressure distribution of a mega-splay fault Nakano M, Nakamura T, Kamiya S, Ohori M, Kaneda Y (2013) Intensive seismic activity system in the Nankai Trough subduction zone: insight into up-dip extent of around the Nankai trough revealed by DONET ocean-floor seismic observations. the seismogenic zone. Earth Planet Sci Lett 396:165–178. https://doi.org/10. Earth Planets Space 65:5–15. https://doi.org/10.5047/eps.2012.05.013. 1016/j.epsl.2014.04.011. Nakano M, Nakamura T, Kaneda Y (2015) Hypocenters in the Nankai Trough Tsuji T, Minato S, Kamei R, Tsuru T, Kimura G (2017) 3D geometry of a plate determined by using data from both ocean-bottom and land seismic boundary fault related to the 2016 Off-Mie earthquake in the Nankai networks and a 3D velocity structure model: implications for seismotectonic subduction zone, Japan. Earth Planet Sci Lett 478:234–244. https://doi.org/10. activity. Bull Seism Soc Am 105:1594–1605. https://doi.org/10.1785/ 1016/j.epsl.2017.08.041. Tsuru T, Miura S, Park J-O, Ito A, Fujie G, Kaneda Y, No T, Katayama T, Kasahara J Noda H, Nakatani M, Hori T (2013) Large nucleation before large earthquakes is (2005) Variation of physical properties beneath a fault observed by a two- sometimes skipped due to cascade-up—implications from a rate and state ship seismic survey off Southwest Japan. J Geophys Res 110:B05405. https:// simulation of faults with hierarchical asperities. J Geophys Res 118:2924– doi.org/10.1029/2004JB003036. 2952. https://doi.org/10.1002/jgrb.50211. Utsu T (1961) A statistical study on the occurrence of aftershocks. Geophys Mag Noda H, Nakatani M, Hori T (2014) Coseismic visibility of a small fragile patch 30:521–605. involved in the rupture of a large patch—implications from fully dynamic Wallace LM, Araki E, Saffer D, Wang X, Roesner A, Kopf A, Nakanishi A, Power W, rate-state earthquake sequence simulations producing variable manners of Kobayashi R, Kinoshita C, Toczko S, Kimura T, Machida Y, Carr S (2016) Near- earthquake initiation. Progress Earth Planetary Sci 1:8. https://doi.org/10.1186/ field observations of an offshore M 6.0 earthquake from an integrated 2197-4284-1-8. w seafloor and subseafloor monitoring network at the Nankai Trough, Obana K, Kodaira S (2009) Low-frequency tremors associated with reverse faults southwest Japan. J Geophys Res 121:8338–8351. https://doi.org/10.1002/ in a shallow accretionary prism. Earth Planet Sci Lett 287:168–174. https://doi. 2016JB013417. org/10.1016/j.epsl.2009.08.005. Wessel P, Smith WHF (1998) New, improved version of generic mapping tools Obana K, Kodaira S, Kaneda Y (2004) Microseismicity around rupture area of the released. Eos 79:579. 1944 Tonankai earthquake from ocean bottom seismograph observations. Yamamoto Y, Takahashi T, Kaiho Y, Obana K, Nakanishi A, Kodaira S, Kaneda Y Earth Planet Sci Lett 222:561–572. https://doi.org/10.1016/j.epsl.2004.02.032. (2017) Seismic structure off the Kii Peninsula, Japan, deduced from passive- Obana K, Kodaira S, Kaneda Y (2005) Seismicity in the incoming/subducting and active-source seismographic data. Earth Planet Science Lett 461:163–175. Philippine Sea plate off the Kii Peninsula, central Nankai trough. J Geophys https://doi.org/10.1016/j.epsl.2017.01.003. Res 110:B11311. https://doi.org/10.1029/2004JB003487. Yokota Y, Ishikawa T, Watanabe S, Tashiro T, Asada A (2016) Seafloor geodetic Park J-O, Fujie G, Wijerathne L, Hori T, Kodaira S, Fukao Y, Moore GF, Bangs NL, constrains on interplate coupling of the Nankai Trough megathrust zone. Kuramoto S, Taira A (2010) A low-velocity zone with weak reflectivity along Nature 534:374–377. https://doi.org/10.1038/nature17632. the Nankai subduction zone. Geology 38:283–286. https://doi.org/10.1130/ G30205.1. Park J-O, Kodaira S (2012) Seismic reflection and bathymetric evidences for the Nankai earthquake rupture across a stable segment-boundary. Earth Planets Space 64:299–303. https://doi.org/10.5047/eps.2011.10.006. Park J-O, Moore GF, Tsuru T, Kodaira S, Kaneda Y (2003) A subducted oceanic ridge influencing the Nankai megathrust earthquake rupture. Earth Planet Sci Lett 217:77–84. https://doi.org/10.1016/S0012-821X(03)00553-3. Park J-O, Tsuru T, Kodaira S, Cummins PR, Kaneda Y (2002) Splay fault branching along the Nankai subduction zone. Science 297:1157–1160. https://doi.org/ 10.1126/science.1074111. Podvin P, Lecomte I (1991) Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools. Geophys J Int 105:271–284. Rani S, Singh SJ (1992) Static deformation of a uniform half-space due to a long dip-slip fault. Geophys J Int 109:469–476. Rice JR (1993) Spatio-temporal complexity of slip on a fault. J Geophys Res 98: 9885–9907. Sakaguchi A, Chester F, Curewitz D, Fabbri O, Goldsby D, Kimura G, Li C-F, Masaki Y, Screaton EJ, Tsutsumi A, Ujiie K, Yamaguchi A (2011) Seismic slip propagation to the updip end of plate boundary subduction interface faults: vitrinite reflectance geothermometry on Integrated Ocean Drilling Program NanTro SEIZE cores. Geology 39:395–398. https://doi.org/10.1130/G31642.1. Sakai S, Yamada T, Shinohara M, Hagiwara H, Kanazawa T, Obana K, Kodaira S, Kaneda Y (2005) Urgent aftershock observation of the 2004 off the Kii Peninsula earthquake using ocean bottom seismometers. Earth Planets Space 57:363–368. Seno T (2005) The September 5, 2004 off the Kii Peninsula earthquakes as a composition of bending and collision. Earth Planets Space 57:327–332. Takemura S, Shiomi K, Kimura T, Saito T (2016) Systematic difference between first-motion and waveform-inversion solutions for shallow offshore earthquakes due to a low-angle dipping slab. Earth Planets Space 68:149. https://doi.org/10.1186/s40623-016-0527-9. Takenaka H, Kennett BLN (1996) A 2.5-D time-domain elastodynamic equation for plane-wave incidence. Geophys J Int 125:F5–F9. Tobin H, Henry P, Vannucchi P, Screaton E (2014) Subduction zones: structure and deformation history. In: Stein R et al (eds) Developments in marine http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Progress in Earth and Planetary Science Springer Journals

The 2016 Mw 5.9 earthquake off the southeastern coast of Mie Prefecture as an indicator of preparatory processes of the next Nankai Trough megathrust earthquake

Free
17 pages
Loading next page...
 
/lp/springer_journal/the-2016-mw-5-9-earthquake-off-the-southeastern-coast-of-mie-Rq8LPz0e9y
Publisher
Springer Berlin Heidelberg
Copyright
Copyright © 2018 by The Author(s).
Subject
Earth Sciences; Earth Sciences, general; Geophysics/Geodesy; Planetology; Biogeosciences; Hydrogeology; Atmospheric Sciences
eISSN
2197-4284
D.O.I.
10.1186/s40645-018-0188-3
Publisher site
See Article on Publisher Site

Abstract

Megathrust earthquakes have occurred repeatedly at intervals of 100 to 150 years along the Nankai Trough, situated in the southwest of Japan. Given that it has been 70 years since the last event, the occurrence of the next devastating earthquake is anticipated in the near future. On April 1, 2016, a moderate earthquake (M 5.9, M 6.5) occurred off the southeastern coast of Mie Prefecture in the source region of the 1944 Tonankai JMA earthquake (M 8.2). In this study, we investigated the influence of the 2016 earthquake on future megathrust earthquakes. We first determined the hypocenter distributions using a precise velocity structure obtained from seismic surveys in the source region. Using data obtained from the DONET ocean-bottom observation network, we found that this earthquake occurred along the plate boundary fault, which is also believed to have slipped during past megathrust earthquakes. We then performed a preliminary numerical simulation to reproduce the occurrence of a moderate earthquake in the middle of a megathrust earthquake cycle. We used a hierarchical asperity model, in which smaller asperities causing moderate earthquakes are embedded in a hyper-asperity that serves as the source region of megathrust earthquakes. The simulation shows that moderate earthquakes, caused by ruptures of a smaller asperity, occur as a result of shrinkage of strongly coupled areas in the hyper-asperity. This result is consistent with the observation that the hypocenter of the 2016 earthquake was located at the edge of a strongly coupled region along the plate boundary. The simulation also reproduced postseismic slip (afterslip/slow slip) along the plate boundary updip of the hyper-asperity, which is consistent with the observations of slow slip events and very-low-frequency earthquakes after the 2016 earthquake. Thus, the occurrence of the moderate earthquake offshore southeastern Mie Prefecture in the middle of the megathrust earthquake cycle implies that the shrinkage of the strongly coupled area along the plate boundary is occurring as a preparatory process for the next megathrust earthquake in the region. Keywords: Hypocenter determination, Crustal structure, Ocean floor observations, Accretionary prism, Earthquake cycle simulation, Hierarchical asperity model * Correspondence: mnakano@jamstec.go.jp R&D Center for Earthquake and Tsunami, Japan Agency for Marine-Earth Science and Technology, Yokohama, Japan Full list of author information is available at the end of the article © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 2 of 17 Introduction seismic activities (Kaneda et al. 2015; Kawaguchi et al. At convergent plate boundaries, large earthquakes occur 2015). These observations revealed that, in recent years, repeatedly with intervals of one to several hundred years. earthquake activities along the Nankai Trough have The Nankai Trough, where the Philippine Sea plate mostly occurred in the oceanic crust or upper mantle of (PSP) subducts beneath the Amur plate (Bird 2003; the subducting PSP (Obana et al. 2004, 2005; Obana and DeMets et al. 2010), is one such region where mega- Kodaira 2009; Mochizuki et al. 2010; Akuhara et al. thrust earthquakes and tsunamis have occurred at inter- 2013; Nakano et al. 2013, 2014, 2015; Akuhara and vals of 100 to 150 years (e.g., Ando 1975; Ishibashi Mochizuki 2014). The most recent notable events during 2004). The 1944 Tonankai (M 8.2; Ichinose et al. 2003) this time have been the 2004 M 7.1, 7.4, and 6.5 w JMA and 1946 Nankai (M 8.4; Murotani et al. 2015) earth- earthquakes off the Kii Peninsula, which occurred in the quakes were two such devastating events, which rup- subducting plate (e.g., Sakai et al. 2005; Seno 2005). tured the eastern and western segments of the trough, Significant improvements in supercomputers have fa- respectively (Fig. 1), causing severe damage in western cilitated intensive numerical simulations of megathrust and central Japan. Future earthquakes with matching in- earthquakes (e.g., Hori 2006; Kodaira et al. 2006; Hori et tensities will almost certainly result in severe damage, al. 2009; Hyodo and Hori 2013). These studies have given the strong ground motion and tsunamis likely to mainly focused on reproducing the slip distributions, be associated with them. Recent studies have indicated rupture segmentations, and recurrence intervals of his- that such megathrust earthquakes would also cause torical earthquakes, based on realistic fault models and massive damage due to long-period ground motion in subduction zone structures obtained from seismic sur- distant regions, including metropolitan Tokyo (e.g., Fur- veys. In contrast, simulations based on simplified fault umura et al. 2008). models can provide us with insights into earthquake In order to clarify the deformation process along the rupture patterns, cycles, and related phenomena such as convergent margin and the potential for such disasters, slow slip events (SSEs) (e.g., Matsuzawa et al. 2010; interdisciplinary studies, including seismic surveys of Noda et al. 2013, 2014; Ariyoshi et al. 2014). Hori and subsurface structures, seismic monitoring, geodetic ob- Miyazaki (2011) simulated occurrences of M7–8 earth- servations, and numerical simulations, have been con- quakes during M9 earthquake cycles to reproduce the ducted. A number of seismic surveys have investigated occurrence of the 1978 M 7.5 earthquake that occurred geophysical and geological structures in detail (e.g., offshore Miyagi before the 2011 Tohoku-oki earthquake Mochizuki et al. 1998; Nakanishi et al. 1998, 2008; Park (M 9.0). They used the hierarchical asperity model, in and Kodaira 2012; Park et al. 2002, 2003, 2010; Kamei et which a smaller asperity causing M7−8 earthquakes was al. 2012). One of the most remarkable structures in the embedded in a hyper-asperity that caused the M9 earth- accretionary prism along the Nankai Trough is the quake; namely, the fundamental rupture mode in the megasplay fault or the out-of-sequence thrust branching subduction zone. Numerical simulations can also indi- from the decollement, imaged as strong reflections in cate the temporal evolution of coupling ratios and stress multi-channel seismic survey (MCS) profiles (Park et al. states along fault interfaces. These parameters are, how- 2002; Moore et al. 2007). Recent studies, based on de- ever, still very difficult to estimate from observations but tailed seismic surveys, have interpreted the megasplay are considered to be crucial for monitoring the prepara- fault as the main plate boundary fault (Bangs et al. 2009; tory processes of megathrust earthquakes (e.g., Kodaira Kamei et al. 2012, 2013; Tobin et al. 2014; Tsuji et al. et al. 2006). Thus, combinations of numerical simula- 2014). The megasplay fault branches in the shallower tions with seismic and geodetic observations could pro- part; one section reaches the ocean bottom and the vide information about the current state of stress and other connects to the decollement beneath the toe of strain along plate boundary faults. the accretionary prism (Tsuji et al. 2014). Both branches On April 1, 2016, an earthquake occurred off the are considered to have slipped during past megathrust southeastern coast of Mie Prefecture (M = 5.9 as per earthquakes (Sakaguchi et al. 2011; Tsuji et al. 2014). the U.S. Geological Survey) (“off-Mie earthquake” here- Extensive observations of seismic activities along the after). It occurred updip of the hypocenter of the 1944 Nankai Trough have been conducted using temporary Tonankai earthquake (Fig. 1a). The hypocenter was be- deployments of recoverable ocean-bottom seismometers neath the DONET stations, and real-time ocean-bottom (OBS) (e.g., Obana et al. 2004; Obana and Kodaira 2009; observations were obtained. The epicenter was within Mochizuki et al. 2010). The recent development of a the slip area of the 1944 Tonankai earthquake (Ichinose permanent seafloor observation network called the et al. 2003) and close to the eastern margin of the 1946 Dense Oceanfloor Network System for Earthquakes and Nankai earthquake slip area (Murotani et al. 2015). Dur- Tsunamis (DONET) (Fig. 1 and Additional file 1: Figure ing the 70-year period since the occurrence of the last S1) has facilitated continuous monitoring of offshore megathrust earthquake, earthquakes with magnitudes Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 3 of 17 Fig. 1 (See legend on next page.) Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 4 of 17 (See figure on previous page.) Fig. 1 Past seismic activities in the study area. a Regional location map (inset) and seismic activity along the Nankai Trough. Orange star represents the epicenter location of the 2016 off-Mie earthquake; results from the Japan Meteorological Agency (JMA), U.S. Geological Survey (USGS), and the Global CMT Project (GCMT) are shown together with mechanism symbols. Epicenters of earthquakes larger than M6 before and after 1976 are represented by diamonds and circles, respectively, from the JMA unified earthquake catalogue. Blue and red contours represent the slip distributions of the 1944 Tonankai and 1946 Nankai earthquakes, respectively (Ichinose et al. 2003;Murotani et al. 2015). Gray triangles represent locations of DONET stations. The ellipsoid labeled with “A” designates the region of interplate seismic activity reported by Akuhara and Mochizuki (2014). Orange arrows represent the subduction direction of the Philippine Sea plate. Plate name abbreviations are as follows: AM, Amur; OK, Okhotsk; PA, Pacific; PS, Philippine Sea. b Occurrence date of earthquakes plotted against their source longitude. Symbols are the same as in a larger than 6 have been very rare along the eastern Nan- that the sources were distributed very close to the plate kai Trough (Fig. 1b). Exceptions were aftershocks of the interface. However, strong horizontal heterogeneities of 1944 and 1946 earthquakes that continued for about seismic wave velocities due to plate subduction could 10–15 years and the sequence of the 2004 intraplate affect the precision of determining the hypocenter loca- earthquakes off the Kii Peninsula. Hence, large interplate tions. Takemura et al. (2016) showed that an accurate earthquakes have not been known to occur until the 3D velocity model reflecting the plate subduction is middle of a megathrust earthquake cycle. Wallace et al. necessary to explain the spatial distributions of P-wave (2016) consider that the 2016 off-Mie earthquake most first motion polarities observed for this earthquake. likely occurred on the plate interface, based on inte- Therefore, a velocity model reflecting the horizontal grated analyses of hypocenter distributions, seafloor heterogeneities is crucial for precise determination of crustal deformations, and tsunami modeling using sea- hypocenters. floor geodetic and seismological data from DONET. In this study, we determined the hypocenters of the After this earthquake, SSEs along the plate boundary, 2016 off-Mie earthquake and its aftershocks using a vel- very-low-frequency earthquakes (VLFEs), and tectonic ocity model that includes horizontal heterogeneities tremors have been observed updip of the hypocentral re- along the subduction zone. We used a velocity structure gion (Fig. 2) (Annoura et al. 2017; Araki et al. 2017; obtained from a wide-angle OBS survey along a line Kaneko et al. 2018; Nakano et al. 2018). Recent studies passing through the source region shown by Wallace et have shown that signals from these slow events radiated al. (2016). We compared the hypocenter distribution from a common source (Kaneko et al. 2018; Nakano et with the MCS reflection profiles to confirm whether this al. 2018). The development of a comprehensive model of earthquake occurred on the plate boundary. Since this earthquake sources that can explain the generation of earthquake was located along the plate boundary fault, these regular and slow earthquakes would help us to we conducted a preliminary numerical simulation of understand the processes of stress accumulation and/or plate boundary earthquakes to understand the implica- release along the plate boundary. It is important to de- tions of this earthquake for the Nankai Trough mega- termine the precise source location of the 2016 off-Mie thrust earthquake cycle. Using the hierarchical asperity earthquake to evaluate its influence on the generation of model, we found that the 2016 off-Mie earthquake oc- the next megathrust earthquake and its preparatory curred due to shrinkage of the strongly coupled area processes. along the plate interface, which is likely to be a prepara- The centroid moment tensor (CMT) solution for the tory process for the next megathrust earthquake. Inter- mainshock shows a thrust-type mechanism consistent estingly, this model can also explain the absence of large with plate boundary slip (Fig. 1a). However, the source earthquakes until the middle of the megathrust earth- centroid location, which is based on regional or teleseis- quake cycle. mic data, is uncertain; depending on the institution that determined the location (Japan Meteorological Agency, Results and discussion U.S. Geological Survey, and the GlobalCMT project), Hypocenter distribution of the 2016 off-Mie earthquake there is a variation of ~ 20 km in the horizontal direc- and aftershocks tion and 7 km (between 12 and 19 km) in depth. Such a We determined the hypocenters of the mainshock and discrepancy may have occurred because of strong struc- aftershocks of the 2016 off-Mie earthquake using a local tural heterogeneities that exist in the dip direction of velocity model obtained in the source region (Fig. 2). plate subduction. Precise determination of the source lo- We used a horizontally heterogeneous velocity structure cation, especially for the depth, is crucial to clarify based on that obtained along a survey line that passes whether this earthquake occurred along the plate bound- through the hypocentral region (see the section ary. Wallace et al. (2016) determined the hypocenter dis- “Methods/experimental”). The hypocenters were distrib- tribution using a 1D velocity structure, which showed uted several kilometers shallower than those based on Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 5 of 17 Fig. 2 (See legend on next page.) Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 6 of 17 (See figure on previous page.) Fig. 2 Hypocenters of the off-Mie earthquake and aftershocks. a Hypocenter distribution of the 2016 off-Mie earthquake and aftershocks (red circles). Pink focal mechanism symbols indicate shallow VLFEs, and magenta dots are locations of tectonic tremors in April 2016 (Kaneko et al. 2018; Nakano et al. 2018). White triangles represent DONET stations used for our hypocenter relocation. Green and blue lines represent location of MCS survey lines shown in Fig. 5. The thin black line represents the location of MCS survey line (COP2), and the thick black segment represents the location of a topographic high (Tsuru et al. 2005). The purple dashed line indicates the boundary of the topographic high in the basin. Red and blue dashed contours represent regions of low- and high-velocity anomalies, respectively, in the overriding plate (Yamamoto et al. 2017). A high-velocity body off the Kii Peninsula is represented by an orange oval (Kodaira et al. 2006). b Velocity structure used for our hypocenter determination obtained along MCS survey line KI03 (Wallace et al. 2016) our routine hypocenter determinations using a 1D vel- but the depth difference from the plate boundary fault ocity structure (Fig. 3). The mainshock depth was 9.2 ± was smaller in our result. The hypocenter location corre- 0.6 km below sea level, about 4 km shallower than that sponds to the shallower edge of the strongly coupled re- obtained by the routine analysis, and about 2 km shal- gion along the plate boundary, as noted by Yokota et al. lower than that obtained by Wallace et al. (2016). The (2016) (Additional file 2: Figure S2). differences can be attributed to the horizontal heteroge- neities of the velocity structure considered in our ana- Earthquake cycle simulation based on the hierarchical lysis. Aftershocks were distributed in a very limited area, asperity model about 10 km downdip of the mainshock, at depths of be- In order to investigate the generating mechanism of the tween 10 and 14 km, with source location errors of off-Mie earthquake in the middle of the megathrust 1.9 km in both horizontal and depth directions. There is earthquake cycle along the Nankai Trough, we con- a seismic gap distinct from the mainshock hypocenter, ducted a preliminary earthquake cycle simulation based which agrees with the results obtained by Wallace et al. on the hierarchical asperity model (Hori and Miyazaki (2016). The aftershock activity continued for 2 weeks, 2010, 2011) that has been used previously to simulate and a burst of reactivation occurred on April 19 (Fig. 4), multiscale earthquakes. In this model, smaller asperities, coinciding with the end of the SSE following the main- causing moderate-sized earthquakes, are embedded in a shock (Araki et al. 2017). larger asperity (“hyper-asperity”), causing the fundamen- The aftershock distribution can be further separated tal rupture mode of the fault system. To simulate the into two clusters: northeastern and southwestern generation of earthquakes on a megathrust, we solved (Fig. 3a). Events in the northeastern cluster were distrib- the quasi-dynamic equation of shear stress along fault uted almost parallel to the strike of the PSP, with a segments (Rice 1993) using the rate- and state- near-identical source depth, while those in the south- dependent friction (RSF) laws to represent frictional western cluster showed a rather scattered distribution in properties across the fault (see the section titled the dip direction. Event sizes were systematically larger “Methods/experimental”). The occurrence of an earth- in the northeastern cluster, and decay times of these ac- quake on the plate boundary is modeled as a perturb- tivities were clearly different. Activity in the northeastern ation from the steady-state relative plate motion. Shear cluster lasted for only a day, whereas, in the southwest- stress builds up due to secular motion of the subducting ern cluster, it continued for much longer (Fig. 4). plate which is finally released as earthquakes when the Figure 5 compares the earthquake depths and reflec- accumulated shear stress exceeds the fault strength gov- tors imaged from MCS surveys along lines passing erned by RSF. through the source region (see the section “Methods/ex- In our simulation, we attempted to reproduce the fol- perimental”). Clear reflections indicative of the plate lowing characteristics observed for the off-Mie earth- boundary are observed in the MCS profiles. The main- quake. First, the earthquake occurred after 70 years of a shock was located at a slightly shallower depth than the quiet period characterized by the absence of M6 or lar- upper reflector, possibly corresponding to the megasplay ger earthquakes after the previous M8 earthquake. Sec- fault. The aftershocks were distributed about 10 km ond, a slow slip and/or afterslip episode was observed landward from the mainshock, at depths about 2 km immediately after this earthquake (Wallace et al. 2016; deeper than the reflector (Fig. 3). This result agrees with Araki et al. 2017). In order to reproduce such compli- that of Wallace et al. (2016), although the separation cated slip behavior during an earthquake cycle, we need depth was comparable to the location errors. to properly define the frictional parameters along the Considering the location errors associated with deter- megathrust, as well as the sizes and locations of asper- mining the hypocenters and reflector depths (±1.1 km), ities. Suitable parameters need to be determined empir- we consider that this earthquake has slipped along the ically and the computational cost for conventional 3D plate boundary fault, as Wallace et al. (2016) reported, simulations is high; hence, we conducted 2D simulations Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 7 of 17 Fig. 3 (See legend on next page.) Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 8 of 17 (See figure on previous page.) Fig. 3 Comparison of routine and relocated hypocenter locations. Yellow and red circles represent routine and relocated hypocenters, respectively. Black and red crosses at the top left of each panel represent the uncertainties of hypocenter locations for the mainshock and aftershocks, respectively. a Map view. Gray triangles represent locations of DONET stations. b Cross-section along line E–F. Black and blue lines represent strong reflections picked from the MCS profile shown in Fig. 5a (temporal evolution of slip on a fault, represented by line updip of the hyper-asperity, and continues for more than segments) for the parameter search. In this study, we 20 days (Fig. 7b). After the occurrence of the M6 event, present one of the simulation results that reproduced the strongly coupled area in the hyper-asperity further the characteristics described above well. shrinks because of the secular aseismic slip in the sur- In the numerical simulation, the frictional parameters rounding area, and the shear stress of the deeper asper- of asperities are set as unstable slip conditions, while re- ity further increases (Fig. 8). When the stress at the gions outside of them are set as stable. Since the hypo- deeper asperity reaches the threshold, a rupture initiated center of the off-Mie earthquake was located at the at the asperity propagates along the entire hyper-asperity shallower edge of the strongly coupled region along the and an M8 earthquake occurs. Nankai Trough (Additional file 2: Figure S2), we defined Using earthquake cycle simulation based on the hier- a smaller asperity, causing an M6-class earthquake, at archical asperity model, we successfully reproduced the the shallower edge of the hyper-asperity. Figure 6 shows occurrence of a moderate earthquake in the middle of the fault geometry and frictional parameters, in which the megathrust earthquake cycle and the SSE (or after- we assume two smaller asperities (dark blue segments) slip) that follows the earthquake. The rupture of the embedded in the hyper-asperity (both dark blue and shallower asperity is caused by shrinkage of the strongly light-blue segments). Based on a simulation using this coupled area along the hyper-asperity, which results frictional parameter distribution, the rupture of a shal- from the secular plate subduction. This finding implies lower asperity placed on a 10-km-long line segment that the occurrence of the moderate 2016 off-Mie earth- along the plate boundary fault, corresponding to a 10– quake in the middle of the earthquake cycle represents 11 km depth range (“shallower asperity” hereafter), shrinkage of a strongly coupled area, which could be a caused the M6 earthquake. Whereas, the rupture of a preparatory process of the next megathrust earthquake deeper asperity, a 30-km-long segment along the plate along the Nankai Trough. We note here that the cycle of boundary fault corresponding to a 13–17 km depth the M8 earthquake and the timing and frequency of the range (“deeper asperity” hereafter), propagates along the M6 earthquake depend on the frictional parameters of entire hyper-asperity to cause an M8 earthquake, the the fault, their distributions, and the slip velocity at the fundamental rupture mode of this fault system. plate interface, all of which were assumed in the Figure 7a shows the spatiotemporal evolution of the simulation. slip rate along the fault. At the beginning of the M8 Aftershocks of the off-Mie earthquake occurred only megathrust earthquake cycle (marked by “M8” at the downdip of the mainshock, with a clear separation of bottom of Fig. 7a), the slip velocity on the fault is very about 10 km (similar to Wallace et al. 2016). The ab- small along the entire asperity, that is, the fault is sence of an updip aftershock is not surprising, because strongly coupled. Outside the hyper-asperity, a secular this region is considered to have conditionally stable aseismic slip occurs during the entire earthquake cycle frictional properties (Wallace et al. 2016). Instead, slow because of the stable sliding frictional properties. The earthquakes as SSEs, VLFEs, and tectonic tremors oc- strongly coupled area on the hyper-asperity gradually curred there, sharing a common updip source region shrinks from the deeper and shallower edges due to the (Araki et al. 2017; Kaneko et al. 2018; Nakano et al. secular aseismic slip, which causes a gradual increase in 2018) (Fig. 2). Both regular and slow earthquakes were shear stress on the shallower and deeper asperities almost entirely absent between the downdip aftershocks (Fig. 8). From the first half to the middle stage of the and updip slow earthquakes, although the occurrence of M8 megathrust earthquake cycle, the rate of increase of an SSE or an afterslip close to the hypocenter region stress is higher in the shallower asperity because of its was implied after the mainshock (Wallace et al. 2016; location closer to the edge of the hyper-asperity. In our Nakano et al. 2018). Our numerical simulation success- simulation, it takes about 120 years for sufficient stress fully reproduced the slow slip updip of the asperity, to accumulate to rupture the shallower asperity (M6 where stable sliding frictional parameters were defined. earthquake) in a 200-year megathrust earthquake cycle. Although the reason for the existence of a gap between The M6 earthquake is followed by postseismic slip the mainshock hypocenter and aftershocks is still elu- (afterslip/slow slip) that occurs mainly in a region of sive, Tsuji et al. (2017) considered that aftershocks oc- velocity-strengthening frictional parameters, that is, an curred close to the location where the plate boundary Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 9 of 17 0 0 7 14 21 Date of April 2016 Fig. 4 Cumulative number and magnitude of mainshock and aftershocks plotted against time. Circles with vertical bars represent the magnitude of the earthquake. The black line represents the cumulative number of aftershocks. Blue and red lines correspond to those occurring in the northeastern and southwestern clusters separated by a dashed line, shown in Fig. 3a decollement soles out into the top of oceanic crust layer in geological structures would cause differences in 2 (TOC). Such heterogeneous structures may cause frictional properties along aftershock faults. The after- stress accumulations, which could explain the very lim- shock region corresponds to the boundary between older ited aftershock distribution. Wallace et al. (2016) sug- (~ 14 Ma) and younger (~ 6 Ma) prisms (Tsuji et al. gested that this gap coincided with the mainshock fault 2017), which would lead to a difference in the degree of that did not trigger aftershocks because stress was heterogeneities in plate boundary structure. already released, or stress was released due to aseismic Previous studies showed that structural heterogeneities slow slip, or both. The migration of VLFE sources updip, exist in overriding and subducting plates around the away from the mainshock hypocenter, supports the latter source region of the off-Mie earthquake (Fig. 2). One of model (Nakano et al. 2018). the characteristic structures off the Kii Peninsula is a The aftershock activity was separated into two clusters, plutonic rock body with high seismic velocity and high which showed different decay times, namely, the p value density, which is coincident with the segment boundary of the Omori–Utsu aftershock decay law (Figs. 3 and 4). of megathrust earthquakes along the Nankai Trough −p The Omori–Utsu law is given as n(t)= K(t + c) , where (Kodaira et al. 2006). A topographic high in the TOC n(t) is the number of events per unit time at t since the can be recognized to the west of the aftershock area in mainshock; K, c, and p are constants (Utsu 1961). Al- the MCS reflection survey obtained along the strike dir- though source properties corresponding to the p value ection of the subducting plate (line COP2 in Tsuru et al. are not well studied, Kisslinger and Jones (1991) ob- (2005)). A boundary of the ocean floor topography (the served that high heat flow causes a high p value. They purple dashed line in Fig. 2a) is also coincident with the considered that the p value represents the stress relax- epicenter distribution of the aftershocks. This region ation time, and the aftershock decay rate is controlled by corresponds to the boundary of low- and high- velocity the temperature at the source depth. In our case, how- anomalies in the overriding plate (Yamamoto et al. ever, the distance of the aftershock clusters was only sev- 2017). The source of the off-Mie earthquake is close to eral kilometers and the temperature would not differ the boundary of the backstop and sediments (Tsuji et al. significantly. Mikumo and Miyatake (1979) conducted 2015) and also corresponds to the location where the de- numerical simulations and showed that the p value be- collement soles out into the TOC (Tsuji et al. 2017). comes smaller as the distribution of fault strength be- Tsuji et al. (2017) also pointed out that the aftershock comes more heterogeneous. Small-scale heterogeneities region corresponds to the boundary between the older Magnitude Cumulative Number Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 10 of 17 a b c d Fig. 5 Hypocenters compared with plate interfaces on sections from a MCS survey. Hypocenter distribution projected on the time-migrated MCS reflection profiles along lines A–B(a) and C–D(b) shown in Fig. 2a. Arrows represent strong reflections indicative of a plate boundary (PB) and younger prisms. These heterogeneous structures subducting or overriding plates, and there is very low ac- concentrate stress and limit the area of the aftershock tivity in the rupture area of the megathrust earthquakes activity. (Obana et al. 2004; Obana and Kodaira 2009; Nakano et In our earthquake cycle simulations, we successfully al. 2015). However, Akuhara and Mochizuki (2014) ob- reproduced the occurrence of a moderate earthquake in served continuous activity of an earthquake cluster with the middle of the megathrust earthquake cycle. The magnitudes smaller than 4 along the plate boundary simulation implies that the M6-class off-Mie earthquake west of the Kii Peninsula (labeled “A” in Fig. 1a). Their was caused by shrinkage of the strongly coupled area hypocenters roughly correspond to the downdip end of along the plate boundary. Present seismic activities in the rupture area of the 1946 Nankai earthquake (Muro- the eastern Nankai Trough are mostly within the tani et al. 2015). These observations imply that the Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 11 of 17 ab Fig. 6 Model setting of 2D earthquake cycle simulation. a Configuration of model plate interface and the schematic illustration of stable or unstable regions. Colors on the model plate interface reflect the assumed distribution of frictional parameters, (a – b)σ : red and blue colors denote stable and eff unstable regions, respectively. The shallow gray segment denotes the steady sliding zone. b Depth profiles of frictional parameters used to simulate an earthquake cycle scenario asperity of the earthquake cluster is continuously loaded distributions on a 3D model would improve our under- because of weak interplate coupling, probably due to the standing of earthquake cycles along plate boundaries. shrinkage of the strongly coupled region. In fact, the slip deficit rate in this region is less than that in the source Conclusions region of the off-Mie earthquake (Yokota et al. 2016). Using a local velocity model, we found that the 2016 Monitoring such seismic activities based on precise hy- Mw5.9 earthquake off the southeastern coast of Mie Pre- pocenters is important for the understanding of the state fecture occurred on the plate boundary along the Nankai of interplate coupling along the plate boundary fault, as Trough. We performed a numerical simulation based on well as geodetic data inversion. To achieve this, a local a hierarchical asperity model, which successfully repro- 3D velocity structure is crucial for accurate hypocenter duced the occurrence of a moderate earthquake in the locations (Akuhara et al. 2013; Yamamoto et al. 2017). middle of a megathrust earthquake cycle. The results of If other asperities exist, their ruptures could occur in the numerical simulation imply that the occurrence of the middle of the megathrust earthquake cycle due to the off-Mie earthquake represents shrinkage of the inter- the decrease in the interplate coupling ratio (Hori and seismically coupled area along the plate boundary and Miyazaki 2011). The timing and magnitude of the earth- serves as an indication of the preparatory processes of quake depends on both the location and size of the as- the next megathrust earthquake. perities. Our simulations can also model the magnitude, timing, frequency of the moderate earthquake, and the Methods/experimental interval of the fundamental rupture mode causing the Hypocenter determinations using a 2.5-D velocity model megathrust earthquake by defining the fault frictional The strong heterogeneity in the seismic velocity struc- parameters appropriately. In our 2D simulations, the ture in the subduction zone prevents precise hypocenter moderate M6 earthquake occurs 120 years after the M8 determinations. We used the 2-D velocity structure ob- earthquake, with a cycle of about 200 years. Both values tained from a tomographic inversion of the first-arrival are much longer than the observed values: the cycle of refraction travel times of a wide-angle OBS survey along megathrust earthquakes in the Nankai Trough is consid- line KI03, which passes through the hypocenter region ered to be 100 to 150 years (Ando 1975; Ishibashi 2004), (Fig. 2) (Wallace et al. 2016). Using this 2D structure, we and the off-Mie earthquake occurred 70 years after the created a 3-D velocity volume, assuming that the struc- previous 1944 Tonankai earthquake. In 3D simulations, ture is identical along the strike of the subducting plate; however, these periods would become shorter with the then, we computed theoretical travel times at each same frictional parameter values because the shrinkage DONET station. Since the PSP strikes about 40° N of the coupled area on the asperity also occurs in the around the hypocenter region (Citak et al. 2012) strike direction of the fault. Studies based on more real- (Additional file 1: Figure S1), the 2D structure is ori- istic fault geometry and frictional parameter ented in this direction (Additional file 3: Figure S3). We Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 12 of 17 ab Fig. 7 Simulation result. a Spatiotemporal evolution of the slip rate V during the resultant earthquake cycles. Colors denote the slip velocity. Black region corresponds to strong coupling on the plate interface. Arrows denote occurrences of earthquakes with denoted magnitude. The vertical color bar in the middle panel denotes the depth distribution of (a – b)σ . b Depth distributions of postseismic slip and slip rate just after the eff occurrence of a M6 earthquake indicated by solid arrows in a. Cumulative slip and slip rate profiles are shown in the left and right side panels, respectively. Colors of lines indicate days after the M6 earthquake. The broken gray line represents the coseismic slip distribution of the M6 earthquake call this the 2.5-D model for convenience because of the time at each station by averaging the residuals between similarity to the computations of wave propagation in the theoretical and observed travel times and redeter- three dimensions using a 2-D velocity model (e.g., Take- mined the hypocenter locations. This procedure was re- naka and Kennet 1996). We did not use the model of peated three times to obtain the final solution. The Kamei et al. (2012) because it did not fully cover the resulting site corrections are within ± 0.09 s, except at study area. KME19, which needed a correction of − 0.17 s. The source For the hypocenter determinations, we used the seafloor location errors are 0.7 km in the horizontal direction and seismic recordings from a broadband seismometer in- 0.6 km in the vertical direction for the mainshock and stalled at each DONET station. We manually picked 1.9 km in both directions for the aftershocks. The P-wave onsets from the vertical component seismograms. root-mean-square (RMS) of travel time residual is 0.02 s. To avoid the effect of possible heterogeneities in the ac- We examined the accuracy of the hypocenter cretionary prism in the strike direction of the PSP, we only depths by changing the dataset or assumed strike of used data from nine stations located within 30 km of the the PSP (see the Additional file 4: Supplementary mainshock source (white triangles in Fig. 2a). Moreover, note). Aftershock depths were several kilometers we did not use S-wave arrival times because of the uncer- shallower when we included distant stations tainty in S-wave velocity, especially in the soft sediments (Additional file 5:FigureS4).Thisdifference couldbe covering the accretionary prism. We used NonLinLoc caused by horizontal heterogeneities of the structure software (Lomax et al. 2000) for the hypocenter determi- along strike (e.g., Park et al. 2010). The increasing nations. This software adopts the eikonal finite-difference RMS travel time residual, including distant stations, scheme of Podvin and Lecomte (1991) to compute implies that theoretical travel times are less accurate first-arrival times at each station in a 3D volume for hypo- at distant stations because of inaccurate velocity center determinations using heterogeneous velocity struc- structures. Stable aftershock depths were obtained tures. We estimated site corrections for the P-wave travel when using S-wave arrival times or changing the Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 13 of 17 along the same line as the velocity structure, the profile in Fig. 5b was obtained along an intersecting line. Hence, the TWTs on this profile do not exactly match those computed from the hypocenter depths due to possible along strike heterogeneity of the subduction zone. We picked TWTs of strong reflections indicative of the plate boundary from the MCS reflection profile (Fig. 5a), converted them to depth values by inverting the process described above, and plotted them on the depth section (Fig. 3b). The uncertainties in picking the reflections are about 0.15 s, corresponding to 1.1 km at reflector depths. The errors in hypocenter depths correspond to a TWT error of 0.31 s for the mainshock and 0.67 s for the aftershocks. Two-dimensional forward simulation of M8-class earthquake cycles along the Nankai Trough Fig. 8 Temporal changes in shear stress on the plate interface. The Here, we attempt to interpret the conceptual relation same time interval as Fig. 7a at the four depths are shown. The blue between M8-class interplate earthquakes (such as the and red lines represent shear stress profiles at the initiation points of 1944 Tonankai earthquake) and the occurrence of a M6 (depth of 11 km) and M8 earthquakes (depth of 17 km), respectively. smaller earthquake (such as the 2016 off-Mie At depths deeper than the hypocenter of the M8 earthquake (black line), the stress accumulation rate is higher than that at the hypocenter earthquake) inside the source region through seismic because of secular loading from stable sliding occurring at a cycle modeling, which assumes the equilibrium of region deeper than the hyper-asperity. At depths shallower than quasi-dynamic shear stresses on the fault and the hyper-asperity (green lines), shear stress is almost constant frictional resistance governed by a rate- and because of stable sliding, except during the postseismic periods state-dependent friction (RSF) law. of M6 and M8 earthquakes. Dashed lines represent the reference stress level of each line, with corresponding colors Governing equations for quasi-dynamic seismic cycle strike of the PSP. Moreover, the mainshock depth was simulation stably located for most cases. Consequently, we con- In quasi-dynamic seismic cycle modeling, the occurrence of sider the obtained hypocenter locations to be stable. interplate earthquakes on plate boundaries has often been We also computed hypocenters by our routine analysis modeled as a perturbation from the steady-state relative based on a 1D velocity structure, using P- and S-wave plate motion by ignoring the contribution of steady (for- readings from all DONET stations (Fig. 3). The main- ward) slip on the whole plate interface. Thus, the interseis- shock hypocenter was located several kilometers east of mic process is expressed as the accumulation of slip deficit the 3D result, at a depth of 13.1 km. Aftershocks were (or backslip) in the seismogenic zone. The shear stress distributed several kilometers west of the 3D result, at a builds up with time as the slip deficit increases, and the ac- depth of about 15 km. We note that this computation is cumulated slip deficit is finally released as earthquakes different from that of Wallace et al. (2016), which only when the accumulated shear stress exceeds the fault resist- uses P-wave readings from stations connected to nodes ance governed by RSF. Below, we briefly describe equations D and E. governing the earthquake cycle simulation used here. At first, by discretizing the target plate interface Comparison of hypocenter depths with reflectors into N small sub-faults, we define physical variables obtained from MCS profiles associated with seismic cycles for each sub-fault: We compared the hypocenter depths with distinct re- shear stress τ,state variable θ representing contact i i flectors in the accretionary prism obtained from MCS state, and slip velocity V . Here, subscript i specifies reflection surveys (Fig. 5). We converted the hypocenter the particular sub-fault (1 ≤ i ≤ N). When the material depths to two-way travel times (TWTs) by integrating containing the fault (or a set of sub-faults) is a the slowness from sea level to the hypocenter depth linear-elastic one, temporal changes in shear stress picked at the epicenter, based upon the velocity struc- can be approximated through the whole seismic cycle ture used for the hypocenter determinations (Fig. 2b). using the quasi-dynamic expression given by Rice While the MCS profile shown in Fig. 5a was obtained (1993) as follows, Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 14 of 17 N rate of friction for a particular slip amount. Because dτ ðÞ t G ¼ − K V −V ðÞ t − V ðÞ t ; these parameters characterize the frictional property of ij pl j i dt 2V j¼1 the fault, and accordingly control the spatial and tem- poral behavior of fault slip, we must pay special atten- ð1Þ tion to the distribution of the frictional parameters on where K represents the static slip response function be- ij the model fault in seismic cycle simulations. tween sub-faults, and V is the relative plate velocity of pl The frictional stress on the plate interface is expressed the target plate interface, assumed to be constant re- as the product of the frictional coefficient μ and effective gardless of sub-faults. V and G represent shear-wave normal stress σ . Changes in fault frictional stress eff velocity and rigidity of the material containing the fault, (d(μσ )/dt) are assumed to balance changes in shear eff respectively. The first term on the right-hand side in Eq. stress (Eq. (1)). We then obtain simultaneous differential (1) denotes the contribution of slip deficit rate distribu- equations with respect to V and θ for 0 ≤ i ≤ N. Under i i tion to the change in shear stress. It should be noted the hypothetical distribution of frictional parameters a , that since only perturbation from the steady-state sub- b ,and L for 0 ≤ i ≤ N, we can simulate the spatiotempo- i i duction is treated, shear stress of sub-fault i is designed ral slip behavior on the target plate interface by integrat- so as not to be affected by sub-fault j if the j-th sub-fault ing the simultaneous differential equations from a is slipping with the same slip velocity as relative plate particular initial condition (and boundary condition if motion, V . Further, because the static slip response pl required). function K is used in the first term, the radiation damp- ij ing term, representing the effect of wave propagation, is Configuration of model plate interface and its discretization added as the second term to better approximate the In order to reproduce the complicated slip behavior dur- stress state through the entire seismic cycle, including ing an earthquake cycle, we need to set appropriate fric- the coseismic stage. tional parameters along the megathrust. To obtain a As previously mentioned, we assume that fault friction suitable set of parameters, we require many computa- governed by RSF resists fault motion during an earth- tions with trial parameter settings. To reduce the com- quake cycle. The friction coefficient μ following RSF is putational cost for the seismic cycle simulation generally defined as, including multiple earthquakes with a larger gap in mag- nitudes compared with the previous studies (Hori and V ðÞ t V θ ðÞ t i  i Miyazaki, 2010, 2011), we conducted 2D simulations in- μ ðÞ t ¼ μ þ a In þ b In i i V L stead of conventional 3D ones. In the 2D simulation, we ð2Þ modeled the plate boundary fault with a set of line seg- ments, as the vertical cross section perpendicular to the The first term μ is a reference friction coefficient at Nankai Trough axis off the Kii Peninsula. The configur- the steady-state sliding with the reference slip rate V , ation of the upper interface of the subducting PSP was while the second and third terms of the right-hand side simply assumed to be a smooth curve with an average on Eq. (2) represent contributions of slip rate (i.e., V) dip angle of about 12° (Fig. 6a). For this vertical cross and state variable (i.e., θ) to the friction coefficient, re- section, we discretized the curved plate interface into spectively. As the state variable θ is an indicator of the 4400 line segments of nearly equal lengths of about contact state at the fault interface, the third term repre- 50 m. We confirmed that this length of line segment is sents the dependency of the frictional coefficient on sufficient to accurately model earthquake cycles with the “contact” state. Several versions have been proposed for distribution of frictional parameters shown later. Each the evolution of the state variable θ to describe state line segment corresponds to a sub-fault; the slip re- variable-dependency on μ. We adopt a widely used ver- sponse function K for this type of sub-fault was calcu- ij sion, known as an “aging law”, given by, lated by the analytic expression given by Rani and Singh (1992), which determines the stress change under the dθ ðÞ t V ðÞ t θ ðÞ t i i i ¼ 1 − ð3Þ plane strain field due to the unit dip slip of each dt L sub-fault in an elastic half space consisting of a Poisson In this study, we simply set V = V . Parameters a, b, solid with G = 30 GPa and V = 3273 m/s. We assumed ∗ pl s and L in Eqs. (2) and (3) are called frictional parameters. V = 5 cm/year for the plate subduction velocity. pl As shown in Eq. (2), a and b are parameters concerning rate- and state-dependencies of friction, respectively. L Distribution of frictional parameters represents a characteristic slip distance necessary for the Colors on the model plate interface in Fig. 6a reflect the friction μ to decrease significantly. In other words, the assumed distribution of frictional parameters, (a − b)σ : eff lower the value of L assumed, the larger the decreasing Red and blue colors denote stable ((a − b)σ > 0) and eff Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 15 of 17 unstable ((a − b)σ < 0) regions, respectively. The first assumed that the shallowest interface (depth > eff blue-colored region roughly corresponds to the depth of 10 km) exhibits stable sliding with (a − b)σ > 0; how- eff the source region with a large slip in simulated Nankai ever, the stable sliding condition did not place a suffi- Trough megathrust earthquakes and is thus termed a cient load on the shallower asperity to trigger an M6 “hyper-asperity.” In the hyper-asperity, we identify two earthquake. Therefore, we added the steady sliding zone segments with L values smaller than those of the other at the top of the plate interface as a boundary condition region (dark blue). The shallower one, identified at a (gray segment in Fig. 6a, depth from 5 to 6 km), where depth of 10–11 km with a length along the plate bound- the plate interface was assumed to be slipping with a ary fault of 10 km (“shallower asperity” hereafter), causes constant slip velocity of plate subduction V . The appli- pl M6 earthquakes. The deeper one is defined at a depth of cation of this boundary condition is realized by exclud- 13–17 km, with a length of 30 km along the plate ing the upper-most sub-faults from the numerical boundary fault (“deeper asperity” hereafter). The deeper computation, because our model only considers the per- asperity is set to model the high slip deficit rate patch turbation from steady sliding on the plate interface. deduced from the inversion analysis of the slip deficit After many trial simulations with the above boundary rate distribution reported by Yokota et al. (2016) condition, we found that the frictional distribution (Additional file 2: Figure S2). The slip deficit rate during shown in Fig. 6b provides the scenario most consistent the interseismic period may reflect a degree of interplate with the above two constraints. In this study, we only coupling. If so, parts with high slip deficit rate may re- present the scenario obtained for the frictional parame- lease larger coseismic slip than parts with relatively low ters shown in Fig. 6b, which has implications for the slip deficit rate. Such heterogeneity in the slip deficit rate Nankai Trough megathrust earthquake cycles. (and possible resultant variation in the coseismic slip distribution) can be modeled by the heterogeneous dis- Additional files tribution of L in the seismic cycle simulation with RSF. Under the constant (a − b)σ condition, small L values Additional file 1: Figure S1. Names of DONET1 stations and nodes. eff Gray contours represent depth distribution of the top of oceanic crust generate a higher slip deficit rate and a larger coseismic layer 2 (Citak et al. 2012). (PDF 1118 kb) slip than large L values. Therefore, we assume a smaller Additional file 2: Figure S2. Comparison of hypocenter distribution L value at the deeper asperity than for the rest of the with the slip deficit rate (SDR). Red circles represent hypocenters hyper-asperity (light-blue region). Similarly, the shal- obtained in this study. SDR is after Yokota et al. (2016). (PDF 184 kb) lower asperity is also set to model the source of the Additional file 3: Figure S3. A schematic image of 2.5-D velocity model used for hypocenter determination. (PDF 165 kb) off-Mie earthquake with a smaller L value. Additional file 4: Supplementary note. (DOCX 35 kb) As proposed by Hori and Miyazaki (2010, 2011), it is Additional file 5: Figure S4. Dependence of hypocenter locations on known that the existence of a smaller asperity, defined (a) stations used, (b) use of S-wave first arrivals, and (c) strike direction of by a region with a much smaller L value relative to its PS plate. Values in brackets represent the RMS travel time residual. Plots surroundings, within the hyper-asperity tends to cause are along the line A–B shown in Fig. 2a. (PDF 82 kb) hierarchical rupture behavior. During the cycle of earth- quakes rupturing the entire hyper-asperity, smaller Abbreviations earthquakes rupturing only the smaller asperity can CMT: Centroid moment tensor; DONET: Dense Oceanfloor Network System for Earthquakes and Tsunamis; MCS: Multi-channel seismic survey; occur without propagating to the surroundings. Accord- OBS: Ocean-bottom seismometers; PSP: Philippine Sea plate; RMS: Root- ingly, both shallower and deeper asperities with smaller mean-square; RSF: Rate- and state-dependent friction; SSE: Slow slip event; L values in our model can produce the complex seismic TOC: Top of the oceanic crust layer 2; TWT: Two-way travel time; VLFE: Very- low-frequency earthquake cycle behavior of a hyper-asperity, depending on the as- sumed gap in L. Here, we simulate earthquake cycle sce- Acknowledgements narios primarily by modifying the gap in L values in the We appreciate comments from Dr. L. Wallace and an anonymous reviewer, hyper-asperity so that the following two behaviors are which greatly improved the manuscript. We also appreciate Dr. Y. Yamamoto for sharing data for plotting Fig. 2. The Japan Meteorological Agency unified satisfied: (i) M8-class earthquake that ruptures the earthquake catalogue was used to plot Fig. 1. All figures were drawn using hyper-asperity only occurs when the deeper asperity be- Generic Mapping Tools (Wessel and Smith 1998). gins to rupture and (ii) the shallower asperity only rup- tures a part of the hyper-asperity and its postseismic slip Authors’ contributions MN conducted the hypocenter determinations and proposed the model. MH (afterslip/slow slip) mainly occurs at the shallow part of performed the numerical simulations of the hierarchical asperity model. MN the shallower asperity. and MH wrote the paper. AN and MY processed the structural data from To satisfy the latter constraint (ii), the shallow side of MCS reflection surveys. TH and TT organized the geophysical interpretations. SKa and KS participated in the data processing for hypocenter determinations. the plate interface must gradually reduce the slip deficit SKo, NT, and YK organized the data acquisitions and the project. All authors during an M8 earthquake cycle so that shear stress on equally contributed to scientific discussions. All authors read and approved the the shallower asperity accumulates. To achieve this, we final manuscript. Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 16 of 17 Competing interests Hori T, Miyazaki S, Mitsui N (2009) A model of earthquake-generation cycle with The authors declare that they have no competing interests. scale-dependent frictional property—preliminary results and research plan for a project of evaluation for coming Tokai, Tonankai, and Nankai earthquakes. J Disaster Res 4:111–117. Publisher’sNote Hyodo M, Hori T (2013) Re-examination of possible great interplate earthquake Springer Nature remains neutral with regard to jurisdictional claims in published scenarios in the Nankai Trough, southwest Japan, based on recent findings maps and institutional affiliations. and numerical simulations. Tectonophys 600:175–186. https://doi.org/10. 1016/j.tecto.2013.02.038. Author details Ichinose GA, Thio HK, Somerville PG, Sato T, Ishii T (2003) Rupture process of the R&D Center for Earthquake and Tsunami, Japan Agency for Marine-Earth 1944 Tonankai earthquake (Ms 8.1) from the inversion of teleseismic and Science and Technology, Yokohama, Japan. Earthquake and Volcano regional seismograms. J Geophys Res 108:2497. https://doi.org/10.1029/ Research Unit, National Research Institute for Earth Science and Disaster 2003JB002393. Resilience, Tsukuba, Japan. Institute of Education, Research and Regional Ishibashi K (2004) Status of historical seismology in Japan. Ann Geophys 47: Cooperation for Crisis Management Shikoku, Kagawa University, Takamatsu, 339–368. Japan. Kamei R, Pratt RG, Tsuji T (2012) Waveform tomography imaging of a megasplay fault system in the seismogenic Nankai subduction zone. Earth Planet Sci Received: 28 December 2017 Accepted: 14 May 2018 Lett 317-318:343–353. https://doi.org/10.1016/j.epsl.2011.10.042. Kamei R, Pratt RG, Tsuji T (2013) On acoustic waveform tomography of wide- angle OBS data—strategies for pre-conditioning and inversion. Geophys J Int References 194:1250–1280. https://doi.org/10.1093/gji/ggt165. Akuhara T, Mochizuki K (2014) Application of cluster analysis based on waveform Kaneda Y, Kawaguchi K, Araki E, Matsumoto H, Nakamura T, Kamiya S, Ariyoshi K, cross-correlation coefficients to data recorded by ocean-bottom Hori T, Baba T, Takahashi N (2015) Development and application of an seismometers: results from off the Kii Peninsula. Earth Planets Space 66:80. advanced ocean floor network system for megathrust earthquakes and https://doi.org/10.1186/1880-5981-66-80. tsunamis. In: Favali P et al (eds) Seafloor observatories. Springer, Heidelberg, Akuhara T, Mochizuki K, Nakahigashi K, Yamada T, Shinohara M, Sakai S, pp 643–662. https://doi.org/10.1007/978-3-642-11374-1_252. Kanazawa T, Uehira K, Shimizu H (2013) Segmentation of the Vp/Vs ratio and Kaneko L, Ide S, Nakano M (2018) Slow earthquakes in the microseism frequency low-frequency earthquake distribution around the fault boundary of the band (0.1–1.0 Hz) off Kii Peninsula, Japan. Geophys Res Lett 45. https://doi. Tonankai and Nankai earthquakes. Geophys Res Lett 40:1306–1310. https:// org/10.1002/2017GL076773. doi.org/10.1002/grl.50223. Kawaguchi K, Kaneko S, Nishida T, Komine T (2015) Construction of the DONET Ando M (1975) Source mechanisms and tectonic significance of historical real-time seafloor observatory for earthquakes and tsunami monitoring. In: earthquakes along the Nankai Trough, Japan. Tectonophysics 27:119–140. Favali P et al (eds) Seafloor observatories. Springer, Heidelberg, pp 211–228. Annoura S, Hashimoto T, Kamaya N, Katsumata A (2017) Shallow episodic tremor https://doi.org/10.1007/978-3-642-11374-1_10. near the Nankai Trough axis off southeast Mie prefecture, Japan. Geophys Kisslinger C, Jones LM (1991) Properties of aftershock sequences in southern Res Lett 44:3564–3571. https://doi.org/10.1002/2017GL073006. California. J Geophys Res 96:11,947–11,958. Araki E, Saffer DM, Kopf AJ, Wallace LM, Kimura T, Machida Y, Ide S, Davis E, IODP Kodaira S, Hori T, Ito A, Miura S, Fujie G, Park J-O, Baba T, Sakaguchi H, Kaneda Y Expedition 365 shipboard scientists (2017) Recurring and triggered slow-slip (2006) A cause of rupture segmentation and synchronization in the Nankai events near the trench at the Nankai Trough subduction megathrust. trough revealed by seismic imaging and numerical simulation. J Geophys Science 356:1157–1160. https://doi.org/10.1126/science.aan3120. Res 111:B09301. https://doi.org/10.1029/2005JB004030. Ariyoshi K, Nakata R, MatsuzawaT,HinoR,HoriT, HasegawaA,Kaneda Y Lomax A, Virieux J, Volant P, Berge C (2000) Probabilistic earthquake location in (2014) The detectability of shallowslowearthquakes by theDense 3D and layered models: introduction of a Metropolis–Gibbs method and Oceanfloor Network system for Earthquakes and Tsunamis (DONET) in comparison with linear locations. In: Thurber CH, Rabinowitz N (eds) Tonankai district, Japan. Mar Geophys Res 35:295–310. https://doi.org/10. Advances in seismic event location, Kluwer, Amsterdam, The Netherlands, pp 1007/s11001-013-9192-6. 101–134. Bangs NLB, Moore GF, Gulick SPS, Pangborn EM, Tobin HJ, Kuramoto S, Taira A Matsuzawa T, Hirose H, Shibazaki B, Obara K (2010) Modeling short- and long- (2009) Broad, weak regions of the Nankai Megathrust and implications for term slow slip events in the seismic cycles of large subduction earthquakes. J shallow coseismic slip. Earth Planet Sci Lett 284:44–49. https://doi.org/10. Geophys Res 115:B12301. https://doi.org/10.1029/2010JB007566. 1016/j.epsl.2009.04.026. Mikumo T, Miyatake T (1979) Earthquake sequences on a frictional fault model Bird P (2003) An updated digital model of plate boundaries. Geochem Geophys with non-uniform strengths and relaxation times. Geophys J R Astr Soc 59: Geosyst 4:1027. https://doi.org/10.1029/2001GC000252. 497–522. Citak SO, Nakamura T, Nakanishi A, Yamamoto Y, Ohori M, Baba T, Kaneda Y Mochizuki K, Fujie G, Sato T, Kasahara J, Hino R, Shinohara M, Suyehiro K (1998) (2012) An updated model of three-dimensional seismic structure in the Heterogeneous crustal structure across a seismic block boundary along the source area of the Tokai–Tonankai–Nankai earthquake. In abstract of AOGS– Nankai Trough. Geophys Res Lett 25:2301–2304. AGU (WPGM) Joint Assembly, Singapore, 13–17 august 2012, abstract no. Mochizuki K, Nakahigashi K, Kuwano A, Yamada T, Shinohara M, Sakai S, OS06-A015. Kanazawa T, Uehira K, Shimizu H (2010) Seismic characteristics around the DeMets C, Gordon RG, Argus DF (2010) Geologically current plate motions. fault segment boundary of historical great earthquakes along the Nankai Geophys J Int 181:1–80. https://doi.org/10.1111/j.1365-246X.2009.04491.x. Trough revealed by repeated long-term OBS observations. Geophys Res Lett Furumura T, Hayakawa T, Nakamura M, Koketsu K, Baba T (2008) Development of 37:L09304. https://doi.org/10.1029/2010GL042935. long-period ground motions from the Nankai Trough, Japan, earthquakes: Moore GF, Bangs NL, Taira A, Kuramoto S, Pangborn E, Tobin HJ (2007) Three- observations and computer simulation of the 1944 Tonankai (M 8.1) and dimensional splay fault geometry and implications for tsunami generation. the 2004 SE Off-Kii Peninsula (M 7.4) earthquakes. Pure Appl Geophys 165: Science 318:1128–1131. https://doi.org/10.1126/science.1147195. 585–607. https://doi.org/10.1007/s00024-008-0318-8. Murotani S, Shimazaki K, Koketsu K (2015) Rupture process of the 1946 Nankai Hori T (2006) Mechanisms of separation of rupture area and variation in time earthquake estimated using seismic waveforms and geodetic data. J interval and size of great earthquakes along the Nankai Trough, southwest Geophys Res 120:5677–5692. https://doi.org/10.1002/2014JB011676. Japan. J Earth Sim 5:8–19. Nakanishi A, Kodaira S, Miura S, Ito A, Sato T, Park J-O, Kido Y, Kaneda Y (2008) Hori T, Miyazaki S (2010) Hierarchical asperity model for multiscale characteristic Detailed structural image around splay-fault branching in the Nankai earthquakes: a numerical study for the off Kamaishi earthquake sequence in subduction seismogenic zone: results from a high-density ocean bottom the NE Japan subduction zone. Geophys Res Lett 37:L10304. https://doi.org/ seismic survey. J Geophys Res 113:B03105. https://doi.org/10.1029/ 10.1029/2010GL042669. 2007JB004974. Hori T, Miyazaki S (2011) A possible mechanism of M9 earthquake generation cycles in the area of repeating M7∼8 earthquakes surrounded Nakanishi A, Shiobara H, Hino R, Kodaira S, Kanazawa T, Shimamura H (1998) by aseismic sliding. Earth Planets Space 63:773–777. https://doi.org/10. Detailed subduction structure across the eastern Nankai Trough obtained 5047/eps.2011.06.022. from ocean bottom seismographic profiles. J Geophys Res 103:27,151–27,168. Nakano et al. Progress in Earth and Planetary Science (2018) 5:30 Page 17 of 17 Nakano M, Hori T, Araki E, Kodaira S, Ide S (2018) Shallow very-low-frequency geology. Elsevier BV, pp 599–640. https://doi.org/10.1016/B978-0-444-62617- earthquakes accompany slow slip events in the Nankai subduction zone. Nat 2.00020-7. Commun 9:984. https://doi.org/10.1038/s41467-018-03431-5. Tsuji T, Ashi J, Strasser M, Kimura G (2015) Identification of the static Nakano M, Nakamura T, Kamiya S, Kaneda Y (2014) Seismic activity beneath the backstop and its influence on the evolution of the accretionary prism in Nankai trough revealed by DONET ocean-bottom observations. Mar Geophys the Nankai Trough. Earth Planet Sci Lett 431:15–25. https://doi.org/10. Res 35:271–284. https://doi.org/10.1007/s11001-013-9195-3. 1016/j.epsl.2015.09.011. Tsuji T, Kamei R, Pratt RG (2014) Pore pressure distribution of a mega-splay fault Nakano M, Nakamura T, Kamiya S, Ohori M, Kaneda Y (2013) Intensive seismic activity system in the Nankai Trough subduction zone: insight into up-dip extent of around the Nankai trough revealed by DONET ocean-floor seismic observations. the seismogenic zone. Earth Planet Sci Lett 396:165–178. https://doi.org/10. Earth Planets Space 65:5–15. https://doi.org/10.5047/eps.2012.05.013. 1016/j.epsl.2014.04.011. Nakano M, Nakamura T, Kaneda Y (2015) Hypocenters in the Nankai Trough Tsuji T, Minato S, Kamei R, Tsuru T, Kimura G (2017) 3D geometry of a plate determined by using data from both ocean-bottom and land seismic boundary fault related to the 2016 Off-Mie earthquake in the Nankai networks and a 3D velocity structure model: implications for seismotectonic subduction zone, Japan. Earth Planet Sci Lett 478:234–244. https://doi.org/10. activity. Bull Seism Soc Am 105:1594–1605. https://doi.org/10.1785/ 1016/j.epsl.2017.08.041. Tsuru T, Miura S, Park J-O, Ito A, Fujie G, Kaneda Y, No T, Katayama T, Kasahara J Noda H, Nakatani M, Hori T (2013) Large nucleation before large earthquakes is (2005) Variation of physical properties beneath a fault observed by a two- sometimes skipped due to cascade-up—implications from a rate and state ship seismic survey off Southwest Japan. J Geophys Res 110:B05405. https:// simulation of faults with hierarchical asperities. J Geophys Res 118:2924– doi.org/10.1029/2004JB003036. 2952. https://doi.org/10.1002/jgrb.50211. Utsu T (1961) A statistical study on the occurrence of aftershocks. Geophys Mag Noda H, Nakatani M, Hori T (2014) Coseismic visibility of a small fragile patch 30:521–605. involved in the rupture of a large patch—implications from fully dynamic Wallace LM, Araki E, Saffer D, Wang X, Roesner A, Kopf A, Nakanishi A, Power W, rate-state earthquake sequence simulations producing variable manners of Kobayashi R, Kinoshita C, Toczko S, Kimura T, Machida Y, Carr S (2016) Near- earthquake initiation. Progress Earth Planetary Sci 1:8. https://doi.org/10.1186/ field observations of an offshore M 6.0 earthquake from an integrated 2197-4284-1-8. w seafloor and subseafloor monitoring network at the Nankai Trough, Obana K, Kodaira S (2009) Low-frequency tremors associated with reverse faults southwest Japan. J Geophys Res 121:8338–8351. https://doi.org/10.1002/ in a shallow accretionary prism. Earth Planet Sci Lett 287:168–174. https://doi. 2016JB013417. org/10.1016/j.epsl.2009.08.005. Wessel P, Smith WHF (1998) New, improved version of generic mapping tools Obana K, Kodaira S, Kaneda Y (2004) Microseismicity around rupture area of the released. Eos 79:579. 1944 Tonankai earthquake from ocean bottom seismograph observations. Yamamoto Y, Takahashi T, Kaiho Y, Obana K, Nakanishi A, Kodaira S, Kaneda Y Earth Planet Sci Lett 222:561–572. https://doi.org/10.1016/j.epsl.2004.02.032. (2017) Seismic structure off the Kii Peninsula, Japan, deduced from passive- Obana K, Kodaira S, Kaneda Y (2005) Seismicity in the incoming/subducting and active-source seismographic data. Earth Planet Science Lett 461:163–175. Philippine Sea plate off the Kii Peninsula, central Nankai trough. J Geophys https://doi.org/10.1016/j.epsl.2017.01.003. Res 110:B11311. https://doi.org/10.1029/2004JB003487. Yokota Y, Ishikawa T, Watanabe S, Tashiro T, Asada A (2016) Seafloor geodetic Park J-O, Fujie G, Wijerathne L, Hori T, Kodaira S, Fukao Y, Moore GF, Bangs NL, constrains on interplate coupling of the Nankai Trough megathrust zone. Kuramoto S, Taira A (2010) A low-velocity zone with weak reflectivity along Nature 534:374–377. https://doi.org/10.1038/nature17632. the Nankai subduction zone. Geology 38:283–286. https://doi.org/10.1130/ G30205.1. Park J-O, Kodaira S (2012) Seismic reflection and bathymetric evidences for the Nankai earthquake rupture across a stable segment-boundary. Earth Planets Space 64:299–303. https://doi.org/10.5047/eps.2011.10.006. Park J-O, Moore GF, Tsuru T, Kodaira S, Kaneda Y (2003) A subducted oceanic ridge influencing the Nankai megathrust earthquake rupture. Earth Planet Sci Lett 217:77–84. https://doi.org/10.1016/S0012-821X(03)00553-3. Park J-O, Tsuru T, Kodaira S, Cummins PR, Kaneda Y (2002) Splay fault branching along the Nankai subduction zone. Science 297:1157–1160. https://doi.org/ 10.1126/science.1074111. Podvin P, Lecomte I (1991) Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools. Geophys J Int 105:271–284. Rani S, Singh SJ (1992) Static deformation of a uniform half-space due to a long dip-slip fault. Geophys J Int 109:469–476. Rice JR (1993) Spatio-temporal complexity of slip on a fault. J Geophys Res 98: 9885–9907. Sakaguchi A, Chester F, Curewitz D, Fabbri O, Goldsby D, Kimura G, Li C-F, Masaki Y, Screaton EJ, Tsutsumi A, Ujiie K, Yamaguchi A (2011) Seismic slip propagation to the updip end of plate boundary subduction interface faults: vitrinite reflectance geothermometry on Integrated Ocean Drilling Program NanTro SEIZE cores. Geology 39:395–398. https://doi.org/10.1130/G31642.1. Sakai S, Yamada T, Shinohara M, Hagiwara H, Kanazawa T, Obana K, Kodaira S, Kaneda Y (2005) Urgent aftershock observation of the 2004 off the Kii Peninsula earthquake using ocean bottom seismometers. Earth Planets Space 57:363–368. Seno T (2005) The September 5, 2004 off the Kii Peninsula earthquakes as a composition of bending and collision. Earth Planets Space 57:327–332. Takemura S, Shiomi K, Kimura T, Saito T (2016) Systematic difference between first-motion and waveform-inversion solutions for shallow offshore earthquakes due to a low-angle dipping slab. Earth Planets Space 68:149. https://doi.org/10.1186/s40623-016-0527-9. Takenaka H, Kennett BLN (1996) A 2.5-D time-domain elastodynamic equation for plane-wave incidence. Geophys J Int 125:F5–F9. Tobin H, Henry P, Vannucchi P, Screaton E (2014) Subduction zones: structure and deformation history. In: Stein R et al (eds) Developments in marine

Journal

Progress in Earth and Planetary ScienceSpringer Journals

Published: Jun 5, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off