Evolution of aseismic slip rate along plate boundary faults before and after megathrust earthquakes

Similar earthquakes that occur in approximately the same location have the potential to reveal the spatio-temporal changes in aseismic slip along plate boundaries. Here we identify similar earthquakes with moderate magnitudes that occurred worldwide between 1989 and 2016 by using seismograms recorded by the Japanese dense seismic network. The slip rate along the plate boundaries estimated from similar earthquakes increased rapidly following M > 8 megathrust ruptures and then gradually decayed over periods of ~10 years, which correlates with after-slip progressing around the source areas. More than 30 years after large megathrust earthquakes, the slip rate begins to show a gradual increase. This gradual increase in slip rate after the decay may be due to an increase in stress levels that accumulate during tectonic loading. The spatio-temporal characteristics of inter-plate aseismic slip can be used to provide a valuable framework for understanding the long-term evolution of slip-rate during megathrust earthquake cycles. The aseismic slip rate, which increases sharply immediately after a giant earthquake and then rapidly decays, increases again gradually toward the next event, according to observations of earthquakes globally using the Japanese dense seismic network.

A n earthquake is a sudden fast slip event that propagates dynamically along a fault and releases elastic energy accumulated due to a locking of the fault during the interseismic period. Growing evidence of geophysical observations has demonstrated that subduction zone faults host a broad spectrum of slip modes, from slow to fast slip 1,2 . This recent view is quite different from the classic notion of allowing only two discrete states of fault: locking and fast slip.
Earthquakes with similar waveforms are typically considered to occur at approximately the same location but at different times, and are referred to as similar earthquakes. Seismic observations made over several decades in Japan have shown that many similar earthquake sequences in the subducting plates are caused by repeated slips in small patches surrounded by creeping aseismic regions along with the plate boundary [3][4][5] . These sequences are also referred to as repeating earthquakes due to the recurring nature of seismic energy release, and they yield information on the spatio-temporal evolution of fault slip rate (i.e., creeping rate).
Recently, similar earthquakes have been detected in many regions around the world, based on similarities in seismogram characteristics and the inter-event distance of each earthquake pair. This includes subducting plate boundaries in Japan 4 , the Kuril-Kamchatka-Aleutian Islands, and the South Sandwich Islands 6 , the Tonga-Fiji-Solomon region 6,7 , Sumatra 8 , and Chile 9,10 . In addition, similar earthquakes have been detected in the plate collision zone in Taiwan 11 , along with transform faults in California 12 , from plate interiors [13][14][15][16][17][18][19] , in mines 20 , and at intermediate-deep depths in subducting slabs 6,21 .
The recurrent nature of similar earthquakes over relatively long periods of time can be used to directly analyze spatio-temporal changes in slip rate along inter-plate boundaries at depth 4 . This is normally difficult to achieve from surface geodetic measurements, especially in offshore regions that contain few geodetic stations. In addition, changes in seismogram characteristics at the time of recurrence can be used to investigate temporal changes in seismic velocity within the Earth 22,23 . In this study, we detected moderate-sized similar earthquakes around the world by calculating the cross-correlation coefficients of bandpass-filtered seismograms recorded by the Japanese dense seismic network. Based on these data, we investigated the spatial distributions and temporal characteristics of inter-plate aseismic slip at a number of subduction zones around the world.

