FDM simulation of earthquakes off western Kyushu, Japan, using a land–ocean unified 3D structure model

Abstract Seismic activity occurred off western Kyushu, Japan, at the northern end of the Okinawa Trough on May 6, 2016 (14:11 JST), 22 days after the onset of the 2016 Kumamoto earthquake sequence. The area is adjacent to the Beppu–Shimabara graben where the 2016 Kumamoto earthquake sequence occurred. In the area off western Kyushu, a M7.1 earthquake also occurred on November 14, 2015 (5:51 JST), and a tsunami with a height of 0.3 m was observed. In order to better understand these seismic activity and tsunamis, it is necessary to study the sources of, and strong motions due to, earthquakes in the area off western Kyushu. For such studies, validation of synthetic waveforms is important because of the presence of the oceanic water layer and thick sediments in the source area. We show the validation results for synthetic waveforms through nonlinear inversion analyses of small earthquakes (~M5). We use a land–ocean unified 3D structure model, 3D HOT finite-difference method (“HOT” stands for Heterogeneity, Ocean layer and Topography) and a multi-graphic processing unit (GPU) acceleration to simulate the wave propagations. We estimate the first-motion augmented moment tensor (FAMT) solution based on both the long-period surface waves and short-period body waves. The FAMT solutions systematically shift landward by about 13 km, on average, from the epicenters determined by the Japan Meteorological Agency. The synthetics provide good reproductions of the observed full waveforms with periods of 10 s or longer. On the other hand, for waveforms with shorter periods (down to 4 s), the later surface waves are not reproduced well, while the first parts of the waveforms (comprising P- and S-waves) are reproduced to some extent. These results indicate that the current 3D structure model around Kyushu is effective for generating full waveforms, including surface waves with periods of about 10 s or longer. Based on these findings, we analyze the 2015 M7.1 event using the cross-correlations between the observed and synthetic waveforms. The result suggests a rupture propagation toward the NNE, with a major radiation about 25 km north of the onset point. Graphical Abstract .

section_content: Japan is a country surrounded by seas where large earthquakes often occur due to the tectonic settings, such as the subductions of the Pacific and Philippine Sea plates. Therefore, it is important to construct three-dimensional (3D) structure models in the oceanic area and validate the synthetic seismic waveforms from suboceanic earthquakes, to study the detailed rupture processes and evaluate strong ground motions (e. g., Okamoto 2002; Nakamura et al. 2014 Nakamura et al., 2015)). The seismic activity off western Kyushu, Japan, that began with an M7. 1 earthquake on November 14, 2015 (JST) (Figs. 1, 2 ) provides a unique opportunity to test the current 3D structure model, since the M7. 1 event was the largest recorded by the Japan Meteorological Agency (JMA) in the region over the last 93 years, and many aftershocks occurred throughout a broad area (Fig. 1 ). The M7. 1 event also generated a tsunami with a height of 0. 3 m that was observed at Nakanoshima Island. More seismic activity occurred on May 6, 2016 (JST), 22 days after the onset of the 2016 Kumamoto earthquake sequence (Fig. 2f ). 

