Research ArticleGEOPHYSICS

Spatial and temporal seismic velocity changes on Kyushu Island during the 2016 Kumamoto earthquake

See allHide authors and affiliations

Science Advances  24 Nov 2017:
Vol. 3, no. 11, e1700813
DOI: 10.1126/sciadv.1700813


Monitoring of earthquake faults and volcanoes contributes to our understanding of their dynamic mechanisms and to our ability to predict future earthquakes and volcanic activity. We report here on spatial and temporal variations of seismic velocity around the seismogenic fault of the 2016 Kumamoto earthquake [moment magnitude (Mw) 7.0] based on ambient seismic noise. Seismic velocity near the rupture faults and Aso volcano decreased during the earthquake. The velocity reduction near the faults may have been due to formation damage, a change in stress state, and an increase in pore pressure. Further, we mapped the post-earthquake fault-healing process. The largest seismic velocity reduction observed at Aso volcano during the earthquake was likely caused by pressurized volcanic fluids, and the large increase in seismic velocity at the volcano’s magma body observed ~3 months after the earthquake may have been a response to depressurization caused by the eruption. This study demonstrates the usefulness of continuous monitoring of faults and volcanoes.


Seismic velocity in crustal rocks changes during earthquakes by several mechanisms related to, for example, fault zone damage, pore pressure and stress state changes (for example, pressure perturbation), and healing processes (18). Seismic velocity variations can often be explained by the presence of cracks; earthquake ruptures and associated abnormal pore pressure generate open cracks in the crust, and the generation of cracks changes the crust’s elastic properties, including seismic velocity (9, 10). Therefore, seismic velocity is a key geophysical property for characterizing dynamic processes and the status of deep faults. Further, previous studies have shown that seismic velocity associated with eruptions in volcanic areas varies greatly, mainly because of increased pore pressure around the magmatic body (2, 1113). In particular, under high pore pressure conditions, seismic velocity shows large variation because even minor increases in pore pressure under these conditions easily generate cracks, thereby influencing seismic velocity (9).

Several studies have observed changes of seismic velocity due to earthquakes or volcanic activity by using multiplet earthquake analysis (14, 15), controlled source experiments (5, 16), and seismic interferometry (13, 7, 8, 11, 17, 18). Although seismic velocity variation between several seismometer pairs has been reported using these approaches, temporal variation of seismic velocity along intraplate faults has yet to be mapped using many seismometer pairs. Here, by applying seismic interferometry to ambient noise data, we reveal the temporal and spatial variation of seismic velocity due to the 2016 Kumamoto earthquake and volcanic activity at Aso volcano in central Kyushu Island, Japan.

The 2016 Kumamoto earthquake took place in the central part of Kyushu Island, where dozens of Hi-net seismic stations (19) are deployed (Fig. 1). This earthquake comprised a series of temblors along the Hinagu-Futagawa fault system (2023). Their source mechanisms typically indicated right-lateral motion, consistent with the north-south tectonic extension of central Kyushu (24). The series began along the Hinagu fault system (oriented north-northeast–south-southwest) with a moment magnitude (Mw) 6.2 foreshock on 14 April (Fig. 1). The largest event of the series (Mw 7.0), the mainshock, occurred 28 hours later (16 April) along the Futagawa fault system (oriented northeast-southwest) (20, 22). The maximum displacement along the fault plane, modeled from surface deformation data, was ~6 m and was observed west of Mount Aso (22), and the depth of the fault rupture was shallower than 10 km. The mainshock rupture was halted at the Aso caldera complex. Two hours after the mainshock, a series of ruptures occurred successively farther northeast, beginning with a Mw 5.8 event near Mount Aso (Fig. 1). High volcanic tremor activity around Mount Aso was observed before the Kumamoto earthquake (25), and the higher activity level continued after the earthquake. Aso volcano erupted on 16 April and 1 May 2016, and these small eruptions were followed by a large eruption on 7 and 8 October 2016 (26). The spatial and temporal variation of seismic velocity derived by seismic interferometry can be used to estimate stress and pressure perturbation associated with the seismic activity. Therefore, we can infer that stress and pressure perturbation due to the 2016 earthquake are a likely cause to the other earthquakes and activate the volcano.