Results
Spatial distribution of similar earthquakes. We used seismic waveform data collected by dense nationwide seismic networks located across the Japanese islands (Fig. 1). To extract similar earthquakes, we calculated the cross-correlation coefficient of the vertical component in bandpass-filtered seismograms at each station (see 'Methods'). The hemisphere map centered on Japan in Fig. 1 shows the distribution of the cross-correlation coefficients obtained in this study. We identified numerous similar earthquakes that occurred along subducting plate boundaries, including the Andaman-Sumatra-Java, the Solomon-Vanuatu-Tonga, and Japan-Kuril-Kamchatka-Aleutian subduction zones. These similar earthquakes are linked to aseismic fault slip along shallow inter-plate boundaries at depths of <100 km. The similar earthquakes are distributed heterogeneously, even within single subduction zones that host megathrust earthquakes. Similar earthquakes were also identified at depths of >100 km in regions where intra-slab earthquakes occur frequently (e.g., the Solomon, Tonga, and Izu-Ogasawara regions). These deeper events might have been detected due to relatively simple source processes, the relatively homogeneous structure near the source, and a nearly identical path to radiated seismic waves.
In contrast, similar earthquakes were not detected in China, although previous studies have identified pairs of similar events in this region 14 . The contrasting results of the present and previous studies may have been caused by a difference in the analyzed seismic phase and frequency band.
Average slip rate along the subducting plate boundary. Based on the spatio-temporal evolutions of similar earthquakes that occurred during the analysis period, we demonstrate that interplate similar earthquake sequences occurred more frequently after megathrust earthquakes. Other sequences of similar earthquakes have a constant recurrence interval of several years to several decades. Pairs of similar earthquakes with the longest recurrence interval (~20 years) occurred in the Aleutian Islands. These sequences may reflect spatio-temporal changes in inter-plate aseismic slip, and we, therefore, attempt to analyze such changes in subduction zones around the Pacific and Indian oceans based on similar earthquakes extracted in this study.
To estimate the aseismic slip history, we assume that the slip rate during the inter-seismic period of similar earthquake sequences is equal to the slip rate of surrounding creeping areas (see 'Methods'). If the moderate-sized similar earthquakes analyzed in this study have recurrence intervals that relate to the plate convergence rate, the recurrence interval at the plate boundaries with low average slip rates will exceed 10 years. The resolvable lower limit of the slip rate that can be estimated from the similar earthquakes in the 27-year interval used in this study is about 20-30 mm yr −1 , depending on the recurrence time and magnitude of the similar earthquakes. Figure 2a shows the average slip rates in subduction zones estimated from similar earthquakes that occurred at depths of <100 km. The average slip rates vary laterally between different subduction zones, and in some cases along the length of the same subduction zone. Figure 2b shows the relationship between the relative velocity of the two converging plates, taking into account back-arc spreading 24 ; the estimated slip rates in each region are shown in Fig. 2a. The estimated slip rates in several areas are higher than the relative velocity of the two converging plates. These areas correspond to regions in which postseismic slip has occurred recently following giant earthquakes, including the 2004 Sumatra-Andaman earthquake (Mw9.1) in area 12, the 2005 northern Sumatra earthquake (Mw8.6) in area 13, the 2006 Kuril earthquake (Mw8.3) in area 2c, the 2007 southern Sumatra earthquake (Mw8.4) in area 14, and the 2011 Tohoku-Oki earthquake (Mw9.1) in area 2f. These areas show large error bars due to the spatio-temporal variation of the slip rates linked to the postseismic slip (Fig. 2b). The increase of slip rate in area 5 near the Kanto District in Japan was induced by the 2011 Tohoku-Oki earthquake 25 . The recurrence intervals of similar earthquake sequences decreased significantly following these recent giant earthquakes, resulting in an increase in the average slip rate.
In contrast, we recognized several areas where the slip rates are much lower than the relative plate velocities, including the Kuril-Kamchatka-Aleutian (areas 1-2), Nankai-Ryukyu (areas 6-7), Manila (area 9), Philippines (area 10), and Solomon (areas 16-20) subduction zones. These observations indicate that the plate boundary faults in these regions are strongly locked, and that elastic strain energy has been accumulating within each of these subduction zones.
Slip rates in some areas are almost equal to the rate of relative plate motion, including the northern extents of the 1952 and 2011 M9 earthquakes, the Izu-Bonin-Mariana region (areas 3-4), and parts of the Solomon subduction zone (areas [16][17][18][19][20]. These regions may host episodic aseismic slip events, including slow slip events and/or creeping movement. The Tonga subduction zone (area 22) has a high rate of relative plate motion and also shows relatively high rates of fault slip. Back-arc spreading has occurred in the Ryukyu (area 7), Mariana (area 4), and Tonga (area 22) areas. The slip rates estimated from similar earthquakes are consistent with the relative velocities of plate motion, which includes the back-arc spreading rate.
These results suggest that the spatio-temporal changes in interplate aseismic slip in many areas around the world can be estimated from similar earthquake sequences. It might not be possible to detect similar earthquakes in areas with slip rates lower than the detection threshold (20-30 mm yr −1 ). Previous studies have suggested that the source regions of large earthquakes show notable slip deficits during the inter-seismic period 4 . This period may be represented by some areas of the Aleutian (area 1), Kuril-Kamchatka (areas 2b, d), and Solomon-Vanuatu (areas 17, 19-21) subduction zones that are currently strongly locked.  Short-term history of slip rate. The average slip rates shown in Fig. 2 indicate slip acceleration in source regions that have recently hosted megathrust earthquakes. However, the similar earthquakes have recurrence intervals of years to decades, which are too long to investigate temporal variations in short-term slip rates. To track such short-term variations, we compute "short-term slip rates" assuming that the slip rate of each recurring similar earthquake is perturbed according to temporal changes observed in the regular seismic activity in adjacent regions (see 'Methods'). Figure 3 shows the short-term histories in three areas near the hypocenters of the recent M8-9 class megathrust earthquakes. Immediately after each megathrust earthquake, we detect clear afterslip in the source regions. The areas with afterslip coincide approximately with those estimated from geodetic data analysis 26-28 , although it is not possible to detect all of the afterslip areas due to the limited occurrence of regional moderate-sized earthquakes. The duration of afterslip following M8 class events was ∼2-5 years. On the other hand, as of February 2016, the areas affected by the 2004 Sumatra-Andaman earthquake (Mw9.1) and the 2011 Tohoku-Oki earthquake (Mw9.1) had higher slip rates than the relative plate motions along each subduction boundary. The great earthquakes have a higher maximum slip rate and longer duration of afterslip. The estimated slip rate was highest in the western part of the epicentral area of the Tohoku-Oki earthquake, namely the down-dip side of the ruptured area. The amount of slip that occurred within an area of 1.5°by 1.0°was 230 cm over 3 years, which is consistent with estimates derived from geodetic data 28 .

Discussion
The long-term trend in slip rates, measured over decadal scales, seems to increase several years before megathrust earthquakes (Fig. 3), indicating an unlocking of the plate boundary fault. Indeed, geodetic measurements made using onshore GNSS networks reveal that long-term aseismic creep on the plate interface had occurred down dip and to the south of the large coseismic slip zone of the 2011 Tohoku-Oki mainshock rupture 29,30 . Foreshock sequences including slow slip transients have been reported over timescales of several weeks immediately preceding the Tohoku-Oki earthquake 31 . In contrast, our result for the area near the epicenter (areas iii and v in Fig. 3) shows that the slip rates increased for about 5-8 years before the mainshock, but then tend to decrease ∼1-3 years before them. Because the recurrence interval of similar earthquakes analyzed in this study (magnitude of ≥4.5) is longer than several years, it is difficult to determine the short-term temporal evolution of slip before megathrust earthquakes. To improve the temporal resolution of slip rate using this type of analysis, we require additional observational data on the waveform similarity of earthquakes with smaller magnitudes. Such an improved catalog would strengthen our understanding of the relationships between spatio-temporal changes in slip rate and megathrust earthquakes.
Nevertheless, the present catalog provides a valuable framework for understanding the long-term evolution of slip rate during the earthquake cycle. Figure 4 shows a time series analysis of stacked slip rates plotted against the elapsed time since the 11 megathrust earthquakes with Mw ≥ 8.3 that occurred after 1923. The slip rates calculated for each subduction zone were smoothed using a 3-year moving average, and normalized to the plate convergence rate in the area (see 'Methods'). The analysis shows that following rupture of each earthquake, the slip rates initially increased rapidly and then gradually decayed over periods of~10 years, which correlates with afterslip. After~30 years, the slip rate begins to increase as we approach the next megathrust rupture (blue shaded curve in Fig. 4). This gradual increase in slip rate during the middle and later stages of the earthquake cycle may be due to an increase in stress levels that accumulate during tectonic loading. Numerical simulations of megathrust rupture show that aseismic sliding at a slip velocity of 30-60% of the plate convergence rate starts midway through the earthquake cycle and continues until megathrust rupture in the afterslip area where many similar earthquakes occur 32 . The local heterogeneity of slip rate within the source region may also contribute to the stress loading of the main locked area.
The absence of similar earthquakes in some areas during the analysis period suggests that such areas are either completely locked or experience quasi-static slip. Source regions of giant earthquakes that occurred before the analysis period, such as the 1952 Kamchatka earthquake, have large slip deficits. In regions where no similar earthquakes were detected but other regular earthquakes continued to occur, the slip rate is assumed to be extremely low. Areas that experienced large coseismic slips during megathrust earthquakes show very low numbers of regular earthquakes. This means that seismic potential can be assessed in previously ruptured areas or in regions that are surrounded by seismogenic areas. According to the time predictable model, an earthquake will occur when the amount of accumulated slip is equal to the slip during the previous large earthquake. Even if the slip predictable model is applied, the amount of slip produced by a potential earthquake can be estimated by tracking the cumulative slip deficit.
We believe that additional similar earthquakes will be detected by extending the analysis periods in areas where similar earthquakes were not identified in this study. A comparable analysis of small-magnitude earthquakes using local seismic networks could also be used to estimate inter-plate slip at a higher spatial and temporal resolution. A comprehensive database of local and global similar earthquakes would enable comparison of interplate aseismic slip characteristics along various plate boundary zones around the world. Furthermore, such a database could be used to extract information on spatio-temporal changes in the seismic velocity structures beneath seismic observation networks in Japan and around the world.

Methods
Data. We used seismic waveform data collected by dense nationwide seismic networks located across the Japanese islands. The seismic networks have been operated consistently over several decades by the National Research Institute for Earth Science and Disaster Resilience (NIED), the Japan Meteorological Agency (JMA), and Earthquake Research Institute (ERI) at the University of Tokyo (Fig. 1). This represents the most suitable seismic observation network in the world in terms of detecting and analyzing the long-term activity of similar earthquakes. Using the worldwide hypocenter catalog of the United States Geological Survey, we focused on events with magnitudes of ≥4.5 at distances of <90°that occurred from 15 September 1989 to 29 February 2016. For each event, we selected waveform data with a high signal-to-noise ratio around the onset of the P-wave arrival. A total of 1340 seismic stations were available and 37,341 target events were analyzed.
Extraction of similar earthquakes. To extract similar earthquakes, we detected earthquake pairs with epicenter separations and differences in back azimuth of <3°. For each earthquake pair, we used the method of Igarashi 4 to calculate the crosscorrelation coefficient of the vertical component in bandpass-filtered seismograms at each station. The time window used to compute cross-correlation coefficients was 80 s from the P-wave onset. We applied multiple bandpass filters to the waveforms, with passbands of 0.1-0.5, 0.2-1.0, 0.5-2.0, and 1.0-4.0 Hz, and computed the cross-correlation coefficient at each filter band. We chose the highest cross-correlation coefficient from those in the analyzed filter bands as a representative value for each station. The median of the correlation values over the nationwide seismic network was adopted as the cross-correlation coefficient of the earthquake pair, following Aster and Scott 33 .
Considering multiple sources of scattering during the propagation of seismic waves, the similarity of waveforms decreases exponentially with an increasing inter-event distance 12,34 . Thus, earthquakes sharing almost identical waveforms are considered to have occurred in nearly the same location. We identified similar earthquakes based only on waveform similarities, because we could not obtain information on the relative hypocenter locations that are needed to determine the overlapping source areas within each earthquake pair.
When the median of the cross-correlation coefficient from an earthquake pair exceeded 0.85, it was identified as a similar earthquake pair. This threshold was selected by comparison with the existing catalogs of similar earthquakes compiled in the previous studies 8 . If several earthquake pairs shared the same event, they were grouped together as within the same sequences of similar earthquakes. A total of 3579 similar earthquakes (from 1125 sequences) that occurred over a 26 year period were extracted. The complete list of similar earthquakes is provided in Supplementary Data 1.
Estimation of slip rate. To estimate the aseismic slip history, we applied the method of Igarashi 4 to our catalog of similar earthquakes. This method assumes that the slip rate during the inter-seismic period of similar earthquake sequences is equal to the slip rate of surrounding creeping areas, which is calculated based on the amount of slip and the recurrence interval of each similar earthquake. The amount of slip during each event was calculated from an empirical relationship between seismic moment (Mo) and slip (d) 35 , where log(d) = 0.17 log (Mo) -2.36. The scalar moment is calculated from magnitude using the relationship between magnitude and scalar seismic moment 36 . The slip rate between events can be obtained by dividing the amount of slip during the next similar earthquake by the recurrence interval according to the slip-predictable model 37 . If an expected previous or next event falls outside the analysis period, we extrapolated by using the lowest value of the average or adjacent slip rate. We assume that sequences occurring in almost the same location indicate aseismic slip in the same area. To display the spatial distribution of slip rates, we calculated the average slip rate and plotted this at the center of each similar earthquake pair. A total of 619 time series of slip rates were obtained in the present study. The list is provided in Supplementary Data 2. The average slip rate, relative plate motion, and their standard deviation in each subduction zone area were calculated from the slip rate for each group obtained here and the relative plate motion at their occurrence.
To track short-term variations, we use the method outlined in the previous section but also assume that the slip rate of each recurring similar earthquake is perturbed according to temporal changes observed in the regular seismic activity in adjacent regions. We calculated the short-term history of the slip rate by dividing the amount of slip in the next similar earthquake by the frequencies of nearby seismic events that occurred during the inter-seismic period of similar earthquake sequences. We selected seismic events with a magnitude of ≥4.5 and hypocentral distances of <30 km from each similar earthquake sequence. Both methods result in the same cumulative slip if the time window is sufficiently long compared with the recurrence interval of each similar earthquake sequence.
To investigate the long-term evolution of slip rate during the earthquake cycle, we calculated the slip rate against the elapsed time since the 11 megathrust earthquakes with Mw ≥ 8.3 that occurred after 1923. The list is provided in Supplementary Table 1. As of 2016, 93 years have elapsed since the oldest megathrust event occurred. We first calculated the slip rates in each plate boundary area smoothed using a 3-year moving average and normalized by the plate convergence rate in the area (pink bars in Fig. 4). Then, the time series for the same period were stacked (red bars in Fig. 4). With the information on the time series of slip rates, we outlined the long-term evolution of slip rate (blue curve in Fig. 4). The curve between 10 and 30 years of the time elapsed from megathrust rupture is somewhat uncertain due to the lack of data. To fill this gap, it is required to track the time-series of slip rate using similar earthquake sequences within the next decades.

Data availability
Seismic waveform data used in this study are available from the Data Management Center of the National Research Institute for Earth Science and Disaster Resilience (NIED) 38,39 (https://hinetwww11.bosai.go.jp/auth/?LANG=en), the Japan Meteorological Agency (JMA) (http://www.data.jma.go.jp/svd/eqev/data/bulletin/ index_e.html, http://www.data.jma.go.jp/svd/eqev/data/daily_map/index.html), and the Earthquake Research Institute, the University of Tokyo (http://tkypub.eri.u-tokyo.ac.jp/ harvest). The hypocenter catalog used in this study is available from the United States Geological Survey (http://earthquake.usgs.gov/earthquakes/search). Additional data related to this paper are available in Supplementary Data 1, Supplementary Data 2, and Supplementary Information.

Code availability
The Generic Mapping Tools 40 created figures are obtained from https://www.genericmapping-tools.org/. The analysis codes used in this study are available from the corresponding author upon request.