Slow slip event in the focal area of the 1975 Kurile tsunami earthquake inferred from unusual long-term seismic quiescence

In subduction zones, slow slip events (SSEs) have been observed in the portion deeper than the downdip edge of seismogenic zone. However, since it is far offshore from geodetic networks on land, shallow SSEs near the trench axis are hardly observed. Despite of less quantitative than seafloor geodetic observations, a method for inferring the shallow SSEs based on seismic quiescence was presented in this study. Unusual decrease in occurrence rate of M ≥ 5.0 earthquakes was found in the southwestern Kurile Islands. The occurrence rate was ∼1.3 events/year between 1977 and around 2004 and no earthquake was observed during 16 years after 2004. The spatial pattern of the seismic quiescence can be explained qualitatively by the Coulomb failure stress change due to shallow SSE and its fault plane is on the upper boundary of the subducting Pacific plate in the focal area of the 1975 tsunami earthquake.


Introduction
Two megathrust earthquakes occurred in 1969 and 1975 at almost the same place in the southwestern Kurile Islands subduction zone.The 1969 earthquake was a typical high-frequency earthquake in this region, whereas the 1975 earthquake was a tsunami earthquake that caused a tsunami which was much higher than predicted from the magnitude determined by the short-period seismic wave (Sapporo District Meteorological Observatory 1976;Abe 1989).The tsunami source area of the 1975 event is closer to the trench axis than that of the 1969 event (Hatori 1975) and the aftershock distribution according to the earthquake catalog of the Japan Meteorological Agency also supports the near-trench location of the 1975 event.Fukao (1979) presented a hypothesis that the shear stress was loaded on the shallower plate interface after the coseismic faulting of the 1969 event and the 1975 event was triggered trenchward due to this steep accumulation of the stress.Ioki and Tanioka (2016) conducted tsunami waveform inversion and found that the large slip area of the 1975 event is closer to the trench axis than that of the 1969 event.This result supports the Fukao's hypothesis.
A slow slip event (SSE) occurred in September and October 2014 at the Hikurangi subduction margin offshore New Zealand (Wallace et al. 2016).By using data from the absolute pressure gauge (APG), they found that the fault plane of SSEs was very close to the trench axis and extremely shallow (within 2 km from the seafloor).In the subduction zone many previous studies reported deep SSEs at the downdip limit of the seismogenic zone, i.e., at depths of 30-40 km (Schwartz and Rokosky 2007).Comparing with the deep SSEs, the Hikurangi SSEs are in the obviously shallow plate boundary with low temperature and low pressure.Moreover, the tsunami earthquakes occurred in the SSEs area in 1947 (Bell et al. 2014).
In the Hikurangi region, SSEs occur in the focal area of the tsunami earthquake, but what about the Kurile Islands?Since there is no APG network, it is not possible to directly measure seafloor deformation around the focal area of the earthquakes in 1969 and 1975.Instead of the geodetic measurement, I inferred SSEs based on the following idea that when SSEs occur, the stress is transferred in the surrounding faults, and consequently the long-term seismic activity is expected to be changed.There are previous studies that estimate SSEs from change in seismic activity (e.g., Ogata 2007;Kumazawa et al. 2010).Immediately after a megathrust interplate earthquake, the extensional stress in the direction of the dip of slab decreases within the descending oceanic plate near the downdip edge of megathrust seismic fault and the downdip-extension-type intraslab earthquakes decrease (Lay et al. 1989;Taylor et al. 1996).In this study, first, decrease in occurrence rate of earthquakes, i.e., seismic quiescence, was searched.Second, the fault model of SSE was constructed so that the spatial pattern of change in the Coulomb failure stress (ΔCFS) matches with that of the seismic quiescence.Finally, I show that an earthquake swarm was observed near the SSE fault and it is consistent with the SSE model.

Data
Earthquakes occurred between 1 January 1964 and 30 September 2019, with the body wave magnitude m b ≥ 5.0, and the depth of hypocenter h ≤ 60 km were selected from an earthquake catalog published by the International Seismological Center (ISC) (Fig. 1).Following Katsumata and Zhuang (2020), the study area is in a circle with a radius of 882 km centered at 46.71°N and 154.33°E.Katsumata and Zhuang (2020) confirmed that earthquakes with m b ≥ 5.0 and h ≤ 60 km were detected and located without fail in the study area between 1964 and 2019.Although the ISC is in the process of rebuilding the catalog (Storchak et al. 2017), I used old data in the present study to maintain the temporal homogeneity of the earthquake catalog.