Fig. 1 Map of central Kyushu Island, Japan, with locations of Hi-net stations (yellow dots).

Hypocenters of large earthquakes and aftershocks associated with the 2016 Kumamoto earthquake are shown by green stars and white dots, respectively. The Aso caldera is outlined by a dashed orange line. The black lines connect the seismometer pairs shown in Fig. 2 and figs. S3 and S4. The base map is a 10-m mesh digital elevation model published by the Geospatial Information Authority of Japan. We drew this figure with Generic Mapping Tools (53).


To retrieve virtual seismograms of seismic waves propagating between pairs of stations (Fig. 2, A to C), we applied seismic interferometry (27, 28) to the ambient noise recorded at 36 Hi-net seismic stations (Fig. 1) from December 2015 to November 2016. We estimated daily variation of seismic velocity by applying a stretching interpolation technique (29) and moving-window cross-spectral (MWCS) technique (1, 11, 30, 31) to the coda wave on the virtual seismogram (see Materials and Methods). Although the relative velocity changes are calculated from both techniques, the results from the MWCS were instable mainly because of tremors at Mount Aso (see Materials and Methods). Therefore, we decided to use the stretching technique for further analyses. By assuming that surface waves were dominant in the coda (32), we used frequency-dependent depth sensitivity of surface waves in our interpretation (Fig. 3; see Materials and Methods). The surface wave velocity is closely related to the S-wave velocity (33); thus, we inferred that the seismic velocity variation estimated in this study (Fig. 2) was induced mainly by the S-wave velocity variation.

Fig. 2 Temporal variation of seismic traces and velocity between three seismometer pairs.

(A to C) Typical seismic traces (cross-correlations) between seismometer pairs: (A) across the fault plane from A to A′ in Fig. 1, (B) far from the fault (from B to B′ in Fig. 1), and (C) across Mount Aso (magmatic body; from C to C′ in Fig. 1). The vertical axis shows travel time in seconds, and the horizontal axis shows dates from December 2015 to November 2016. The red arrows indicate the date of the mainshock and volcanic eruption (16 April), and the yellow arrows indicate the date of the large eruption of Aso volcano (7 and 8 October). (D to F) Typical seismic velocity variation between station pairs derived from the seismic traces displayed in (A) to (C). Background color indicates the cross-correlation coefficient obtained by trace stretching; black curves show daily variations of the estimated velocity changes Embedded Image with respect to changes before the Kumamoto earthquake, defined at the maximum value of the coefficient; and dashed black curves indicate the SD of the velocity change estimation (see Materials and Methods). White dashed lines show the time window (30 days) influenced by the mainshock and the largest eruption.

Fig. 3 Depth-dependent S-wave sensitivity kernels (partial derivatives) for fundamental mode Rayleigh waves with respect to S-wave velocity, which constitute a proxy for depth resolution.

The sensitivity kernels were computed by using the DISPER 80 program (51) for a one-dimensional (1D) layered model near Aso volcano estimated by Nishida et al. (52). Sensitivity kernels were normalized by the maximum amplitude at 0.8 Hz.

To measure the seismic velocity change associated with the 2016 Kumamoto earthquake, we determined velocity changes relative to the averaged change over 4 months before the mainshock (Embedded Image; Fig. 2, D to F, and figs. S1 and S2; see Materials and Methods). To stabilize the results, we estimated the daily variation of the velocity change by moving the 30-day window. Although the temporal resolution of our estimated velocity variation is not high enough to correlate the variation with minor events (for example, small aftershocks), we were able to identify the first-order dynamic variation (for example, pressure variation) of the seismogenic faults and volcano during the entire earthquake sequence. Because the surface wave velocity has dispersion characteristics, the depth sensitivity is related to the frequency range of the analyzed ambient noise (Fig. 3; see Materials and Methods). Our results using the frequency range of 0.1 to 0.9 Hz were sensitive to S-wave velocity variations from the surface to a depth of ~10 km. Using the velocity changes measured between all pairs of stations with an interstation distance of <40 km (fig. S2), we mapped seismic velocity changes in the area around the seismogenic faults and Aso volcano, central Kyushu Island (Fig. 4 and movie S1; see Materials and Methods) (2).

