An estimate of tidal and non-tidal modulations of plate subduction speed in the transition zone in the Tokai district

Abstract

Non-volcanic tremors and slow slip events in subduction zones have been found to be triggered by small external stress disturbances, as demonstrated by the synchronization of temporal variations in tremor rate with diurnal and semi-diurnal tides. Therefore, long-term variations in tremor rate might be predicted by amplitude modulations of diurnal and semi-diurnal tides at decadal time scales. Given that tremors and slow slips are shear slip on the plate boundary, their long-term variations must be associated with fluctuations in plate subduction speed below the seismogenic zone. In previous work, we showed a good correlation between long-term seismicity and empirically calculated tremor rate based on observed tidal levels in the Nankai region, western Japan. Here, we present an improved method of modeling long-term slip rate fluctuation based on the calculation of Coulomb stress due to ocean and solid earth tides on the plate interface. We also include the effects of non-tidal ocean variations, such as the Pacific Decadal Oscillation and the Kuroshio Current, employing an ocean model developed by the Japan Meteorological Agency. We apply the method to the Tokai district, where the effects of the Kuroshio Current are large, and demonstrate the importance of considering non-tidal effects. Our calculated slip rate fluctuations could amount to 1 mm/year in decadal scales, and periods with faster rates partly corresponded to variations in seismicity. Slow slip events in the study region weakly corresponded to times of higher stress.

Findings

Introduction

Recent seismological and geodetic observations have revealed diverse shear slip patterns on plate subduction interfaces (Ide et al. 2007a; Ide et al. 2007b), including non-volcanic tremors and slow slip events (SSEs) (Obara 2002; Shelly et al. 2006, 2007; Schwartz and Rokosky 2007; Beroza and Ide 2011; Vidale and Houston 2012). These events occur in the transition zone at depths of 30–40 km, where the frictional property on the plate interface changes from unstable to stable (Scholz 2002). In the seismogenic zone located just above the transition zone, the plate interface is locked partially, and stress accumulated due to plate motion releases with unstable slip. Below the transition zone, slip is always stable, and the oceanic plate is subducted at a constant velocity. In the transition zone, the presence of high-pressure fluids due to dehydration of the subducted slab reduces the effective normal stress (e.g., Kodaira et al. 2004; Peacock 2009). In such cases, a conditionally stable regime can exist, and even small stress disturbances can trigger tremors and SSEs, which are caused by remote earthquakes and diurnal and semi-diurnal tides (e.g., Miyazawa and Mori 2006; Nakata et al. 2008; Rubinstein et al. 2008; Lambert et al. 2009; Houston 2015).

Stress disturbances due to tremors and SSEs can trigger a major rupture on the plate boundary by providing a final push (Dragert et al. 2004; Beroza and Ide 2011), and interactions between these slow events and the major rupture zone have been actively studied using numerical simulations (e.g., Matsuzawa et al. 2010; Segall and Bradley 2012). To perform such simulations, the frictional properties on the plate boundary must be known. One advantage of observing the tidal responses of tremors is that a portion of the frictional properties in the transition zone at large spatial scales can be inferred, which is difficult to determine using laboratory experiments. The obtained properties can also be employed in numerical simulations of earthquake cycles by assimilating observational data (e.g., Hori et al. 2014), which allows for more precise long-term forecasts of large earthquakes in each anticipated region. In addition, tidal response measurement is useful for precisely monitoring small episodic stress changes (Nakata et al. 2008).