Method
In this study, the PMAP method (Katsumata and Zhuang 2020) was used for investigating the seismic quiescence.Since Katsumata and Zhuang (2020) described the PMAP method in detail, the outline of the method is briefly explained here.Suppose that N earthquakes were observed during the period T years and let t i (1 ≤ i ≤ N) be the origin time of the i-th earthquake.Assuming that the earthquake occurrence follows the stationary Poisson process, the probability P N (0) that no earthquake occurs between t i and time t (t i < t < t i+1 ) is defined as follows: In the actual calculation, the procedure is as follows.First, the study area is divided into grids and N earthquakes are selected around a node.Second, P N (0) is calculated and search for the minimum value while changing N from 5 to 40.Finally, the P-value is defined at node (x, y) at time t as: where x ranges from 143 to 163°E with an interval of 0.1°, y ranges from 40 to 54°N with an interval of 0.1°, and t ranges from 1977.7 through 2019.7 with an interval of 0.1 years.When selecting N earthquakes around the node (x, y), select the earthquake inside a circle with a radius of r km centered on the node.The higher the seismicity, the smaller the r.Since r designates the spatial resolution in search for the seismic quiescence, this circle with radius r is defined as the resolution circle.If N = 5 and r > 50.0 km, the P-value was not calculated at the node.

Results
The P-value was calculated in the study area at 1,426,348 nodes, i.e., 3388 spatial × 421 temporal nodes.From the P-value maps calculated every 0.1 years from 1977.7 to 2019.7,four representative times were shown in Fig. 2. It is obvious that small P-values are rare (see Additional file 1: Fig S1).For example, there are more than 30,000 nodes with a P-value of 0.01 or less, and they are frequently observed.The number of nodes with a P-value (1) (2) P x, y, t = min 5≤N ≤40 P N (0) of 0.00002 or less has decreased to 295 and is divided into three clusters.As a result, I found two nodes 1 and 1' with P = 1.57× 10 -6 , which is the smallest value among the P-values calculated (see Additional file 1: Table S1).Since the nodes 1 and 1' are close each other, it is appropriate to recognize them as the one seismic quiescence.Hereafter, only the node 1 will be discussed.The number of earthquakes included in the resolution circle of the node 1 is N = 36 and the distribution of epicenters is plotted on the map in Fig. 3.The seismic quiescence area is defined as an area in the resolution circle.The 36 earthquakes occurred between 2 March 1978 and 10 February 2004 and the occurrence rate is 1.4 events/year.No earthquake was observed between 10 February 2004 and 30 September 2019 and this period was recognized as the seismic quiescence period (see Additional file 1: Figs.S2,  S3).The epicenters in the seismic quiescence area were compared with the coseismic slip distribution of the past great earthquakes presented by Ioki and Tanioka (2016).In the case of the 1975 tsunami earthquake, the epicenters in the seismic quiescence area 1 are distributed around the downdip edge of the subfaults that have a large slip during the main shock.On the other hand, the The second smallest P-value was 4.79 × 10 -6 observed at the node 2 and the third smallest P-value was 1.67 × 10 -5 observed at the node 3 (see Additional file 1: Table S1).Comparing with the coseismic slip distribution of the largest aftershock of the 1963 Kurile earthquake presented by Ioki (2013), both two seismic quiesce areas 2 and 3 are included in the focal area of the 1963 aftershock.It is noteworthy that both the 1963 aftershock and the 1975 earthquake were tsunami earthquakes.

Statistical significance of seismic quiescence
The statistical significance of seismic quiescence was estimated by a numerical simulation using earthquake catalogs created by the ETAS model.The simulation procedure is as follows.First, a synthetic earthquake catalog including background and cluster activities is produced by assuming the ETAS parameters obtained in this study.Second, after the declustring, the P-values were calculated, which is the same analyses as those applied to the actual earthquake catalog.Finally, the minimum P-value is searched among the P-values calculated.This simulation procedure is repeated 1000 times and the distribution of the minimum P-value was obtained (see Additional file 1: Fig. S4).Observing the entire Kurile Islands for 42 years, the seismic quiescence of P ≈ 0.0001 is not unusual according to the simulation result.Whereas log 10 P = − 5.804 for the node 1 and the number of cases with log 10 P ≤ − 5.804 is 39, therefore the rate of by-chance is 39/1000 = 3.9%.Since the rate of by-chance is smaller than 5% at the node 1, the seismic quiescence 1 is not likely to occur by chance.On the other hand, the rate of by-chance is rather large at the nodes 2 and 3, and the seismic quiescence 2 and 3 is less significant statistically than the seismic quiescence 1.