Fig. 4 Spatial and temporal variation of seismic velocity in central Kyushu during the 2016 Kumamoto earthquake.

The surface wave velocity changes Embedded Image within each time window (30 days) relative to the averaged pre-earthquake value are displayed (see Materials and Methods). Each panel shows the central date within the 30-day window: (A) 8 March 2016, (B) 1 May 2016, (C) 14 May 2016, (D) 1 June 2016, (E) 1 August 2016, and (F) 1 October 2016. Warm colors indicate regions where seismic velocity was decreased. White dots are Hi-net stations.


Before the 2016 Kumamoto earthquake, we could not detect large velocity variation in most of central Kyushu (movie S1). Small seismic velocity fluctuations detected around Mount Aso before the earthquake may have been related to volcanic activity, including a small eruption. During the 2016 earthquake, seismic velocity around the Hinagu-Futagawa fault system and Mount Aso decreased greatly (Figs. 2, D and F, and 4B and movie S1). Because the depth range in which the variation of seismic velocity could be estimated (<10 km in depth) includes the depth of the fault plane with maximum slip (~5 km depth) (22), the velocity variation could be related to fault zone damage and stress and pressure changes around the deep rupture fault. Around the Hinagu-Futagawa fault system (Fig. 2D), seismic velocity was reduced by 0.3 to 0.4%. This velocity reduction is similar to that observed during the 2012 Costa Rica earthquake (Mw 7.6) in a subduction zone setting (~0.6%) (6) but is larger than the reduction during the 2004 Mw 6.0 Parkfield earthquake (<0.1%) (1). These velocity reductions around the seismogenic fault also depend on the interval of seismometer pairs. Because the interval between Hi-net seismometers displayed in Fig. 2D is ~20 km (Fig. 1), the seismic velocity reduction within the narrow fault damage zone may have been larger (>0.4%). Laboratory experiments that used a similar lithology (34) demonstrated that S-wave velocity at an effective pressure of ~50 MPa (corresponding to a depth of ~5 km) decreases by ~0.36% if the stress perturbation is ~1 MPa (35). Such a velocity reduction related to crack opening could thus account for the observed velocity reduction around the fault zone. Because seismic velocity around the stress accumulation area at the fault edge also decreased, the crack opening due to formation damage or stress perturbation (increase in pore pressure) could be a primal mechanism for the velocity reduction. Even far from the seismogenic faults and volcano, seismic velocity decreased slightly during the earthquake (Figs. 2E and 4 and movie S1), suggesting that pressure perturbation and surface damage due to earthquake vibration occurred widely in central Kyushu.

After the 2016 Kumamoto earthquake (from mid-April to November), seismic velocity changes around the Hinagu-Futagawa fault system reflected gradual healing (Fig. 4, C to F, and movie S1). Open cracks temporarily generated in the fault zone by the earthquake could have closed by several mechanisms, allowing seismic velocity to recover (36, 37). The short-term velocity recovery could be explained by crack closing as a result of a reduction in pore pressure. However, in most of the area around the seismogenic faults, healing was not complete and the pre-earthquake state was not attained during our observation period (7 months after the earthquake; Fig. 2D). Such a longer-term velocity reduction around the seismogenic faults could be caused by stress change or damages related to fault displacement. Seismic velocity temporarily decreased even after the earthquake, suggesting that aftershocks continued to generate fault zone damage and stress state changes. Especially, the seismic velocity at the eastern and western edges of the Hinagu-Futagawa fault system, where the aftershocks were concentrated (2022), did not recover (fig. S1), and the seismic velocity in these areas continued to fluctuate several months after the mainshock (September to November 2016; movie S1). Far from the seismogenic faults and the volcano (Fig. 2E), seismic velocity slightly decreased or fluctuated after the earthquake, possibly owing to temporal formation weakening in a wide area (>50 km from the faults and volcano), such as pore pressure increase due to vibration. Although the velocity recovering process is influenced by several mechanisms, the spatiotemporal variation of seismic velocity estimated in this study could be important information for estimating fault strength recovery processes and future earthquake activity.