The area is located in the northern part of a back-arc basin, the Okinawa Trough, and lies southwest of the Beppu-Shimabara graben (BSG), where the 2016 Kumamoto earthquake sequence occurred. Both the Okinawa Trough and the BSG are spreading actively, with rates of 10 mm/year for the northern part of the Okinawa Trough (Nishimura et al. 2004 ) and 14 mm/year for the BSG (Tada 1985). Based on these tectonic and geophysical backgrounds, it has been argued that the BSG is possibly the landward extension of the Okinawa Trough (Tada 1984 (Tada, 1985;; Takayama and Yoshida 2007). We note here that the seismicity in that region, shown in detail in Fig. 2, also supports this view. The 24-h aftershocks of the 2015 M7. 1 event have a linear distribution with a length of about 70 km, which can be considered to be approximately equal to the length of the fault (Fig. 2b ). The seismic activity then spread toward the northeast, namely toward the BSG (Fig. 2c ). The northern aftershocks did not follow in an extension of the line formed by the 24-h aftershocks, but formed a bend toward the BSG in the final distribution (Fig. 2d ), implying some connection between the northern Okinawa Trough and the BSG. 

In order to improve our understanding of these seismic activities, including the large (M ~ 7) earthquake, tsunami generation, and tectonic settings, it is necessary to study the earthquake sources and the strong motions due to the earthquakes off western Kyushu. For this purpose, we need to analyze waveform data for the earthquake source parameters using synthetic Green tensor waveforms. This analysis, however, also requires the validation of the synthetic waveforms themselves due to the presence of an oceanic water layer and thick sediments that can cause large effects on the excitations and propagations of the seismic waves, in the source area. 

Therefore, in this study, we validate the synthetic waveforms computed for the current 3D structure model and study the structural effects on the estimated earthquake parameters in the areas off western Kyushu. For these purposes, we select small earthquakes (~M5) that can be approximated with point sources, generate synthetics
section_name: Introduction

section_content: We analyze five events that occurred in the area off western Kyushu (Table 1 ; Fig. 3 ), as the ground motions from these events are recorded with a sufficient signalto-noise ratio. We use waveform data observed at eight K-NET stations and one F-net station operated by National Research Institute for Earth Science and Disaster Resilience (NIED), as shown in Fig. 3a ; both records are processed and filtered to yield the ground velocity components.
section_name: Data

section_name: Methods

section_content: As described above, first, we estimate "best" point source parameters in order to avoid biases in the synthetic waveforms. We apply a nonlinear inversion method (Okamoto and Takenaka 2009) to the moderate-sized events, nos. 2-5 (referred to as EV2 to EV5 hereafter), to determine the moment tensor and the source time function of a point source simultaneously. The source time function is represented by a series of unit triangles with widths of 2 s. The synthetic waveform S i (x, t) for the i-th station located at x and at time t is where M ℓ (ℓ = 1,. . ., 5) denotes the five elementary basic moment tensors (Kikuchi and Kanamori 1991), G ℓ i (x, t; x 0, t 0 ) are the corresponding Green tensor wave- forms computed for a point source with the unit pulse placed at x 0 with origin time t 0, �τ is the origin time (1)

correction, A k is the amplitude of the k-th unit pulse, and �ξ = 1 s. The sixth (isotropic) elementary moment tensor is not used. In order to avoid spurious oscillations in the source time function, nonnegative conditions A k = α 2 k (k = 1,. . ., N k ) are imposed, and α k values are used as the inversion parameters. The nonlinear part of the inversion is solved using an iterative algorithm (Marquardt 1963). As in Okamoto and Takenaka (2009), we minimize the square residual F, defined below, in the inversion: where D i (x i, t) and S i (x i, t) are the i-th observed and synthetic waveforms, respectively, T i is the length of the dataset, N W is the number of waveforms, T = N W i=1 T i, and w i is weight, which is fixed as unity in this analysis. Note that time shifts between data and synthetics are allowed only through the "origin time correction term (�τ)" in Eq. 1. Next, we perform a space-time grid search to infer the best point source position x 0 and origin time correction �τ that minimizes the square residual (Eq. 2). For additional details, see "Inversion Procedure" section in Additional file 1.
section_name: Waveform inversion and space-time grid search

section_content: In the inversion analysis, we apply two different passbands to the observed waveforms: One is for "shortperiod" data, with a passband of 4-40 s, and the other is for "long-period" data with a passband of 10-40 s. As discussed further in this paper, for shallow earthquakes in the oceanic area, long-period waveforms (longer than about 10 s) can be reproduced well compared to shortperiod waveforms (less than about 10 s). However, for shallow earthquakes, the surface waves dominate the long-period components, and important information related to the body waves (P-S times, for example) that are required to deal with the trade-off between the origin time and the source location can be obscured. Therefore, we apply two passbands to the same raw data in order to incorporate information based on both the body waves and surface waves. We call the resultant moment tensor solution the first-motion augmented moment tensor (FAMT) because the first motions of the P-and S-waves are distinct in the short-period waveforms.
section_name: First-motion augmented moment tensor (FAMT)

section_content: We construct a land-ocean unified 3D structure model because we should consider the effects of the oceanic water layer and the thick sediments (Fig. 3b ) when simulating the seismic wave propagation from the shallow suboceanic earthquakes, as they strongly affect the excitation and propagation of seismic waves, especially surface waves. We compile the land-ocean topography (Kisimoto 2000) and the Japan Seismic Hazard Information Station V2 model (Fujiwara et al. 2012) for the subsurface structure, as well as the Japan integrated velocity structure model (JIVSM; The Headquarters for Earthquake Research Promotion 2012) for the depths of the Conrad and Moho surfaces and the velocity structure in the crust. We extrapolate the Conrad and Moho surfaces of the JIVSM toward the west (down to 128°E) based on the results of the seismic survey in the nearby area (Iwasaki et al. 1990 ) because the JIVSM is defined longitudinally from 129°E in this region. An example of a vertical cross section of the 3D model is shown in Fig. 3b ; the oceanic water layer (indicated by blue color) and land topography are incorporated in the model.
section_name: Land-ocean unified 3D structure model

section_content: The Green tensor waveforms are computed using a reciprocal method (Okamoto 2002 ): A single force is applied to the station location, and the response strains in the source region are stored as the components of the reciprocal Green tensor. We apply a 3D HOT-FDM scheme (Nakamura et al. 2012 ) that implements the correct fluid-solid boundary scheme (Okamoto and Takenaka 2005) to simulate wave propagation in the land-ocean unified 3D model. We use a program that incorporates a multi-graphic processing unit (GPU) acceleration (Okamoto et al. 2010 (Okamoto et al., 2013) ) and conduct the simulation on the TSUBAME 2. 5 supercomputer at the Tokyo Institute of Technology, Japan. Anelastic attenuation is introduced using the τ-method (Blanch et al. 1995) with a correction term (Carcione 2001). The grid and time intervals are 100 m and 0. 005 s, respectively; based on the shortest wavelength of 6 × x = 600 m (Moczo et al. 2000) and the minimum S-wave velocity of 650 m/s in the 3D model, we define the maximum frequency as 1. 1 Hz. We use six relaxation mechanisms to introduce the viscoelasticity with approximately constant Q in a frequency band from 0. 01 to 1. 1 Hz. The FDM grid size is 2880 × 3520 × 930 (NS × EW × depth; Fig. 3a ), and the subdomain size assigned to a GPU is 320 × 320 × 310. Using 297 GPUs (i. e., 99 nodes of the TSUBAME 2. 5), 10,850 s (including times for I/O) is required on average to perform a simulation for 24,000 time steps. For some stations, because of limited computational resources, only the Green tensor waveforms that compose the vertical component waveform at the station location are computed. Note that as per the reciprocal method, a single full 3D simulation is required to compute the Green tensor waveforms that compose a single component of the displacement (or the ground velocity) at a single station location.
section_name: FDM simulation

section_name: Results and discussion

section_content: First, we show the results of the grid search analysis of EV2. For the origin time correction, the minimum residual is found for a correction time of 1 s (Additional file 1: Figure S1 ). Figure 4a shows the horizontal slices of the residual distribution of the space grid search for the best origin time correction of 1 s. Although there are some (relatively weak) local minimums, the global minimum is well defined at a depth of 8. 1 km and about 8 km landward from the JMA-specified epicenter. In Fig. 4b, we compare the observed waveforms with the synthetics computed for the best point source parameters. The first 13 pairs from KGS003U to KYK U are the short-period waveforms that are added to estimate the FAMT solution. 

The synthetic waveforms reproduce the observed waveforms well, especially for long-period components. Even for the short-period components, the peaks and troughs in the waveforms (up to the first swing of S-wave) are well reproduced at stations KGS003, KGS004, KGS036, and KYK, which are distant from the P-wave nodes. In order to show the differences between the 1D and 3D structure models, we compute synthetic waveforms using the 1D model for the Kyushu region (e. g., Takenaka et al. 2006 ; Additional file 1: Figure S2 ). The earthquake parameters for 1D waveforms are the same as those of the FAMT for EV2 (see "1D structure model" in Additional file 1 for details). Typically, the amplitudes and the phases of the large later arrivals (surface waves) are not reproduced well by the 1D model (Fig. 4c for five selected waveforms and Additional file 1: Figure S3 for all waveforms). Also note that the residual (F) for 1D waveforms (1. 54) is larger than that for 3D waveforms (0. 55).
section_name: Grid search analysis of EV2

section_content: To determine the degree to which we can reproduce the characteristics in the observed waveforms, we compare bandpass-filtered waveforms using different passbands for EV2. Note that at some stations, the triggering times were late and the first motions were not recorded. Thus, in the comparisons we use 13 components that include the P-wave first motions, of which two typical components are displayed in Fig. 5. All 13 components are Here we plot two typical components as an example (for all components see Additional file 1: Figure S4 ). FAMT parameters are used to generate the synthetic waveforms. The passbands and the normalized residuals (F in Eq. 2) for all the thirteen components used in this analysis are attached to each trace. The symbols "P" and "S" at the top traces denote the onset of P-and S-wave packets, respectively -18. 2, 187. 2 86. 0, 198. 2, 282. 8 shown in Additional file 1: Figure S4. As a measure of the fit of the synthetics to the data, we also indicate the total residuals (F) computed using the 13 components (Fig. 5 ). The source parameters determined by the FAMT analysis are used (Table 2 ). For the long-period passband of 10-15 s, the synthetics reproduce the observed full waveforms well, both in phase and in amplitude (Fig. 5 and Figure S4(e) ). As noted before, even for shorter periods of 4-6 and 6-10 s, the first body wave (up to the first swing of the S-wave) is reproduced to some extent at stations far from the P-wave nodes (KGS003, KGS004, KGS036, KYK: third and fourth traces in Fig. 5 and Figure S4(c ) and (d)). However, even at these stations, for passbands of 4-6 and 6-10 s, the later surface waves in the observed and synthetic traces could be out of phase with each other and/or the amplitudes of the synthetic waveforms could deviate from those of the observations. Thus, at these passbands, the residuals computed for all 13 components increase. For shorter passbands of 2-4 and 1-2 s, even the first part of the waveforms is difficult to reproduce. These results indicate that the current 3D structure model around Kyushu is effective for generating full waveforms, including surface waves with periods about 10 s or longer. For shorter periods, down to around 4 s, the first part of the waveforms (i. e., P-and S-waves) is reproduced well to some extent and can be used for waveform analysis. Figure 5 (bottom) and Figure S4 (f ) show broadband (1-40 s) waveforms. The residual for the broadband case is larger than that for the case of the 10-15 s band, but smaller than those of other cases with shorter-period bands. This implies that the longperiod components are slightly larger in the EV2 strong motion records, which would be a favorable situation for 3D modeling.
section_name: Validating waveforms

section_content: Figure 6 displays the results of the grid search analyses for EV3 to EV5. In Table 2, we summarize the estimated source parameters. The source time functions are shown in Figs. 4 and 6. The nodes of EV3 and EV4 are oriented in NNE-SSW directions (≈200°N to 213°N), which are slightly different from the near north-south trend of the node of EV2. The difference suggests slightly different fault orientations between the aftershock area of EV1 (Fig. 2b ) and the area of 2016 seismic activity (Fig. 2f ). The estimated FAMT locations shift landward from the JMA-specified epicenters by about 13 km on average (Fig. 7 ). The spatial shifts are not likely the result of the difference between the hypocenter and the centroid because the shifts are systematic and large considering the size of the events (~M5) and because we incorporate the first-motion phases in the analysis. Thus, we consider the shift to reflect the difference in the models used in the hypocenter analysis and the 3D structure used in this study. That is, some biases may be incorporated in the earthquake parameters if a simple 1D structure is used in the analysis of earthquakes in oceanic regions. 

We note here that the FAMT solutions obtained in this study have some common features with the F-net moment tensor solutions: Both solutions have axes trending northwest-southeast to north-south (Additional file 1: Figure S5 ). In addition to the moment tensors, the FAMT analysis determines the source locations and the source time functions. On the other hand, the lateral locations are fixed to the JMA epicenters, and the source time functions are not estimated in the F-net solutions. These parameters (detailed locations and source time functions) can be obtained using short-period waveforms, which requires 3D waveform modeling. This feature is the advantage of the FAMT analysis proposed in this paper.
section_name: FAMT solutions of EV3 to EV5

section_content: Using the Green tensor waveforms computed and validated in this study, we analyzed the source of the M7. 1 earthquake (EV1). Since there could be biases in the aftershock parameters that we refer to in constructing fault model (Figs. 7, 8a), as a preliminary step, we apply a cross-correlation analysis that does not require an inversion procedure or detailed constraints on the rupture model. We use a rectangular grid for the fault and compute the following averaged cross-correlation R ij (τ) between the observed and synthetic waveforms, where the (i, j) pair denotes a grid point on the fault grid, N W is the number of components, D (k) (k) is set to 120 s after the origin or the end time of the observed waveform. Fifteen long-period components (passband of 10-40 s) from all nine stations in Fig. 3a are used in this study. The base width of the triangular source time function used to generate each S ij (k) (t) is determined by a parameter study in which we employ the correlation analyses for base widths from 1 to 10 s with an interval of 1 s. We select the base width of 6 s, which provides the highest correlation coefficient. 

Referring to the best double couple of the GCMT solution and aftershock distribution (Fig. 8a ), we select a near north-south trending, approximately vertical nodal plane (dip = 86, slip = -165, strike = 192), and approximate it with a vertical rectangular grid whose grid spacings are 5. 1 and 3 km in the horizontal and vertical directions, respectively (15 × 6 grid points in total). We place the fault grid points by referring to the hypocenter specified by JMA. Note that the locations of the grid points deviate slightly from the aftershock region if we align the grid points based on the strike of the nodal plane (Fig. 8a ). The hypocenters of both the M7. 1 event and the aftershocks must be refined in the future analysis considering the biases in the aftershock locations that we have shown in this paper. The moment tensor of the GCMT solution is used to generate S ij (k) (t). We plot the distribution of R ij (τ) for different values of τ in Fig. 8b : The origin (range = 0 km, depth = 0 km) is placed at the JMA epicenter, and the range is measured toward the NNE. At the lag time τ = 0 s, the reddish area with high cross-correlation values is near the JMA hypocenter (0, 17 km): The maximum value (0. 39) occurs at a grid point with a range of 5. 1 km and a depth of 18. 1 km (Fig. 8b, top). Then, the reddish area of high cross-correlation values moves NNE over time: The "grand" maximum value is obtained for τ = 4 s at a grid point with a range of 25. 5 km and a depth of 15. 1 km (Fig. 8b, middle ). These distributions of cross-correlation roughly indicate a "smeared" space-time extent of the major radiation that developed the main (large) surface wave packets in the observed waveforms rather than the rupture front of the propagating source; this is because we use components with passband of 10-40 s, whose lower (shorter) bound (10 s) is much longer than the period of around 1 s that is usually used for backprojection analyses of the rupture front. Although in definition the cross-correlation value does not reflect the absolute amplitude, we still regard high values as roughly corresponding to major radiations because high values are obtained when the synthetic waveforms fit well with the observed large-amplitude surface waves (Fig. 8c ). Therefore, we think that the rupture initiated at the southern part of the fault near the JMA hypocenter and propagated NNE, with major radiation 25 km NNE of the hypocenter. Note that the GCMT solution about 22 km NNE of the JMA epicenter and at depth of 12 km is near the point of the maximum cross-correlation shown by a cross in Fig. 8b, middle. The depth of the maximum correlation (15. 1 km) is between the GCMT depth (12 km) and the United States Geological Survey body-wave moment tensor (17. 0 km). The aftershocks also distribute down to deep parts of the assumed fault (Additional file 1: Figure S6 ). We also note that the reddish area of large correlation extends from deep to shallow parts of the fault (Fig. 8b, middle), which implies a complex, distributed radiation.
section_name: The 2015 M7.1 earthquake off western Kyushu (EV1)

section_content: We validated synthetic waveforms using a nonlinear waveform inversion analysis of small earthquakes (~M5) off western Kyushu. We used a land-ocean unified 3D structure model, 3D HOT-FDM (Nakamura et al. 2012), and a multi-GPU acceleration (Okamoto et al. 2010 (Okamoto et al., 2013) ) to simulate wave propagations. We estimated the first-motion augmented moment tensor (FAMT) solution based on both long-period surface waves and shortperiod body waves. The synthetics computed for the FAMT solutions reproduce well the observed waveforms with periods of 10 s or longer. However, for shorter periods the later surface waves are not reproduced well, while the first part (P-and S-waves) is reproduced well to some extent. These results indicate that the current 3D structure model around Kyushu is effective for generating full waveforms, including surface waves with periods about 10 s or longer. For shorter periods, the first part of the waveforms (i. e., P-and S-waves) can be used for waveform analysis. We also found that the FAMT solutions systematically shift landward by about 13 km on average from the JMA epicenters. Based on these findings, we analyzed the 2015 M7. 1 event using the cross-correlations between the observed and synthetic waveforms. In this analysis, we used long-period waveforms with passband of 10-40 s because we need to use full waveforms. The result suggests a rupture propagation toward the NNE with a major radiation about 25 km north of the onset point.
section_name: Conclusions

section_name: Acknowledgements

section_content: The authors are grateful to the JMA and GCMT project for providing the earthquake parameters; to Kiyoyuki Kisimoto, NIED, and the Headquarters for Earthquake Research Promotion for providing the structure models; to Minoru Takeo for providing the 1D synthetic seismogram program; and to NIED for providing the waveform records of K-NET and F-net. Comments by Kimiyuki Asano (editor) and two anonymous reviewers were helpful in improving the manuscript. This study is partially supported by KAKENHI ( 26282105 ) and JHPCN (jh160029-NAH).
section_name: Acknowledgements

section_name: Abbreviations

section_content: FDM: finite difference method; 3D: three-dimensional; JMA: Japan Meteorological Agency; FAMT: first-motion augmented moment tensor; NIED: National Research Institute for Earth Science and Disaster Resilience; GCMT:
section_name: Abbreviations

section_content: Additional file 1. Additional description of the waveform inversion procedure in FAMT analysis and 1D structure model used in the section "Grid search analysis of EV2" and in Fig. 4. Figure S1. Residuals (F in Eq. 2) versus origin-time correction for EV2 (Table 1 and Fig. 3a ). We employ grid searches by assigning different values of origin-time correction and plot the minimum residuals in each grid search. The "grand" minimum value is obtained when we add 1 s (i. e, a delay of 1 s) to the JMA origin time.
section_name: Additional file

section_content: TO compiled the structure models to construct the 3D structure model; coded the FDM program; performed the FDM simulation, inversion analysis, and the cross-correlation analysis; and wrote the draft of the manuscript. HT and TN took part in discussions of the 3D structure model and the interpretation of the results of the inversion analyses. TH proposed the analysis method using waveform cross-correlation. All authors contributed to the discussion and conclusions. All authors read and approved the final manuscript.
section_name: Authors' contributions

section_content: 1 Department of Earth and Planetary Sciences, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro, Tokyo 152-8551, Japan. 2 Department of Earth Sciences, Okayama University, 3-1-1 Tsushima-Naka, Kita-ku, Okayama 700-8530, Japan. 3 Earthquake and Tsunami Research Division, National Research Institute for Earth Science and Disaster Resilience, 3-1, Tennodai, Tsukuba, Ibaraki 305-0006, Japan. 4 International Institute of Seismology and Earthquake Engineering, Building Research Institute, 1 Tachihara, Tsukuba, Ibaraki 305-0802, Japan.
section_name: Author details

section_content: The authors declare that they have no competing interests.
section_name: Competing interests

section_content: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
section_name: Publisher's Note