Inferences about slow slip event near the trench axis
Hypocenter relocation using HYPODD For constructing SSE models based on the seismic quiescence, it is very important to use the reliable location of hypocenters.Therefore, the hypocenters in the seismic quiescence area 1 were relocated by a double-difference earthquake location method.Since the seismic quiescence 2 and 3 have low statistical significance, only the seismic quiescence 1 is considered here.The location method was developed by Waldhauser and Ellsworth (2000), known as HYPODD.There are previous studies that used HYPODD to determine hypocenters using data from worldwide seismograph networks (Waldhauser and Schaff 2007;Pesicek et al. 2010).A subroutine of HYPODD calculates traveltimes assuming a horizontally layered velocity structure.In this study, the subroutine was replaced suitably by the iaspei-tau package (Kennett and Engdahl 1991;Snoke 2009) which is a program to calculate traveltimes assuming a radially stratified velocity structure of the Earth.The arrival times of P-wave in the ISC bulletin were used, which were observed at 273 seismographic stations within 5000 km from the seismic quiescence area 1 (see Additional file 1: Fig S5).The maximum separation of event pairs was 30 km and then 4678 double-difference data sets were produced from the hypocenter and the origin time in the ISC bulletin.The standard 1-D Earth model iasp91 was assumed as the P-wave velocity structure.As a result, the rms of residuals of the 36 earthquakes reduced 21% from 1.67 to 1.32 s and the hypocenters shifted slightly (see Additional file 1:  (Ioki and Tanioka 2016).Twelve squares between 148 and 152°E are subfaults of the largest aftershock of the 1963 Kurile earthquake (Ioki 2013).The subfaults in grey have a coseismic slip larger than 1 m.Three circles numbered as 1, 2, and 3 are the seismic quiescence areas on Additional file 1: Table S1.Red crosses show the epicenter of earthquake that occurred during times other than the seismic quiescence period within each area.The solid lines along the Kurile Islands indicate the trench axis (Bird 2003).b The coseismic slip of the 1969 Kurile earthquake (Ioki and Tanioka 2016) was compared with the seismic quiescence.The notation of symbols is the same as a HYPODD.The fault parameters of the two models are listed on Additional file 1: Table S2.As mentioned in the introduction, many previous studies reported deep SSEs at the downdip limit of the seismogenic zone and some authors suggested the existence of shallow SSEs near the trench axis.We compare the deep and the shallow SSE models in the present study.The model 1 represents the shallow SSE on the plate boundary near the trench axis.
The fault area of the model 1 corresponds to the subfaults 6 and 7 of the 1975 tsunami earthquake defined by Ioki and Tanioka (2016) and they found a large coseismic slip in these subfaults.The model 2 represents the deep SSE on the plate boundary which is an extension of the 1969 earthquake fault in the deep direction.The fault slip was assumed to be 1 m for the models 1 and 2. Since the purpose is to compare positive and negative spatial patterns of ΔCFS, any amount of slip can be used here.The frictional coefficient was assumed to be 0.4 for the models 1 and 2. Earthquakes that are affected by ΔCFS are called as receivers.Christova (2015) conducted a stress inversion using focal mechanisms of intraslab earthquakes shallower than 60 km within the Pacific plate in the southern Kurile Islands and found the σ 1 -axis with strike = 127° and dip = 20°, and the σ 3 -axis with strike = 288° and dip = 69°.The σ 1 and σ 3 are the maximum and minimum principal stresses, respectively.In this region, the σ 3 is predominant in the downdip direction along the subducting plate.When SSEs take place, σ 1 and σ 3 either decrease, increase, or remain unchanged, depending on the location.To estimate the change in σ 1 and σ 3 , I assumed that the σ 1 -and σ 3 -axes match the P-and T-axes of the focal mechanism of the receiver, respectively.Consequently, the focal mechanism of the receiver for calculating ΔCFS was assumed to be strike = 31°, dip = 65°, and rake = 93°.Assuming the auxiliary plane of the focal mechanism, the result is almost the same.The centroid moment tensor (CMT) solution for 16 out of 36 earthquakes are shown in Additional file 1: Fig S7 .These CMT solutions were determined by the global CMT project (Dziewonski et al. 1981;Ekström et al. 2012).

Estimation of the SSE fault model
For the model 1, the area of the seismic quiescence 1 corresponds to that of negative ΔCFS (Fig. 4).This result means that, in the seismic quiescence area 1, earthquakes has been occurring at almost constant rate before 2004, this area experienced negative ΔCFS due to the SSE, and no earthquake has occurred since 2004.This characteristic pattern matches not only on the plan view but also on the cross-sectional view.For the model 2, as in the model 1, the area of the seismic quiescence 1 corresponds to that of negative ΔCFS.Although the other negative ΔCFS area appears around the deep edge of the SSE fault, no seismic quiescence was detected there.Therefore, the seismic quiescence 1 is better-explained by the model 1 rather than the model 2.
Earthquake swarm, increased seismicity, and seismic quiescence following the 1975 tsunami earthquake Some previous studies reported that SSEs were accompanied by earthquake swarm (e.g., Vallée et al. 2013;Ozawa et al. 2019), and there are attempts to estimate SSEs based on the swarm activity (Llenos and McGuire 2011;Marsan et al. 2013).According to the earthquake catalog of the Japan Meteorological Agency, I found that 43 earthquakes with M ≥ 2.5 occurred in February 2003 in a rectangular area near the deep edge of the SSE fault of the model 1 (Fig. 5).After that, the seismicity in the rectangular area was gradually decayed and 15 earthquakes occurred again at one month in December 2019.The northwestward migration of epicenters was observed between 2004 and 2006, showing in Fig. 5e.Considering the swarm activity, it is reasonable to think that the SSE started around February 2003.However, if so, there needs to be a reason why the start of seismic quiescence was delayed by one year.One possibility is that the rupture velocity of the SSE was very slow, and it takes a year that the stress falls below the threshold of earthquake occurrence.
An increase in seismicity rate was systematically searched based on the same ISC data as describded in the Section "Data" of this paper.As a result, I found that the seismicity rate increased in the red color area, where ΔCFS is positive, at the same time of decrease in 2004 (see Additional file 1: Fig S8).The seismicity rate is 0.590 events/year before 2004.0 and 0.754 events/year after 2004.0,indicating 1.28 times increase.Although this increase is not significant statistically, this result is not inconsistent with the SSE hypothesis presented in this paper.
If the seismic quiescence was caused by the SSE near the trench axis, the 1975 tsunami earthquake should be also followed by long-term seismic quiescence.To confirm this hypothesis, I searched for seismic quiescence systematically based on the same ISC data as described in the Section "Data" in this paper.As a result, I found that no earthquake occurred around the deeper edge of the seismic fault after the 1975 tsunami earthquake for approximately 5 years from 1975.5 to 1979.8 (see Additional file 1: Fig S9).However, note that this quiescence is not significant statistically.

Discussion
In general, the trench axis of subduction zone is far from GNSS networks on land, therefore it is difficult to detect SSEs that occurred near the trench axis.There are only a few cases that SSEs near the trench axis have been revealed by geodetic observations on the seafloor (Davis et al. 2015;Wallace et al. 2016;Araki et al. 2017).Instead of the direct measurement on the seafloor, Yamashita et al. (2015) argued the earthquake swarm as evidence for SSEs near the trench axis.Compared to the geodetic observation, since the quantitative relationship between ΔCFS and the occurrence rate of earthquakes is not clear, estimation of SSEs based on the seismic quiescence cannot quantitatively constrain some model parameters, e.g., the amount of fault slip.Despite of poor quantitative estimation, it might be an advantage that it is simple and inexpensive.
From physical point of view, there is a fundamental question whether shallow SSEs are allowed to occur near the trench axis.Yoshida and Kato (2003) conducted a numerical simulation based on rate and state dependent friction laws with assuming a simple model when the unstable area of (a-b) < 0 and the conditional stable area of (a-b) > 0 are adjacent along the upper boundary of the subducting plate.When the seismic slip occurs in the unstable area, the seismic slip occurs in the conditional stable area almost at the same time.In the unstable area, no aseismic slip is observed during the interseismic period.On the other hand, in the conditional stable area, episodic aseismic slips are observed during the interseismic period when the same amount of stress is accumulated as the stress drop during the main shock.The timing of SSE occurrence is consistent with the SSE model proposed in this study.If the 2003 earthquake swarm is the beginning of the SSE, In this area, the Pacific plate is subducting beneath the North American plate at a rate of 0.08 m/year (DeMets et al. 1994), and thus, the slip deficit of 2.24 m at most accumulates in 28 years.This amount of the slip deficit is comparable with the slip of 1.6-2.2m during the 1975 tsunami earthquake estimated by Ioki and Tanioka (2016).A model in which there is a conditional stable area near the trench axis and an unstable area adjacent to it in the deeper part along the plate boundary is presumed to be quite common in subduction zones (Lay et al. 2012).Therefore, it is likely that SSEs occurring near the trench axis are common phenomena in subduction zones.

Conclusions
In the present study, the background seismicity including neither aftershock nor earthquake swarm has been investigated in the Kurile Islands subduction zone and unusual long-term seismic quiescence was found.I argue that the stress release associated with the shallow SSE near the trench axis is the most likely model to explain the seismic quiescence.Many previous studies found such kind of the long-term seismic quiescence which preceded megathrust earthquakes and has been interpreted as a reliable precursor (e.g., Katsumata 2011).According to a retrospective forecast experiment based on the seismic quiescence, when the alarm is issued to the 40% area of the whole area, the prediction rate is about 80% and this means that most of the seismic quiescence area will not experience subsequent large earthquakes (Katsumata and Nakatani 2021).I suggest that the seismic quiescence found in this study is also not to be a precursor of great earthquakes, but to be the one caused by SSEs, which frequently occur near the trench axis during the interseismic period.

Fig. 1
Fig.1Earthquakes considered in the present study.a, c Earthquake distributions before and after the declustering process, respectively.Earthquakes were selected from the ISC catalog and occurred from 1 January 1964 to 30 September 2019, with m b ≥ 5.0 and depth ≤ 60 km.The study area is in a large circle with a radius of 882 km centered at 46.71°N and 154.33°E, and the solid line along the Kurile Islands indicates the trench axis(Bird 2003).b, d Space−time plots before and after the declustering process, respectively

Fig. 2
Fig. 2 Time slices of the P−value distribution at year of (a) 2010.0,(b) 2015.0,(c) 2018.0, and (d) 2019.7.The nodes are not colored if the radius of the resolution circle is larger than 50 km.The smaller the P − value, the more significant the seismic quiescence.Three seismic quiescence areas are found and numbered as 1, 2, and 3 in (c) and listed on Additional file 1: TableS1.The solid line along the Kurile Islands indicates the trench axis(Bird 2003)

Fig. 3
Fig. 3 Seismic quiescence areas and the past large earthquakes.a Eight squares between 146 and 149°E are subfaults of the 1975 tsunami earthquake(Ioki and Tanioka 2016).Twelve squares between 148 and 152°E are subfaults of the largest aftershock of the 1963 Kurile earthquake(Ioki 2013).The subfaults in grey have a coseismic slip larger than 1 m.Three circles numbered as 1, 2, and 3 are the seismic quiescence areas on Additional file 1: TableS1.Red crosses show the epicenter of earthquake that occurred during times other than the seismic quiescence period within each area.The solid lines along the Kurile Islands indicate the trench axis(Bird 2003).b The coseismic slip of the 1969 Kurile earthquake (Ioki and Tanioka 2016) was compared with the seismic quiescence.The notation of symbols is the same as a

Fig. 4
Fig. 4 The Coulomb failure stress change (ΔCFS) caused by the fault motion of slow slip event (SSE).Rectangles indicate the SSE faults of (a) model 1 and (c) model 2. Crosses depict 36 earthquakes in the seismic quiescence area 1, which have been relocated by using HYPODD.a, c are plan views.Colors represent ΔCFS at a depth of 40 km.Solid lines along the Kurile Islands indicate the trench axis (Bird 2003).6 and 7 in (a) are the subfault number defined by Ioki and Tanioka (2016).b, d are vertical cross−sectional views along the broken line AB on the plan view.Gently curved black lines in (b) and (d) represent the upper boundary of the subducting Pacific plate (Nakanishi et al. 2004).Bold lines on the upper boudary of the subducting Pacific plate show the SSE faults assumed in this study

Fig. 5
Fig. 5 Earthquake swarm activity.a Thin lines show a northeastern corner of the SSE fault assumed in Fig. 4a.A rectangle drawn by a bold line indicates the area where the earthquake swarm activity was observed.Dots within the rectangle are the epicentres determined by the Japan Meteorological Agency between 2000 and 2022.b The cumulative number of earthquakes within the rectangle in (a).(c) Magnitude−Time plot of the earthquakes within the rectangle in (a).d Space−Time plot along AB in (a).e Space−Time plot along CD in (a).Broken lines in red color show the migration of earthquakes from 2004 to 2006