The largest seismic velocity reduction (0.7 to 0.8%), which occurred soon after the 2016 earthquake (Figs. 2F and 4B and movie S1), was observed in the western part of the Aso caldera, where there is a magmatic body at a depth of <10 km (26). The large velocity reductions were observed at several receiver pairs across the Aso magmatic body (fig. S1B). During the 2011 Tohoku-oki earthquake, similar seismic velocity reduction was also observed in a volcanic area far from the seismogenic fault (2). The earthquake shaking could mobilize fluids around the magmatic body, and this fluid could inflate the shallower formation and open cracks. Under high pore pressure conditions around a magmatic body, seismic velocity is sensitive to stress or pore pressure changes because thin cracks are easily generated (9, 10). The stress perturbation around the magmatic body during the 2016 earthquake was estimated to be ~1.8 MPa (25), and these stress changes are sufficient to open a preexisting crack (38). Open cracks generated around the magma chamber by the 2016 earthquake caused the seismic velocity to decrease greatly. Although it is possible that strong earthquake shaking weakened the volcanic body without fluid pressure, rapid velocity recovering in the Mount Aso area (Fig. 2F) could demonstrate that the velocity reduction is mainly associated with fluid pressure.

Seismic velocity at the Aso caldera recovered rapidly and was faster than the pre-earthquake velocity after July 2016 (Figs. 2F and 4, E and F, and movie S1). High volcanic activity after the earthquake, including eruptions on 16 April and 1 May 2016, may have caused depressurization, leading to the increased seismic velocity, a phenomenon observed at other volcanoes (12). Furthermore, seismic velocity around the Aso caldera again increased greatly after the large eruption on 7 and 8 October (Fig. 2F).

To reveal the vertical velocity change (that is, the change with depth), we calculated the temporal velocity variations in different frequency ranges for two seismometer pairs (Fig. 5): (i) one across the seismogenic Futagawa fault (A-A′ station pair) and (ii) the other across the Aso caldera (C-C′ station pair). The frequency dependence displayed by these two seismometer pairs had different characteristics. As the frequency became higher, the changes in seismic velocity across the seismogenic fault became greater (Fig. 5A), characteristics similar to those reported across the fault of the Mw 6.0 Parkfield earthquake (37). This result suggests that the damage was intensive around the fault zone from the surface to the depth of the largest fault rupture at a depth of ~5 km. In contrast, the seismic velocity variations around the Aso caldera were similar in both low-frequency (0.1 to 0.6 Hz) and high-frequency (0.4 to 0.9 Hz) ranges (Fig. 5B). A similar trend was observed in the case of the Mw 6.5 San Simeon earthquake (37). This result indicates that the depth range of the velocity change was deeper at the Aso caldera than in the fault zone. If S-wave velocity at the Aso caldera decreased only at shallower depths (<5 km), then the resulting velocity variation in the lower frequency range would be smaller than that in the higher frequency range, because Rayleigh waves at low frequencies have low sensitivity at shallower depth (Fig. 3). Therefore, S-wave velocity at depths of about 3 to 10 km, where lower frequency Rayleigh waves have larger sensitivity (Fig. 3), should decrease. Because this depth range roughly agrees with the location of the magma chamber (26), the seismic velocity may have been reduced around the magma chamber during the earthquake. Another possible cause of the velocity variation in the volcanic area is rainfall (31). However, the seismic velocity indicated by the higher frequency component, which reflects velocity at shallow depth, in the Mount Aso area did not show any clear effect that could be attributed to precipitation (C-C′ station pair; Fig. 5B). We therefore concluded that the seismic velocity variation revealed in this study (Figs. 2 and 4) was not greatly influenced by rainfall.