On the other hand, some studies have reported the presence of seasonality and possibility of decadal time-scale periodicity in the seismicity of large earthquakes in Japan (e.g., Ohtake and Nakahara 1999; Heki 2003; Tanaka 2014). These statistical behaviors indicate the presence of wide-area stress disturbances, and the effects of these disturbances on the seismogenic zone have been investigated. Currently, only Heki (2003) has successfully explained the seasonality in northern Japan quantitatively by using a snow-load model. The effects due to seasonal variations in atmospheric pressure are too small to trigger earthquakes (Ohtake and Nakahara 1999). Heki (2004) mentioned that ocean loading models have uncertainties due to steric (i.e., irrelevant with mass changes) components and that in situ observations of ocean bottom pressure (OBP) are necessary for its correct estimation. Kataoka (2007) estimated seasonal variations in OBP by using the Estimating the Circulation and Climate of the Ocean (ECCO; http://www.ecco-group.org) model and concluded that the amplitudes are insufficient for triggering earthquakes. Currently, such external stress disturbances due to climatic effects have not been incorporated in the above numerical simulations of earthquake cycles.

However, in the abovementioned studies, which attempted to explain the periodicities in seismicity, the high sensitivity of the transition zone to external stress disturbances was neglected. When considering their effects on the transition zone, it is possible that smaller stress disturbances trigger tremors and SSEs, which trigger earthquakes. However, the long-term variations in the occurrence of tremors and SSEs were not studied until recently. Pollitz et al. (2013) applied a hydrologic-loading model to explain the periodicity of slow events and the departures from periodicity at Vancouver Island, which lies in the Cascadia subduction zone. Ide and Tanaka (2014) (abbreviated as IT2014 hereafter) reported the possibility that tides could trigger tremors at the seasonal scale by applying an empirical relationship between observed tremor rates and tide-gauge records at Seto Inland Sea. Furthermore, IT2014 estimated a multi-decadal variation in tremor rate by focusing on long-term modulations of tidal amplitudes (Section 4.2.2 in Pugh 1987). The predicted tremor rate was dominated by an 18.61-year periodic variation, which was generated by modulations of diurnal and semi-diurnal tides. Long-term changes were well correlated with background seismicity during the period spanning from 1965 to 2010 in the Nankai Trough. Historical large earthquakes in this region also tended to occur when the predicted tremor activity with the 18.61-year period was high (IT2014).

However, in IT2014, tidal and non-tidal effects were not separated, and steric components were included because the model was based on observational tide-level data. In this study, the effects of non-tidal variations in OBP on the slip rate in the transition zone were calculated for the first time. For this estimate, a slip response model was applied to the Tokai district as an example. Fluctuations of the slip rates were demonstrated, and the importance of including non-tidal effects in comparisons with variations in seismicity was noted.

Methods

Method for computing slip rate

IT2014 presented the idea that tremors triggered by tides promote plate subduction in the tremor zone (Ide et al. 2007b), which activates seismicity of shallower earthquakes. In contrast to Beeler and Lockner (2003), who reported that diurnal tides cannot trigger earthquakes based on data from laboratory experiments, tremors were found to show tidal sensitivities at such high frequencies in the Nankai Trough (Yabe et al. 2014), which can indirectly trigger earthquakes. In the current study, this idea is generalized under the assumption that plate subduction speed in the transition zone, including in slow slip areas, changes due to external forces. However, to prove that fluctuations in subduction speed below the seismogenic zone can trigger earthquakes, quantitative simulations, such as those made by Hirose and Maeda (2013), would be necessary to calculate the 3-D stress distribution on the plate interface. The stress behaviors would strongly depend on given frictional parameters, and enormous efforts are needed to reveal the parameter sensitivity. In this paper, as a first step in producing such simulations, the influences of external forces, which serve as inputs, are estimated using the empirical model from IT2014, with the number of parameters being kept small.

The Tokai district is located on the east end of the Nankai Trough, where the Philippine Sea Plate subducts beneath the Eurasian (Amurian) Plate (Fig. 1a), which is magnified in Fig. 1b. The anticipated source area of a megathrust earthquake is surrounded by the thick black curve. The black dots and the green circle denote the epicenters of tremors (Idehara et al. 2014) and the area where a large slip occurred due to a long-term SSE from 2000 to 2006, respectively.

Fig. 1

The location of the Tokai district where the slip rate was estimated. a The blue curves denote the depth contours of the Philippine Sea Plates for 20, 30, and 40 km. The green stars denote the epicenters of historical large earthquakes with M > 7.5 (Set C in Tanaka 2014). The light blue circle represents a stationary cold-water mass, which appears when large-scale meandering of the Kuroshio Current occurs. b The thick black curve denotes the anticipated source area of a megathrust earthquake. The black dots and the green circle denote the epicenters of the tremors (Idehara et al. 2014) and the area where a large slip occurred due to a long-term slow slip event from 2000 to 2006, respectively. The solid earth and ocean tides along the red vertical line were compared in Fig. 2a. The slip rate at site H1 was mainly estimated. The ocean bottom pressures at H2–6 are shown in Fig. 3b. c, d The epicenters and magnitudes of small earthquakes with M > 3.5 in the Tokai area, for which the background seismicity from Fig. 4a was inferred. The events took place in an area within 150 km from the Suruga Trough in the direction of subduction (surrounded by the red diamond in (c))

Although IT2014 considered only the response of the tremor zone, our model considered the following three regions on the plate interface: I: seismogenic zone with A − B < 0, including the source area of a megathrust event; II: shallower (but deeper than the seismogenic zone) regions with A − B > 0, including the slow slip area; and III: deeper regions with A − B > 0, including the tremor zone. As in IT2014, slip responses can be expressed as

$$ {V}_i={V}_i^0 \exp \left(\varDelta {\mathrm{CFS}}_i/\left[{A}_i-{B}_i\right]\right) $$

(1)

where i = 1–3 denotes the index of the three regions, and V and V ⁰ denote the slip rate and the average slip rate of each region when external stresses are zero, respectively. ∆CFS represents the Coulomb failure stress (Scholz 2002) on the plate interface (∆CFS > 0 promotes slip), which results from solid earth and ocean tides (Tsuruoka et al. 1995) and from non-tidal variations in OBP.

The slip behavior of Eq. (1) takes the same form as in the steady state (e.g., Nagata et al. 2012); the validity of neglecting temporal variations in the state variable depends on the characteristic distance associated with the decay of the state variable, which has not been constrained well in geologically relevant settings and scales. The estimated slip rates depend on the assumed A–B and V ⁰, which can be high or low. Nevertheless, comparing results including and excluding the non-tidal effects reveals which is more dominant in a relative sense when variations in the state variable can be regarded as sufficiently slow.

Method for computing stress

Stress tensors on the plate interface were computed using Green’s functions from Okubo and Tsuji (2001). PREM (Dziewonski and Anderson 1981) was employed for density and elastic structure. The tidal potential of Tamura (1987) and the SPOTL software (Agnew 2012), which incorporates the regional ocean tide model of GOTIC2 (Matsumoto et al. 2000), were used to calculate the hourly solid earth and ocean tides, respectively. Non-tidal variations in OBP were estimated using the data assimilation model of the Meteorological Research Institute of Japan (Usui et al. 2006). The model input data included daily sea level atmospheric data, such as wind stress, rainfall, and thermal fluxes of JRA-25/JCDAS (Onogi et al. 2007). Temperature and salinity profiles and sea surface height (SSH) based on in situ and satellite observations are assimilated. A reanalysis of OBP from 1980 on grid points of 0.1 × 0.1 is available. The spatial resolution of the model is set sufficiently high to reproduce the variations in the Kuroshio Current. Thus, it is expected that a more realistic OBP is obtained than by using the ECCO, which has lower spatial resolution. The same Green’s function as for the ocean tides was used to calculate stress tensors.

∆CFS was calculated from the obtained stress tensors while considering the geometry of the plate boundary (Baba et al. 2002; Nakajima and Hasegawa 2007; Hirose et al. 2008). The frictional coefficient μ in ∆CFS was set to 0.2 throughout the current study. A similar value was obtained in Cascadia (Houston 2015). Hereafter, the ∆CFS caused by solid earth and ocean tides, the sum of these two ∆CFS _S and the non-tidal effects are abbreviated as ∆CFS _S , ∆CFS _O, ∆CFS _S+O , and ∆CFS _NT , respectively.

Loading effects on asperities

The temporal variations in the stress in the seismogenic zone (region I) can be represented as

$$ \frac{d}{dt}{\tau}_1={\displaystyle {\sum}_{i=1}^3{K}_{1i}\left({V}_{pl}-{V}_i\right)+\frac{d}{dt}\varDelta {\mathrm{CFS}}_1} $$

(2)

where K _1i denotes the elastic Green’s function (e.g., Hirose and Maeda 2013). Of the four terms on the right-hand side of Eq. (2), the second term is the most important. ${V}_1^0$ appears in the first term after applying Eq. (1). ${V}_1^0$ represents the slip rate in the seismogenic zone and is much smaller than ${V}_2^0$ and ${V}_3^0$ (${V}_1^0=0$ if fully locked). Moreover, the shear stress in the seismogenic zone generated by unit slip on region II is larger than that from region III (K ₁₂ ≪ K ₁₃ < 0) because the regions I and II are contiguous. In addition, the amplitudes of ∆CFS₂ near the coast are larger than ∆CFS₃, which is located below land. For example, Fig. 2a shows the amplitudes of ∆CFS _S and ∆CFS _O on Jan. 1, 1980, along the N-S line from Fig. 1b. The effects of ocean tides were observed to decay rapidly as the distance from the coast increased, whereas the effects of solid earth tides were almost constant. ∆CFS _O and ΔCFS _S+O almost agree in phase, implying that the sum of the effects of solid earth and ocean tides also decreases toward the north.

Fig. 2

Changes in stress and slip rate estimated at H1 on the plate interface caused by tides and non-tidal variations in ocean bottom pressure (OBP). a A comparison in peak-to-peak amplitudes in ΔCFS between solid earth and ocean tides on Jan. 1, 1980, along the red line in Fig. 1b. b, c The red and blue lines denote ΔCFS caused by ocean tides and that caused by the sum of solid earth and ocean tides, respectively. d, e Hourly, daily, and monthly slip rates calculated using Eq. (1) with ΔCFS being caused by the sum of solid earth and ocean tides shown in (c). f ΔCFS caused by the non-tidal variation in OBP in Fig. 3b. The spatial distribution of OBP near the coast was considered in computing ΔCFS at H1. g, h Hourly, daily, and monthly slip rates calculated using Eq. (1), with ΔCFS being caused by the sum of solid earth and ocean tides and the non-tidal (NT) variation in (f). i, j The same as (f, h), but the horizontal axis denotes months, and the results from 1980 to 2013 are cyclically superimposed

Model parameters

At the site marked by H1, ΔCFS was calculated using dip, strike, and rake angles of 17°, 250°, and 140°, respectively. The slip rate was calculated using Eq. (1) for ${V}_2^0=4\;\mathrm{cm}/\mathrm{year}$. By observing short-term tidal responses, Yabe et al. (2014) found that A − B in Eq. (1) was 2–3 kPa in almost all of the tremor zones in the Nankai Trough. By contrast, A − B in the slow slip area was set to −17 kPa by Hirose and Maeda (2013) to simulate SSE recurrence and the main ruptures in the Tokai area, which were based on a rate- and state-dependent law. However, the chosen A − B magnitude is not unique because slip behaviors also depend on the characteristic distances discussed above. In this study, the magnitude was set to a similar value as in Yabe et al. (2014) (A − B = 2 kPa) because the tremor zone and the slow slip area were adjacent to each other. The case of A − B in the slow slip area being negative (conditionally stable) will be discussed later.

The slip rates were converted to rates of stress change at the lower end of region I, assuming K ₁₂ = −8 MPa/m in Eq. (2). The value of K ₁₂ was calculated by applying the method of Okada (1992) to a 2-D plate interface with a dip angle of 17°, assuming that K ₁₂ represented the response at H1 (depth 25 km) to a uniform slip on region II, occupying depths from 25 to 35 km. The strike component of the slip was neglected, and isotropic elasticity was assumed.

Earthquakes and SSEs in the Tokai district

The computed slip and stress rates were compared in the “Results and discussion” section using the following data. Figure 1c, d shows the seismicity of shallower earthquakes with M > 3.5. In this area, most events occurred off of the plate interface. The principal cause of these events is considered to be stress accumulation due to plate motion, and it is expected that a faster slip rate enhances inter-plate coupling and corresponds to higher seismicity. The long-term background seismicity corresponding to these events was obtained by the declustering method from Ide (2013).

Mogi (1969) noted the seasonality of great earthquakes in the Nankai Trough. In the Tokai district, earthquakes occurred on Feb. 3, 1605, Oct. 28, 1707, and Dec. 23, 1854. The occurrence months of these events were compared with seasonal variations in the slip rates. The background seismicity of the shallower earthquakes did not show seasonal variations at a statistically significant level.

Long-term SSEs occurred during the early 1980s, late 1980s, the first half of the 2000s (e.g., Kobayashi 2014) and from Nov. 2013 to the time of this writing (Figs. 65 and 66 in Geographical Information Authority of Japan 2015). Recurrence intervals of ~10 years can be reproduced by numerical simulation (Hirose and Maeda 2013). However, why these slips occurred during these periods has not been explained. Slip and stress rates were also compared with the occurrences of these events.

Results and discussion

The effects of tidal modulations on ∆CFS and slip rate

Figure 2b shows the hourly ∆CFS caused by tides between 1980 and 1981 at H1. At the inland H1, ∆CFS _O (red) was an order of magnitude smaller than ΔCFS _S+O (blue). The half amplitudes of ΔCFS _S+O were 0.4–1.2 kPa, which were almost the same as those estimated for the Nankai Trough by Nakata et al. (2008). In the time-series plot of ΔCFS _S+O , fortnightly and semi-annual modulations of the amplitudes are visible. The fortnightly modulations result from neap and spring tides. Figure 2c shows the ΔCFS _O and ΔCFS _S+O between 1980 and 2013 at H1. At this horizontal scale, the semi-annual, 4.425- and 18.61-year modulations were dominant. The variation in amplitudes of ΔCFS _S+O over an 18.61-year cycle was between 10 and 20 %.

Figure 2d, e shows the slip rates obtained from the above value of ΔCFS _S+O . In the hourly slip rate of (d), the non-linear amplification was remarkable, whereas the decadal variations in the slower slip rates near 2.5 cm/year were smaller (i.e., the minimum slip rates were flatter) than those in the larger slip rates of approximately 7 cm/year. The peak-to-peak differences in the maximum slip rates were approximately 5 and 10 mm/year for the 4.425- and 18.61-year periods, respectively, corresponding to 13 and 25 % of the assumed reference rate of 4 cm/year. Comparing the hourly rate in (d) with the daily and monthly rates in (e), longer-period variations were more emphasized for the slip rates that were averaged over longer durations. In the monthly rate of (e), the difference in the amplitude over 18.61 years was 0.3 mm/year. In addition, the average slip rate in the monthly rate of (e) was approximately 0.9 mm/year faster than the assumed reference rate. In contrast with the case of ordinary tidal deformations, the non-linear amplification due to the exponential in Eq. (1) produces bias.

The effects of non-tidal ocean loads on ∆CFS and slip rate

To confirm the validity of the ocean model, Fig. 3a shows the monthly running average of the calculated OBP, which is superimposed on Fig. 3 from Matsumoto et al. (2006). Overall, our model, which is represented by the blue curve, reproduces the inter-annual tendencies of the observed OBP, except for the lower peaks in the model at approximately 2002.8 and 2002.9.

Fig. 3

Non-tidal variations in ocean bottom pressure (OBP). a A comparison under equivalent water thickness to OBP. The blue line denotes the OBP estimated by the ocean model used in this study, which was developed by the Japan Meteorological Agency (Usui et al. 2006). The OBP was estimated at a longitude of 146°E and a latitude of 39.2°N. The red and black lines represent in situ observations from the OBPR-B site and another ocean model, respectively (Fig. 3 in Matsumoto et al. 2006). For easier comparison, the other plots were removed from the original figure in Matsumoto et al. (2006). b The daily OBP calculated by the JMA ocean model at the sites marked H2–6 in Fig. 1b. The orange curve on H2 denotes a 5-year running average. The superimposed pink bars denote periods of large-scale meandering of the Kuroshio Current, as reported by the JMA. It can be observed that the OBPs along the coast in the Tokai district increase (decrease) with the initiation (end) of the events. c The same as (b), but the OBP values from 1980 to 2013 are plotted cyclically against months

Figure 3b shows the daily OBP from 1980 to 2013 corresponding to H2–6 from Fig. 1b. On the whole, the OBPs at these sites were similar, except for their amplitudes. In a closer examination of the time-series plots, the variation patterns nearer the coast seemed to slightly differ from those farther from the coast. At H2, 4, and 6, the amplitudes were relatively large, and the increases in OBP from 1988 to 1992 and 2000 to 2005 were clearer. At H3 and 5, which had smaller amplitudes, the increase corresponding to the latter period began in 2002.

Decadal variations in the SSH along the Japanese coast at mid-latitudes mainly result from westward propagations of Rossby waves, which are excited by changes in wind stress curl over the central North Pacific (e.g., Qiu and Chen 2010). Satellite altimeter data and the results from a numerical model together indicated that the SSH at a longitude of 140°E and latitude of 32°–34°N increased from 1989 to 1994 and 2001 to 2005 (Fig. 6b in Qiu and Chen 2010). Figure 3b shows that these high SSH anomalies correspond to high OBPs along the coast, which is more clearly visible in the 5-year running average of the OBP at H2.

With respect to OBP, shorter-term variations due to the large-scale meandering of the Kuroshio Current were also observed. These periods are marked by pink bars (http://www.data.jma.go.jp/kaiyou/data/shindan/b_2/kuroshio_stream/kuroshio_stream.html). The meandering began with the appearance of a stationary cold-water mass off the Kumano-nada (Okada 1978) (Fig. 1a). Figure 3b shows that OBP increases when such meandering occurs.

Figure 2f shows the computed daily ΔCFS _NT at H1. When OBP was positive, ΔCFS _NT was negative. The variation pattern in ΔCFS _NT was similar to that at H2, 4, and 6, reflecting the ocean loads nearer the coast more strongly. The amplitudes of ΔCFS _NT in the decadal scales were 50–100 Pa. It was also observed that ΔCFS _NT rapidly increased after the meandering events of the Kuroshio Current. A meandering event that started in 1935 and lasted for approximately 10 years was previously reported (Okada 1978). The 1944 Tonankai and 1946 Nankaido earthquakes occurred shortly after the end of this meandering episode, which was consistent with the increase in ΔCFS _NT .

Figure 2g, h shows the slip rates obtained for the sum of ΔCFS _S+O and ΔCFS _NT . The 4.425-year modulation was still visible in the hourly slip rates because ΔCFS _NT was relatively small (Fig. 2g). In the long-term averages of the slip rates (Fig. 2h), decadal variations in non-tidal effects were emphasized more, and their variation patterns were similar to those of the OBP, as shown in Fig. 3b. Moreover, when comparing the monthly slip rate in Fig. 2h with the rate that was caused by only the tides (Fig. 2e), the fluctuation of the former was approximately twice as large as the latter. The 18.61-year change in ΔCFS _S+O mainly consists of modulations of diurnal and semi-diurnal oscillations, and the long-term bias is almost constant with time. Thus, the contributions of tides to long-term slip rates are effectively smaller, despite their larger amplitudes over shorter periods (1–2 kPa). In contrast, slow changes in ΔCFS _NT increase/decrease the decadal-scale slip rate effectively, even though the added non-tidal effects were only about 1/20th of the tidal ones (~0.05–0.1 kPa).

Incidentally, the fluctuation of the slip rate in region III, where the ocean effects are much smaller than in region II, is mostly caused by solid earth tides (Fig. 2a). Consequently, its long-term variation pattern becomes similar to Fig. 2d, e, except that the amplitudes are slightly decreased.

A comparison of slip rate with seismicity of shallower earthquakes

Figure 4a shows the seismicity of the shallower earthquakes that occurred between 1982 and 2010. The annual averages of the stress rates obtained from the slip rates in Fig. 2d,g are superimposed.

Fig. 4

A comparison between slip rates and seismicity. a The red solid line denotes the annual average of the slip rate in Fig. 2g, where the tidal and non-tidal effects are considered. The result is given in units of the rate of stress change by multiplying the slip rates by an elastic constant, based on Eq. (2). A 5-year running average is also superimposed to emphasize the decadal variations. The blue curve denotes the annual average of the rate of stress change, computed from the slip rate shown in Fig. 2d, where only the tidal effects are considered. The yellow circles denote the background seismicity, inferred by the method of Ide (2013). The employed earthquakes are shown in Fig. 1c,d. The superimposed pink bars denote periods of large-scale meandering of the Kuroshio Current, as in Fig. 3b, during which the seismicity seemed to decrease. When excluding non-tidal effects, the increase in seismicity during the mid-1990s cannot be explained. b The horizontal and vertical axes denote months and slip rate, respectively. The red and blue symbols denote the monthly averages in Fig. 2h,e, respectively, from 1980 to 2013, which are plotted cyclically against months. The green line represents the frequency of great historical earthquakes in the Tokai area (1605, 1707, and 1854), for which Mogi (1969) noted the seasonality. c The red and blue curves denote the stress changes obtained by integrating the rates in Fig. 2g with respect to time. The purple bars show the periods when four long-term slow slip events occurred in the Tokai area (Kobayashi 2014; Geographical Information Authority of Japan 2015). d The red and blue curves denote the slip rates for cases including and excluding non-tidal effects, respectively, when A − B is assumed to be negative (see text). The purple bars show the periods of the SSEs

When only the tidal effects were considered, the stress rate was maximized between approximately 1988 and 2006, whereas three upper peaks appeared during 1985, 1997, and 2007 when the non-tidal effects were included. When including the non-tidal effects, the stress rate in 1997 was lower than the other two peaks, occurring approximately in 1985 and 2007, which could not explain the high seismicity observed from 1993 to 1999. However, the decreased seismicity in approximately 1990 and 2005 agreed well with the lower stress rate peaks in these years, which were caused by meanderings of the Kuroshio Current.

The overall agreement in phase between stress rate and seismicity is a reminder of the fact that earthquakes are more likely to be triggered near a phase with a maximum stress rate, when the periodic stress disturbance is relatively small compared to the secular stress accumulation rate (Heki 2003). The epicenters of the inland events in Fig. 1 were located away from the plate interface, where the stress change caused by the fluctuating slip rate on the plate interface might be attenuated.

Seasonality of large earthquakes

Figure 3c is a plot of the OBP during the period from 1980 to 2013, which is cyclic with respect to months. The time-series plots at H2–6 indicate that the OBP tends to be smaller from August to November. Figure 2c shows the corresponding daily ΔCFS _NT , which peaks around September. Figure 4b shows monthly slip rates for the cases including and excluding the non-tidal effects. When considering only the tidal effects, the slip rates are highest in January and July. When adding the non-tidal effects, the slip rates are larger around September.

The decreasing OBP and increasing slip rate in September caused landward and seaward displacements, respectively, at the coastal GNSS stations. Both of these effects on displacements are expected to be less than 1 mm/year and to cancel each other out, which would be difficult to distinguish given the precision of the GNSS data.

The occurrence months of the great earthquakes were closer to the periods when stress was maximized (i.e., 3 months later than the phase of maximum slip rates), although the number of events is too small to find a clear relationship (Fig. 4b).

The above example of seasonal scale also suggested the importance of considering non-tidal effects when discussing the relationship between slip rates and seismicity.

Comparing slip rate with long-term SSEs

Compared to the stress rate in Fig. 4a, SSE initiation seems to be delayed by 3–4 years. According to Heki (2003), when a stress disturbance is larger, earthquakes are more likely to be triggered under maximum stress. SSEs occur on plate interfaces and therefore should be more strongly affected by slip rate fluctuations on the plate interface than shallower inland events. To confirm this, Fig. 4c shows the stress change obtained by integrating the hourly slip rate in Fig. 2g with respect to time. The stress became relatively large from 1986 to 1991, 1998 to 2002, and 2007 to 2013. Except for the first event, which possessed a smaller moment magnitude (Kobayashi 2014), the three events seemed to initiate 2–4 years after the computed stress reached a maximum. The cause behind the delayed onsets of the SSEs is not clear. More events are necessary to confirm whether a weak correspondence between stress and SSE initiation continues in the future.

Can SSEs be interpreted by pure responses to external forces?

In this section, a hypothetical case of region II having a negative A − B is considered. In this case, the slip rate represented by Eq. (1) for i = 2 could be interpreted as slow slip. To prevent the slow slip from developing into a rapid rupture, a frictional law, such as in Matsuzawa et al. (2010), is implicitly assumed.

For example, Fig. 4d shows slip rates for cases including/excluding non-tidal effects when A − B = −2 kPa. When considering only tidal effects, the long-term rates possessed a similar pattern to those in the case with a positive A − B (Fig. 4c) because the long-term slip rates became faster when the diurnal and semi-diurnal tidal amplitudes were larger. When adding non-tidal effects, the rate became larger from 1988 to 1991, 1999 to 2005, and 2010 to 2013. The SSEs in 1988, 2000, and 2013 occurred during or near periods with faster slip rates. For the earliest event, which occurred from 1980 to 1982, the correlation was unclear due to the lack of pre-1980 data.

However, the increase in slip rate from 2000 to 2005 was approximately 0.5 mm/year, which was more than an order of magnitude smaller than the observed slip rate for the 2000 event (Fig. 14 of Suito and Ozawa 2009). Thus, a direct response to external tidal and non-tidal stresses cannot explain the SSE when considering the magnitude of the frictional parameter that was determined by observing tremors (2 kPa). To interpret the SSE as a direct response of the slow slip area to external stresses using our model, A − B must be even smaller.

If the value of A − B is correct, then to interpret the SSE as a direct response, a mechanism is necessary in which the slip rate is amplified. One possibility may involve the effect of random noises on the system. Once the slip rate begins to increase, the effects of local stress disturbances caused by microfractures and episodic slow events to promote the slip are more easily increased by non-linear effects. If these disturbances occur randomly in time and their amplitudes are larger than the fluctuations caused by tidal and non-tidal effects, so-called stochastic resonance occurs, which would emphasize the effects of smaller disturbances with slower time scales. To confirm this, however, an earthquake cycle simulation study that incorporates external forces and noises is necessary.

A note on the periodicity of large earthquakes and the PDO

A possible periodicity in historical large earthquakes in Japan

Thus far, the importance of non-tidal modulations in estimating slip fluctuations has been demonstrated, and the purpose of this study has been attained. In this section, some decadal time-scale phenomena, which seem to be related to tidal and non-tidal modulations, are introduced to provide material for future studies to help elucidate decadal triggering mechanisms of large earthquakes.

Tanaka (2014) reported a possible historical periodicity for earthquakes with M > 7.5 in northeast Japan. The epicenters of 22 events are shown in Fig. 1a. The reasons for choosing these events and the risks of using incomplete historical data are described in detail by Tanaka (2014). The periodicities were investigated for 8.85 and 18.61 years associated with extremely weak tides. Figure 3b,e in Tanaka’s paper shows relatively strong and weak correlations for the 8.85 and 18.61 year periods, respectively. If the principle of superposition holds, then the effects of these two periods should be combined when constructing a histogram to test periodicity. However, to define periods with greater stress, a stress threshold above which earthquakes are more likely to be triggered must be set for combined disturbances consisting of larger and smaller peaks. Such a threshold would depend on the amplitude of each disturbance, which is unknown. To avoid such uncertainties, the two periods were treated separately.

Why do larger earthquakes show stronger correlations?

The hypothesis that fluctuations in subduction speed in the transition zone affect the seismicity of shallower earthquakes seems to hold for earthquakes of all magnitudes. However, smaller events tend to show weaker correlations with long-term tidal periods (Fig. 3c, g in Tanaka 2014).

The lower correlations shown by smaller earthquakes could be qualitatively explained if the characteristic distance discussed in the “Method for computing slip rate” (denoted by D _c) depends on the scale of irregularities, such that smaller (larger) irregularities have shorter (longer) characteristic distances (e.g., Matsu’ura et al. 1992). The characteristic time of the decay of the state variable can be approximated as D _c/V, where V is the slip rate. It follows that large-scale irregularities have longer characteristic times when the slip rate is the same. Therefore, the effects of slower stress disturbances caused by decadal-scale modulations of the state variable would persist for longer on larger asperities. Tidal and non-tidal variations have long wavelengths, which change the stress over an asperity, corresponding to a megathrust event.

Han et al. (2014) showed that tremors and slow slips can trigger small earthquakes. At small spatial scales, local stress changes caused by such slow events occur more frequently and are more dominant than long-term stress disturbances caused by tides and non-tidal effects, which may result in less clear correlations with decadal-scale modulations.

Can the PDO be related to the occurrence of large earthquakes?

Although there are insufficient evidences to answer this question, a working assumption could be made. When the PDO index is positive, negative SSH anomalies are created in the eastern North Pacific. These anomalies propagate west over 3–4 years and decrease the SSH along the Japanese coast at mid-latitudes. Recent negative SSH anomalies were observed in approximately 1998 and 2008 at a longitude of 140° (Fig. 6 in Qiu and Chen 2010). These periods corresponded to lower OBPs and faster rates of stress change in region I, when A − B in region II is positive (Fig. 4a). The 8.85-year peak in seismicity found in historical large earthquakes corresponds to T = 2011.7 − 8.85n (year), where n is an integer number (Fig. 3c in Tanaka 2014). The two recent periods (2011.7 and 2002.9) were delayed by 4–5 years from periods with lower OBPs. The stress changes reached maximums approximately 2.5 years after 1998 and 2008, when the period of disturbance was 10 years. If differences of 1–2 years are allowed, it may be possible that the peak of the seismicity corresponded to a period when stress was at a maximum.

Scafetta and Mazzarella (2015) reported that the PDO index and a world historical earthquake record had common spectral peaks at approximately 5.7, 8.9, and 20 years. Yasuda (2009) found that 18.61-year tides have been correlated with the historical PDO index since the eighteenth century. Therefore, the 20-year peak found by Scafetta and Mazzarella (2015) might be related to the 18.61-year tides. The mechanisms generating the other spectral peaks in the PDO index have not been revealed. The long-term stability of the phase of the periodic variation in seismicity along the Japan Trench could indicate that not only the 18.61-year but also the 8.85- and 4.425-year tidal modulations are connected with the generation of PDO. However, 8.85-year tidal modulations are weaker than 4.425-year modulations (e.g., Haigh et al. 2011). Thus, the fact that approximately 8.85-year variations are more remarkable than 4.425-year variations in the seismicity and the PDO index is mysterious.

The dominance of the 8.85-year variations in seismicity could be explained if the characteristic distance mentioned above depended on the scale of irregularities (i.e., longer disturbances survive longer). For the PDO index, it is expected that oceanographic studies will elucidate the origin of the approximately 9-year variation to confirm whether the periodicity of historical large events is real.

Conclusions

The long-term fluctuations in subduction speed due to tidal and non-tidal modulations in the Tokai slow slip area were estimated using an empirical slip response model, under the assumption that the effects of stress disturbances on the slow slip area are dominant in changing inter-plate coupling and that variations in the state variable are small. Non-tidal effects on the long-term slip fluctuations are larger than tidal effects in areas covered by the ocean, such as Tokai, and non-tidal contributions should not be neglected when discussing temporal coincidences between obtained slip rates and seismicity.

Absolute amplitudes of estimated slip rates depend on frictional parameters. Assuming that the effective normal stress is as low as in the tremor zone, long-term slip rate fluctuations could reach 1 mm/year, suggesting the possibility that the parameters could be constrained by crustal deformation data. However, the fluctuation amplitudes were more than an order of magnitude smaller than the slip rates for past SSEs, and it was not possible to interpret the slip rate fluctuation of an SSE unless the effective normal stress was even smaller.

To determine whether changes in shear stress resulting from fluctuations in slip rate can trigger slow events and earthquakes, our model, which was used to evaluate external stresses, is oversimplified. For example, it is necessary to consider interactions between the transition zone and locked zones based on a rate- and state-dependent frictional law (e.g., Matsuzawa et al. 2010; Shibazaki et al. 2011; Hirose and Maeda 2013; Hori et al. 2014).

To perform such simulation studies, it is important to include both non-tidal and tidal effects. In this study, OBP was evaluated with an ocean model to assimilate observed atmospheric and oceanographic data. In situ observations of OBP (e.g., Ito et al. 2013) have occurred in only recent history, such that wide-area annual variation patterns have not been adequately revealed. On the other hand, a regional-scale OBP was revealed by a satellite gravity mission (e.g., Petrick et al. 2014). Assimilating observations with such different space-time resolutions will be key toward improving the ocean model. In situ and satellite observations, climatology, and oceanography will be helpful in providing more accurate estimates of external stress disturbances for use as inputs in numerical simulations.