Fig. 5 Frequency dependence of temporal variations of seismic velocity between two typical seismometer pairs.

(A) Across the fault plane from A to A′ in Fig. 1. (B) Across Mount Aso from C to C′ in Fig. 1. Curves show daily variations of the estimated velocity change Embedded Image with respect to changes before the Kumamoto earthquake. Vertical dashed lines show the time window (30 days) influenced by the mainshock. The frequency-dependent velocity variation is likely associated with surface wave dispersion (37), supporting our assumption that the coda waves were dominated by surface waves.

Here, we report mapping results showing the continuous temporal variation in seismic velocity across the fault and volcano during the Kumamoto earthquake. Past studies have used similar approaches for velocity estimation, but higher spatial resolution of our estimated seismic velocity variation in a broad area allowed us to identify the spatial distribution of the damage zone or stress state. This approach can be applied anywhere that permanent seismometers are deployed, such as Hi-net in the Japanese islands (19) and arrays in the United States (39). Velocity changes due to slow slip or foreshocks identified through this crustal-scale monitoring could thus be used in the estimation of future earthquakes. The denser deployment of seismometers around fault zones and volcanoes (17) would allow local anomalies related to earthquake and volcanic activities to be accurately resolved.


Computing daily cross-correlations

We divided 1 day of continuous ambient noise data into 30-min segments with a 50% overlap (40). We then corrected the instrumental response and applied a 0.1- to 0.9-Hz band-pass filter to each segment. To remove earthquake data, we rejected any segment in which the maximum absolute value exceeded a specified threshold (40, 41). We computed the mean and SD of the absolute value in each window. We applied four threshold values: 5, 7, 9, and 11 times the SD of the mean. To keep a certain number of segments in each day, we gradually increased the threshold value until more than half of the segments satisfied the threshold. The daily cross-correlations were computed by stacking each segment of power-normalized cross-correlations (cross-coherence) between receiver pairs in the frequency domain (42, 43).

Estimating seismic velocity change by stretching interpolation

To estimate seismic velocity changes, we used a stretching interpolation technique (3, 29, 44). The method elongates the time axis and finds the most similar trace to the reference traceEmbedded Image(1)Embedded Image(2)where ε is a stretching parameter, fref represents the reference trace, fcur represents the current trace, and C(ε) is the correlation coefficient between the reference and the current trace. The parameter ε corresponds to a relative time shift Embedded Image and relates to a velocity change Embedded Image as followsEmbedded Image(3)

We estimated the relative velocity change ϵ that maximizes C(ε) by adopting a grid search algorithm and a ternary search algorithm. A positive value of ε indicates a decrease in velocity compared to the reference trace.

In our study, the stretching interpolation technique was applied to the coda of cross-correlations because the coda is more sensitive to velocity change than the direct wave (29) and less sensitive to variations in noise sources (6). The reference trace was defined as the coda of cross-correlations of a 1-year stack, and the current trace was defined as the coda of cross-correlations of a 30-day stack. Following Meier et al. (29), we selected a 100-s window for the coda. We calculated C(ε) for the range −0.025 < ε < 0.025 with a step size of 0.0005 and picked ε with the maximum value of C(ε). By applying the ternary search method, we estimated the seismic velocity change ε with a resolution of ~5 × 10−7. We further divided the 100-s window into six smaller time windows of 50 s each and applied the stretching interpolation technique in each window to compute SDs of estimates of ε (29). The measured velocity change was considered to represent the velocity change in the middle of the 30-day window. By moving the 30-day window used for the current trace, we estimated the daily variation of the velocity change.

Cross-correlations have causal (positive-time) and anticausal (negative-time) sides, corresponding to a wave traveling from station D to station D′ and from station D′ to station D (for example, fig. S3). Therefore, the stretching interpolation technique can be applied to both positive- or negative-time parts. However, we used only the positive- or negative-time part to reduce the effects of strong volcanic tremors from Aso volcano (45, 46), as described later.

The estimated velocity change is the relative velocity variation of the reference trace, including ambient noise before and after the 2016 Kumamoto earthquake. To measure the velocity change caused by the earthquake, we subtracted averaged velocity variations before the earthquake from the estimated velocity change ε as followsEmbedded Image(4)

We note that the averaged velocitychange ε0 was obtained by averaging the estimates of velocity changes before the earthquake (from 15 December 2015 to 29 March 2016). In Figs. 2 (D to F) and 4, we display the relative velocity variation Embedded Image to clarify the velocity changes associated with the 2016 earthquake.

Estimating seismic velocity change by MWCS analysis

We also tested an MWCS analysis (1, 11, 30, 31) to estimate seismic velocity change (fig. S4). Similar to the stretching interpolation, we computed temporal velocity change between the coda of cross-correlations of a 1-year stack and the coda of cross-correlations of a 30-day stack. The frequency range in this analysis was 0.1 to 0.9 Hz. We applied MWCS analysis for the time window between the starting time of the window for the stretching interpolation analysis and a lapse time of 60 s (2). In our analysis, we use a moving time window of 10 s with an overlap of 80% (2, 6). We basically followed the MWCS processing procedures of Clarke et al. (30). However, we did not apply phase unwrapping of the cross-spectrum between the current and reference traces in each moving time window, because it was difficult to obtain stable results for several receiver pairs after applying phase unwrapping. Therefore, the maximum velocity change we could measure in this approach was limited as lapse time and frequency increase. For example, at a lapse time of 60 s and a frequency of 0.80 Hz, the maximum velocity change is ~1%. This value is mostly within the range of estimated velocity variation by the stretching interpolation (fig. S4).

Influence of volcanic tremors from Aso volcano

The ambient noise data used in this study include persistent volcanic tremors from Aso volcano (45, 46). In the cross-correlations, we observed strong signals associated with volcanic tremors from Aso volcano (for example, see fig. S3). The direct wave from Aso was detected as a low-frequency signal. Because the tremors at Aso were spatially localized, they were not symmetric in the cross-correlations. However, even if the noise is localized, the effect of noise directivity can be reduced by using the coda wave (47).

The disadvantage of the stretching interpolation technique is that it has the potential to induce an apparent velocity change if the frequency component of the ambient noise varies (48). Therefore, if the frequency component of the tremors at Aso varied, then an apparent velocity change would be introduced in the estimated velocity change. In our analysis, to reduce the effects of the tremors at Aso, we applied the stretching interpolation technique to either the causal or the anticausal part of the coda, choosing that in which the signals of the tremors at Aso were weak. For this purpose, we first estimated the domain (causal or anticausal) of the direct waves from Aso based on the geometry among Aso volcano and the two stations used for the cross-correlations, and then selected the opposite window for the stretching analysis (for example, positive-time in fig. S3).

The cross-correlation coefficient in the stretching analysis across the Aso volcano is lower after the large eruption on 7 and 8 October (Fig. 2F). The lower coefficient in the stretching analysis demonstrates that the ambient noise characteristics around Aso, such as source locations and frequency components of noise, were changed by the large eruption. Because the seismometer pair (C-C′) crosses the Aso noise source, we could not efficiently remove the influence of ambient noise variation at Aso volcano.

On the other hand, temporal change in seismic velocity measured by the MWCS is likely less affected by temporal change in the amplitude spectrum of ambient noise because only phase information are used in velocity variation measurements (48). In the MWCS, however, we need to assume a constant time shift in each moving time window. Although longer time window in the MWCS reduces the fluctuation in time shift measurements due to noise, the assumption of a constant time shift is getting more erroneous when the window length is increased (49).

To evaluate the reliability of measurements of temporal velocity change by the stretching interpolation and the MWCS for our data, we compared temporal velocity variations of several receiver pairs using both techniques (fig. S4). The velocity change derived from the MWCS is similar to that from the stretching interpolation between A and A′ (fig. S4A). For this receiver pair, we observed high stability in velocity change estimation during linear regression in the MWCS (fig. S4D). However, the velocity changes derived by the MWCS are different from those by the stretching interpolation for the receiver pairs of B-B′ and C-C′ (fig. S4, B and C). Between B and B′, the measured time shifts are not well fitted to the straight line in the negative-lag time where direct waves from Aso are included (fig. S4E). Between C and C′, the measured time shifts do not show a clear linear trend in positive- and negative-lag times, indicating the difficulty of linear regression to measure velocity change (fig. S4F). This probably indicates that both positive- and negative-lag times are affected by the tremors at Aso, because Aso volcano is located between C and C′ (Fig. 1).

To stabilize velocity change estimation in the MWCS, we usually use both positive- and negative-lag times. However, because virtual seismograms for either positive- or negative-lag time or both around Aso volcano are influenced by tremors at Aso (fig. S3), using both lag times decreased the stability of velocity estimation (that is, lack of a clear linear trend in fig. S4, E and F). On the other hand, we could apply the stretching interpolation by selecting either positive- or negative-lag time where the influence of tremors at Aso was small. Thus, we obtained more stable results using the stretching interpolation in the central Kyushu Island.

Proxy of depth resolution

To evaluate the depth range associated with the estimated surface wave velocity changes (Figs. 2 and 4), we calculated depth-dependent sensitivity kernels for fundamental mode Rayleigh waves with respect to S-wave velocity (Fig. 3) (50). The sensitivity kernels were computed by the DISPER 80 program (51) for a 1D layered model near Mount Aso estimated by Nishida et al. (52). A similar evaluation method is described by Chaves and Schwartz (6). Note that the sensitivity kernels for higher frequencies might not be accurate because Nishida et al. (52) analyzed Love and Rayleigh waves below 0.2 Hz and shallow sediments were not included in their model.

Mapping changes in seismic velocity

We obtained seismic velocity changes between pairs of seismic stations (fig. S2). For each station, we estimated velocity changes between that station and other stations within a distance of <40 km. The averaged velocity changes were assigned at each station as in the study by Brenguier et al. (2). We then used a linear interpolation technique to obtain a map of the velocity changes (see Fig. 4 and movie S1).


Supplementary material for this article is available at

fig. S1. Temporal variation of seismic velocity between each receiver pair in relation to earthquake activity.

fig. S2. Spatial and temporal variation of seismic velocity in central Kyushu during the 2016 Kumamoto earthquake.

fig. S3. Asymmetry of cross-correlation between station D and D′ (1-year stack) due to tremors from Aso volcano.

fig. S4. Comparison of temporal changes of seismic velocity derived from the stretching interpolation and the MWCS.

movie S1. Animation showing spatial and temporal variations of seismic velocity on Kyushu Island during the 2016 Kumamoto earthquake.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We are grateful to the National Research Institute for Earth Science and Disaster Prevention (NIED) for providing us with the Hi-net data. Funding: This study was supported by the Japan Society for the Promotion of Science (JSPS) through Grants-in-Aid for Scientific Research in Innovative Areas (JP15H01143 and JP17H05318) and a Grant-in-Aid for Young Scientists B (JP16K18332), and by the Ministry of Education, Culture, Sports, Science and Technology of Japan (MEXT) through the World Premier International Research Center Initiative. Author contributions: All authors conducted the data analysis. T.I. and T.T. conceived the approach. T.T. interpreted the results and wrote the draft. All authors approved the draft. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Seismic data required to evaluate the conclusions in the paper are available from NIED. The earthquake catalog used in Fig. 1 and fig. S1 was produced by the Japan Meteorological Agency in cooperation with MEXT. The coastline used in figs. S1 and S2 and movie S1 is based on the National Land Numerical Information Coastal Lines data,, provided by the Ministry of Land, Infrastructure, Transportation and Tourism of Japan.

Stay Connected to Science Advances

Navigate This Article