Subsurface velocity structure and site amplification characteristics in Mashiki Town, Kumamoto Prefecture, Japan, inferred from microtremor and aftershock recordings of the 2016 Kumamoto earthquakes

In order to investigate the seismic velocity structure in the region of concentrated severe damage during the 2016 Kumamoto earthquake sequence, we conducted continuous seismic observations in the central area of Mashiki Town, Kumamoto Prefecture. During 4 days of observations at eight temporary sites, 2 months after the mainshock, recordings from 30 aftershocks (1.7 ≤ Mj ≤ 4.3, 1.9 km ≤ depth ≤ 13.5 km) were obtained. The aftershock data showed that site amplifications at approximately 1 Hz are dominant in a zone where almost no buildings were damaged along the Akitsu riverside, whereas site amplifications at higher than 3 Hz are observed in the heavily damaged zones. Our data also showed that the peak acceleration and velocity amplitudes, as well as seismic intensities for the small events in the less damaged zone, are clearly larger than those in the damaged zones, implying that the damage distribution is inconsistent with site response based on linear site amplifications. The estimated phase velocities of Rayleigh waves using the aftershock and microtremor data show dispersive characteristics in the lower frequency range (0.26 ≤ f ≤ 1.27 Hz), but the values are substantially smaller than those derived from the P–S logging model at the nearest KiK-net strong-motion observation station KMMH16. The derived microtremor horizontal-to-vertical spectral ratios and earthquake radial-to-vertical (R/V) spectral ratios show common distinct peaks at around 0.4 Hz, which are possibly related to the response of deep sedimentary layers beneath the area. The refined velocity structure model that better accounts for both the phase velocity and common dominant peak indicates that the values of S wave velocity (Vs) above the bedrock layer (Vs = 2700 m/s) are smaller than those inferred from the logging model and the depth to the bedrock layer could be much deeper (about 600 m) in comparison with the logging model (234 m). The derived R/V spectral ratio at station KMMH16 also shows a distinct peak at 0.4 Hz, suggesting that there is no large difference of deep sedimentary structure between the observation area and station KMMH16.


Introduction
The 2016 Kumamoto earthquake sequence including the largest foreshock (M W 6.1 on April 14) and the mainshock (M W 7.1 on April 16) occurred in the central part of Kumamoto Prefecture, southwestern Japan, and caused severe damage in areas between Kumamoto City and Minamiaso Village near the Futagawa and Hinagu fault zones (see the locations in Fig. 1a). Significant surface ruptures associated with the earthquake sequence appeared in Mashiki Town, which is located at the junction of the two major fault zones (Shirahama et al. 2016;Sugito et al. 2016;Toda et al. 2016) and seismic intensity of 7, the maximum seismic intensity scale of the Japan Metrological Agency (JMA), was measured during both larger events. The results of damage investigations (e.g., National Institute for Land and Infrastructure Management and Building Research Institute 2016; Sugino et al. 2016b;Tomozawa et al. 2017;Yamada 2017;Yamada et al. 2017a, b) showed that heavily damaged buildings were concentrated in a narrow region between Kumamoto Prefectural Route 28 (R28) and the Akitsu River running from east to west in the central part of Mashiki (Fig. 1b).
Field surveys after the mainshock also showed that the cause of the varying damage distribution could be explained by different site conditions of this area. According to the landform classification map for flood control published by the Geospatial Information Authority of Japan, the distribution of damaged buildings corresponds well to the river terrace between R28 and the Akitsu River, whereas almost no buildings collapsed in the floodplain and old river channel along the Akitsu riverside (e.g., Yamada 2017;Yamada et al. 2017a, c). These findings indicate that the building damage pattern is strongly affected by the site conditions associated with local subsurface geology. The results of single-point  (M j ) and the color corresponds to the hypocentral depth. Two stars and small black dots indicate epicenters of the largest foreshock (M j 6.4) and the mainshock (M j 7.3) of the 2016 Kumamoto earthquake sequence, and the aftershocks that occurred from April 14 through June 9, 2016, respectively. b Black triangles show the locations of temporary observation sites deployed in Mashiki (S1-S8). Gray triangles indicate the temporary aftershock observation stations deployed by Yamanaka et al. (2016) in the same area (MK01, MK02, MK05 and MK09) and a white triangle shows the location of seismic intensity meter deployed by JMA. Gray squares show the locations of the seismic intensity observation station inside Mashiki town office (MTO) and the KiK-net strong-motion observation station KMMH16 (same location as the Hi-net station N.MSIH) operated during both the largest foreshock and the mainshock. The calculated JMA seismic intensity values at the temporary sites and station KMMH16 (ground surface) during an event that occurred at 22:08 JST on June 12, 2016 (EQ18: M j 4.3, depth = 7 km) are also denoted in rectangles. Colored circles indicate the ratio of the totally collapsed wooden buildings identified during a field survey (Yamada et al. 2017a) microtremor observations in this area (e.g., Kagawa et al. 2017;Kawase et al. 2017;Nagao et al. 2017) showed a clear correlation between the patterns of peak frequency in the horizontal-to-vertical (H/V) spectral ratios and the landform characteristics. Yamada et al. (2017c) performed further microtremor observations with miniature arrays at 108 sites in Mashiki and clarified the differences of Rayleigh wave dispersion characteristics between the damaged and undamaged zones. They found that an extremely low S wave velocity (Vs) layer (< 100 m/s) exists at very shallow depth beneath Mashiki. The estimated thickness of this layer is less than several meters in the severely damaged zones, where the dominant frequencies of microtremor H/V spectral ratio are higher than 2 Hz. In contrast, in the floodplain of the Akitsu River where almost no buildings were collapsed, the low-velocity layer is thicker than 10 m and the dominant H/V frequencies were around 1 Hz.
The predominant frequency of soil layers is fundamental information in seismic microzonation and a reasonable indicator that relates the local site condition to earthquake damage. Past observations of damaging earthquakes in Japan showed that the strong ground motions in the frequency range between 0.5 and 1 Hz greatly affect wooden buildings (e.g., Kitagawa and Hiraishi 2004;Sakai et al. 2008). The natural frequency of low-rise wooden buildings is between 2 and 10 Hz (e.g., Sugino et al. 2016a), but the frequency may shift to lower values (around 1 Hz) when they are subjected to strong shaking. Therefore, larger damage was expected along the Akitsu river where the frequencies around 1 Hz are dominant. However, the damage investigations showed very different results with relatively light damage close to the river and more severe damage associated with sites of higher natural frequencies that were located farther from the river. Besides the local differences in site conditions, there are several other possible factors that could account for the damage concentration in Mashiki; the construction years of buildings (e.g., Mori et al. 2017;Sugino et al. 2016b;Yamada 2017), types of foundations in the buildings (e.g., Mori et al. 2017), and/or surface fault ruptures (e.g., Suzuki et al. 2016). There were no strong-motion records of the larger events in the heavily damaged zones so that aftershock data recorded in both the damaged and undamaged zones are used to help understand the difference of seismic response.
Recent studies have tried to explain this unusual site amplification in Mashiki using nonlinear site response analyses, assuming that frequencies between 0.5 and 1 Hz were dominant in the heavily damaged zones and lower frequencies were dominant in the less damaged zone. Yamada et al. (2017c) calculated the seismic response for both damaged and undamaged zones of Mashiki using shallow seismic velocity structure models obtained from microtremor observations with typical nonlinear properties of sand and clay. Shingaki et al. (2017) and Yoshimi et al. (2017) determined dynamic properties of surface soils including volcanic ash clay, using actual soil samples in Mashiki, and Nakagawa et al. (2017) performed time-history response analysis for the mainshock at two sites in Mashiki using the obtained dynamic deformation properties. For analyzing the ground motions of the mainshock, information from the borehole (252 m in depth) of the closest KiK-net strong-motion observation station, KMMH16 (see the location in Fig. 1b) was used. The station is operated by the National Research Institute for Earth Science and Disaster Resilience (NIED) and located about 1 km to the northeast of the severely damaged area. Velocity models obtained from the borehole logging are well constrained for the shallower depths (less than a few tens of meters) but have lower resolution for the deeper structure. Therefore, it is important to evaluate the deep sedimentary structure in the damaged area.
To investigate the differences of ground motion characteristics in the heavily and less damaged zones as well as station KMMH16, we conducted aftershock and microtremor (ambient noise) observations in and around the damaged zones of Mashiki (see the observation locations in Fig. 1b). Figure 1b also shows the locations of aftershock observation stations deployed by Yamanaka et al. (2016) (gray triangles: MK01, MK02, MK05 and MK09) and a seismic intensity station deployed by JMA (white triangle) just after the mainshock. Our aftershock observation sites were installed in the heavily damaged zones, while the existing temporary aftershock stations were located outside of this zone. In this paper, we show the differences of the ground motion characteristics between the contrasting zones and propose a velocity structure model of the entire area for seismic response analysis.

Field observations
The aftershock and microtremor observations were conducted from June 10 to June 13, 2016, 2 months after the mainshock. Dimensions of the observation area are approximately 500 m in the east-west direction and 750 m in the north-south direction, which corresponds to the area of our damage and microtremor surveys (Yamada et al. 2017a, c). This area is a fluvial terrace of the Akitsu River and the altitude increases as a function of distance toward the north from the river. We installed eight accelerometers in and around the damaged zones at intervals of approximately 150-200 m (S1-S8: see Table 1 and Fig. 1b); six sensors between R28 and the Akitsu River and two sensors near the former main building of the Mashiki town office (MTO). A seismic intensity meter managed by the local government was installed at MTO. We used JEP-6A3 (Mitsutoyo Corp.) over-damped accelerometers at sites S1-S4, and SMAR-6A3P (Mitsutoyo Corp.) over-damped accelerometers at the other sites. Each sensor was connected to a 24-bit data logger DATAMARK LS-8800 (Hakusan Corp.), with a sampling frequency at 200 Hz. All the sensors were installed on the ground surface, outside of buildings and power was supplied from external batteries. Both types of sensors have the same response characteristics, except different sensitivities (0.112 V/m/s 2 for JEP-6A3 and 0.224 V/m/s 2 for SMAR-6A3P). The pretests confirmed that all of the sensors have the correct frequency response to an input signal between 0.2 and 10.0 Hz for both horizontal and vertical components.
We distributed the sensors in areas with different building damage ratios; sites S1 and S2 were located in a zone with no collapsed buildings near the Akitsu riverside, and the other sites were located in the damaged zones. Figure 1b also shows the ratio of the totally collapsed wooden buildings that correspond to damage grade D5 (story failure) in the damage scale proposed by Okada and Takai (1999). The sensor at site S6 was deployed at a location where 77% of the wooden buildings were collapsed. The distance between site S7 and the station MTO was about 30 m so that we will discuss the ground motion differences between the mainshock observation station and the other temporary stations in the latter part. We excluded the data at site S4 in the following analyses due to a possible problem in the sensor cable, although signals related to earthquake ground motions were recorded.
We also used data from the surface and borehole components of KiK-net station KMMH16 and co-located Hinet station M.MSIH. We confirmed that the frequency responses from 0.14 to 20 Hz were very consistent for the two instruments in the borehole (see an example in Additional file 1: Fig. S1).

Aftershock data on July 12, 2016
Earthquake ground motion data from the aftershocks were manually extracted from the continuously recorded waveforms, based on the unified earthquake catalogue provided by JMA. We selected events with maximum JMA seismic intensity in the Kumamoto Prefecture higher than or equal to 1. Ground motion data from 30 events with clear P-and S wave signals were obtained (see the event information in Table 2 and the hypocenter distribution in Fig. 1a). Epicenters of the aftershocks are in the regions correspond to the middle and northern sections of the Hinagu fault zone, the west part of the Futagawa fault zone and the north side of Mt. Aso. The JMA magnitude (M j ) of the selected events ranged from 1.7 to 4.3 and the depths were relatively shallow (between 1.9 and 13.5 km).
The largest event with M j 4.3 (M W 4.2; EQ18 in Table 2) occurred in the southern Yatsushiro City (see the hypocentral location in Fig. 1a) at 22:08 JST on June 12. JMA seismic intensity 5-lower was measured nearby Sakamoto, Yatsushiro City. This event was one of the three events from April 20, 2016, to the present time (June 2018) with maximum JMA seismic intensity larger or equal to 5-lower in Kumamoto Prefecture. Epicentral distance from Mashiki is about 40 km, and JMA seismic intensity of the event was 1 at MTO. Figures 2 and 3 show the observed accelerograms and Fourier amplitude spectra of the M j 4.3 recordings at the observation sites S1-S8 except S4, together with those for the surface and borehole components of KiK-net station KMMH16. In Fig. 2, the accelerograms recorded at sites S1 and S2 clearly show the long period components and long duration with large amplitude in the horizontal components. The Fourier amplitude spectra at sites S1 and S2 have dominant peaks around 1 Hz in the horizontal components which are not seen at other sites (Fig. 3). Characteristics of the site amplifications Figure 4 shows the spectral ratios of the S wave (20.48 s segment from the S wave onsets) between our surface observation sites and the borehole station M.MSIH, for all the detected events. A cosine taper was applied to both sides of the S wave segments with a 5% width of the segment time length. Assuming that the input ground motions beneath our observation area and KMMH16 station are similar and the S wave is dominant in the segments, the spectral ratios represent the approximate borehole to surface S wave transfer functions of subsurface layers at each site. The shapes of the derived spectral ratios look slightly different between the NS and EW components at some sites possibly due to the complex velocity structure, but there are no significant differences in the dominant peak frequencies. The spectral ratios clearly show dominant peaks around 1-2 Hz at sites S1 and S2 with amplitudes greater than 20, while other sites show peaks at higher frequencies (> 3 Hz). This is consistent with the results of past microtremor surveys which showed site amplifications at around 1 Hz near the Akitsu River (Kawase et al. 2017;Nagao et al. 2017;Yamada et al. 2017c). Next, we obtained H/V spectral ratios at each observation site from the ambient noise sections of the recordings. The derived probability density functions for our continuous data show the measurable frequency range of microtremor down to 0.2 Hz (see an example in Additional file 1: Fig. S2). We selected 1-h continuous waveform data during a quiet period of the 4-day continuous records (between 0:00 and 1:00 JST on June 11). The three-component recordings were divided into 40.96 s window segments with 50% overlap, and fast Fourier transforms were applied to each segment. Then, the H/V ratios were computed and smoothed using a Parzen window with bandwidth of 0.1 Hz, and all the ratios were stacked together. Note that we removed segments with large amplitude irregular noise in advance and the total number of segments used for the analysis was more than 140 for each site. Figure 5 shows the derived microtremor H/V ratios, as well as the spectral ratios of the north-south and east-west horizontal components with the vertical component (NS/V and EW/V, respectively) at all observation sites. The derived spectral ratios indicate predominant peaks in the frequency around 1-2 Hz at sites S1 and S2, where there was no building damage. This feature is consistent with the surface/borehole spectral ratios derived from our aftershock data.  is also comparable with the results of single-station microtremor H/V spectral ratios (Yamada et al. 2017c).
We also derived earthquake radial-to-vertical (R/V) spectral ratios using the data of event EQ18 (M j 4.3 and M W 4.2) that contain surface waves in the coda. After rotation of the two horizontal components (NS and EW) to radial and transverse components with respect to each source-sensor azimuth, we selected 20.48 s segments from 10 s after the S wave arrivals. The spectral ratios of radial-to-vertical motions that would correspond to Rayleigh wave ellipticity were computed at each site. Here, the same smoothing filter as the microtremor H/V calculations above was applied. Figure 5 also shows the derived earthquake R/V ratios at each observation site. They show peak frequencies similar to the H/V spectra and the predominant peaks (ratios of 30 or greater) around 1-2 Hz are also seen at sites S1 and S2. The peaks at 1-2 Hz are also recognized at sites S3, S5, S6 and S8, but the ratios are substantially smaller (between 5 and 20) compared to those in the less damaged zone and other predominant peaks with larger values are seen at the higher frequency range (> 2 Hz).
Besides the predominant peaks at frequencies higher than 1 Hz, there are distinct peaks at 0.4 Hz for all the spectra, especially at S5-S8. This peak can also be observed in the microtremor H/V spectra with smaller amplitudes. The common dominant peaks of the microtremor H/V and earthquake R/V ratios at this lower frequency are possibly related to the response of deep sedimentary layers beneath the area, which will be discussed in a later section.

Rayleigh wave phase velocity
We estimated dispersion characteristics of Rayleigh wave phase velocities in the lower frequency range using continuous microtremor and event data. We applied the two-site SPAC (2sSPAC) method (e.g., Morikawa et al. 2004;Hayashi et al. 2013) for the microtremor data, assuming that the microtremor wave-field around the area is spatially and temporally uniform. We selected four sensor pairs around the heavily damaged zones (S6-S8: 167 m, S5-S8: 201 m, S3-S6: 238 m, and S3-S5: 349 m) for the data processing. We used the vertical component of microtremors recorded during the entire observation period (from June 10, 18:00 to June 13, 13:00). At first, 1-h continuous data were divided into 87 segments (time window of 40.96 s) without overlap. After removal of segments containing irregular noise or earthquake signals, fast Fourier transforms were performed on the datasets and the complex coherence functions between the selected site-to-site pairs were calculated. The coherence functions from all the segments were stacked and an hourly stacked coherence function was derived. If an hourly derived coherence function showed small values around 1 Hz in comparison with many other traces, indicating small coherence between the two datasets, those data were excluded. Finally, all the hourly derived coherence functions were summed and the averaged function (with the variance) was regarded as the SPAC coefficients for the sensor pair. The Rayleigh wave phase velocities were estimated from the derived SPAC coefficients using a fifth-degree polynomial that approximates the inverse function of the zero order Bessel function of the first kind. The frequency range we used for the determination of phase velocity depended on the sensor distance; it is determined so that the wavelength corresponding to the estimated phase velocity was 2-7 times the sensor distances. The estimated phase velocity dispersion curves from the four sensor-to-sensor pairs were combined to derive a single dispersion curve. Since each SPAC coefficient curve has the variance so that the dispersion curve has been derived by the inverse-variance weighted method. The coherence functions of microtremors are not always stable in the lower frequency band due to the low detection performances and the final accepted frequency range for the 2sSPAC analysis was between 0.98 and 1.27 Hz (Fig. 6a).
We also estimated the Rayleigh wave phase velocities from aftershock data (event EQ18) based on the semblance analysis of the S-coda wave (Fukumoto et al. 2004), assuming that the coda wave is dominated by surface waves. The phase velocities were estimated at frequencies from 0.2 to 0.9 Hz with an interval of 0.01 Hz, using narrow-band velocity waveforms (0.9 and 1.1 times the target frequency for the low-frequency and the highfrequency cutoffs, respectively) in the vertical component. The length of the time window for the analysis is twice the reciprocal of the target frequency, and the window was shifted every 0.5 s throughout the dominant wave train. For each time window, the apparent velocity as well as incident angle was calculated. The derived apparent velocities with sufficient semblance values (≥ 0.9) were averaged, and the phase velocity for each frequency was determined. Figure 6a shows the estimated dispersion curve for the Rayleigh wave phase velocity using both microtremor and event data. The error bars indicate that the standard deviations are low for the entire frequency range, except around 0.7 Hz. The estimated dispersion curves of Rayleigh wave phase velocity from the miniature microtremor arrays (Yamada et al. 2017c) at around sites S1 (array A059), S6 (array A026), S7 (array A001) and station KMMH16 (array A000) are also shown in the figure. The estimated phase velocities vary with the sites but are less than 1000 m/s at the target frequencies of the miniature microtremor surveys (> 1.7 Hz), which is consistent with our phase velocity results. For the semblance analysis, we obtained dispersive characteristics of phase velocities in the frequency range between 0.24 and 0.86 Hz. The estimated phase velocities from microtremor and event data also show similar continuous characteristics for the common frequency range, suggesting that our estimations are reasonable.

S wave velocity structure model beneath KiK-net station KMMH16
We refined the P-S logging model at the closest KiK-net station KMMH16 (see Table 3 and Fig. 6b) by using our dispersion curve of Rayleigh wave phase velocity. We performed velocity structure inversions using a genetic algorithm (GA) technique (Yamanaka and Ishida 1996) to obtain a velocity structure model that can explain the estimated Rayleigh wave dispersion curve.
There are several previous velocity structure models in this area. The P-S logging model at KiK-net station KMMH16 has information down to 252 m depth. Goto et al. (2016) tuned the KiK-net logging model (hereinafter called GOTO model; see Table 3)  Comparisons of estimated phase velocity dispersion curves, subsurface structure models and the corresponding theoretical H/V ratios. a Circles show the estimated phase velocities of Rayleigh waves using event data (white) and microtremor data (black), as well as miniature microtremor array data (Yamada et al. 2017c) for nearby station KMMH16 (blue) and the aftershock observation sites (S1: green, S6: orange, and S7: red). b Vs structure models used in this study. The P-S logging model at the nearest KiK-net strong-motion observation station KMMH16, GOTO model , J-SHIS model, SIP2017 model (Ogawa et al. 2017;Senna et al. 2017a, b) and our estimated model are shown in black, green, orange, blue and red traces, respectively. c Theoretical dispersion curves of the Rayleigh wave phase velocities calculated with models in (b). d Theoretical H/V ratios calculated for models in (b). Gray shaded area indicates the frequency range in which the dominant peaks of microtremor H/V ratios are detected surface/borehole spectral ratio, and used 75% of the Vs in the logging model for the 6th-12th layers due to the limited sensitivity to the deeper structure (Table 3). For deeper structures, the Seismic Hazard Information Station of NIED (J-SHIS: http://www.j-shis.bosai .go.jp) provides a nationwide velocity model based on the borehole logging, reflection and gravity surveys.
In the J-SHIS model, the velocity structure beneath the target area consists of six layers; three subsurface layers (600 m/s ≤ Vs ≤ 2100 m/s) along with three upper crustal layers (3100 m/s ≤ Vs ≤ 3400 m/s). Figure 6c shows the theoretical dispersion curves obtained from the KiK-net logging data, GOTO and J-SHIS models, respectively. The KiK-net P-S logging model clearly overestimates the Rayleigh wave velocity for almost the entire frequency range, suggesting the actual Vs values are much lower than those of the logging model. The GOTO model explains the dispersion curve reasonably well at the higher frequency range (3-10 Hz) that corresponds to the target frequency range of miniature microtremor survey at the site (Yamada et al. 2017c), but still underestimates the lower frequency range even after the 75% reduction of Vs. The observed dispersion curve in Fig. 6c suggests that the hard rock layer (Vs > 3000 m/s) as in the J-SHIS model is necessary to explain the low-frequency range of the dispersion curve.
We constructed an initial model for the GA in the following procedure. We directly use the GOTO model for the top five layers, since the theoretical dispersion curve explains the observed dispersion curve at the higher frequency range. We also use the GOTO model as the starting model for the 6th-12th layers in the inversion. Due to the uncertainty of the deeper structure, we use the thicknesses and Vs values of the J-SHIS model for the layers of bedrock (Vs ≥ 2700 m/s). In the preliminary analysis, we found that the Vs = 2700 m/s layer in the P-S logging model (234 m) is too shallow to explain our dispersion curve. Therefore, we assumed that the layer would be deeper than the depth of the borehole. Based on this assumption, we modified the properties of the 13th layer in the GOTO model to those of the 12th layer, and added an extra layer (Vs = 1470 m/s) underneath the borehole with an unknown depth.
In total, we have 18 layers for our initial velocity model: the 1st-5th layers of the GOTO model, 6th-13th layers of the GOTO model for which Vs values will be determined, 14th layer with Vs = 1470 m/s with unknown thickness, and four bedrock layers of the J-SHIS model with Vs = 2700-3400 m/s at the bottom (Table 3). Since there is a trade-off between Vs and thickness of layers, we fix the thickness for the 1st-13th layers and the Vs for the 14th layers just above the bedrock. Other observations such as receiver functions might reduce trade-offs between the velocity and thickness of each layer; however, the small events (M3 and smaller) in our dataset were not appropriate for this method. Therefore, we validate the reliability of the inversion results focusing on the predominant frequency of the theoretical microtremor H/V, as well as the travel-times of body waves between the ground surface and borehole sensor depth in comparison with the GOTO model. The determined parameters for the inversion are Vs values of eight layers and thickness of the 14th layer; the range of the parameter search for Vs is ± 50% of the initial value, and 0-1024 m for the thickness. Since there is no available information on the density profile at station KMMH16, we assigned the parameters based on the relation between Vp and density (Gardner et al. 1974). In the GA optimization, the number of generations that corresponds to the total number of iteration is 500 and the number of population in each generation is 30. Crossover and mutation rates were set to 0.5 and 0.01, respectively. We repeated this optimization process twice by setting the output of the first process as an input of the second process, so that the most optimal values are inside of the range of the parameter search. We used a subroutine program package DIS-PER80 (Saito 1988) to compute the theoretical dispersion curves for the generated models. The optimal velocity profile is shown in Fig. 6b and Table 3 with the phase velocities of Rayleigh wave, and theoretical H/V ratios are shown in Fig. 6c, d, respectively. The estimated values for Vs at the 6th-13th layers are substantially smaller than those of the logging model and GOTO model. The estimated depth of the bedrock layer is about 600 m. The comparisons of H/V ratios and phase velocities show that the agreements between the theoretical and estimated phase velocities are substantially improved. Also, the H/V spectral peak around 0.4 Hz is reproduced well by our model, which does not appear in the theoretical H/V spectra of the P-S logging model and GOTO model. Based on the quarter wavelength law, this peak reflects the velocity contrast of the bedrock, suggesting that the bedrock depth can be much deeper than the logging data (234 m).

Earthquake damage and site amplification
We discuss the relationship between the damage distribution and the largest aftershock (EQ18) during our observation period. For various ground motion measures, Fig. 7a, b shows the ratios of the peak ground acceleration (PGA) and peak ground velocity (PGV), and the differences of the JMA seismic intensity values between our observation sites and the KiK-net KMMH16 surface station (actual seismic intensity values are shown in Fig. 1b), respectively. Figure 7b also shows the ratios of the totally collapsed wooden buildings within 50 m radius at each site. Note that site S7 was located inside the town office where there were no houses nearby and the damage ratio was not available. This also shows that the ground motions of EQ18 are larger than those at other sites at sites S1 and S2 where there were no collapsed buildings around the sensors sites during the largest foreshock and mainshock. Figure 8 shows the acceleration response spectra for the EW component using the data of EQ18 at sites S1, S2, S6, S7 and KMMH16 together with the geometrical mean of the spectra (sites S1-S3, S5-S8). This clearly indicates that the different site response between the minimally damaged zone (S1 and S2) and the severely damaged zones (S6) as well as site nearby MTO (S7) and station KMMH16, and that the larger response at around 1 Hz can be seen only in the less damaged zone. Collapsed ratio b S1 S2 S3 S5 S6 S7 S8 S1 S2 S3 S5 S6 S7 S8 S1 S2 S3 S5 S6 S7 S8 NS EW UD NS EW UD Intensity difference Collapsed ratio Fig. 7 Comparisons of ground motions between each temporary observation site and KiK-net station KMMH16 (ground surface). a Circles show the ratios of three-component peak ground acceleration (PGA) and peak ground velocity (PGV) values. b Open circles show the differences of JMA seismic intensity values. Black circles indicate the ratio of the totally collapsed wooden buildings within 50 m radius from each observation site (Yamada et al. 2017a). Note that there are no wooden buildings near site S7 and the damage rate was not obtained As we expected from microtremor geotechnical surveys, our aftershock data show that ground motions around 1 Hz were amplified on the floodplain of the Akitsu River, while the amplification of this frequency range is much less in the damaged zone. This observation is inconsistent with the damage distribution of the mainshock, which shows severe damage away from the river. We need to consider other mechanisms, such as nonlinear site amplification, to explain the strong ground motions in Mashiki.

Comparisons of the velocity model with existing models
We determined a subsurface Vs structure model above the seismic bedrock with the derived dispersion curves of Rayleigh-wave phase velocities from our 4-day observation data in the entire area and miniature microtremor array at KiK-net station KMMH16. We assumed that the deep sedimentary layers for the damaged area are the same as in the vicinity of station KMMH16. If the assumption is correct, the dominant peak around 0.4 Hz in the earthquake R/V spectral ratios that corresponds to the velocity contrast at the seismic bedrock, should be seen in the KMMH16 data. Figure 9 shows R/V spectral ratios at station KMMH16 computed from the 19 shallow crustal earthquakes (> M4) that occurred prior to the 2016 Kumamoto earthquake (see the earthquake lists in Table 4). The length of the time window and the bandwidth for smoothing is the same as the R/V analysis. This averaged spectrum also shows the peak at 0.4 Hz, which suggests that the velocity contrast should exist over a relatively wide area, at least across the area of our aftershock observations and the KiK-net KMMH16 station.
We would like to compare our velocity model to other existing models. In the nationwide J-SHIS model, the depth of the bottom of the subsurface layers is 602 m, which is consistent with our estimation. Since the shallower part of the structure with low Vs (< 600 m/s) is not included in the J-SHIS model, the shallow velocity was greatly overestimated and the computed dispersion curve of Rayleigh waves beneath Mashiki has a large discrepancy with our observed dispersion curve (Fig. 6c). Also, the peak in the H/V spectrum at 0.4 Hz could not be explained (Fig. 6d). A group in NIED conducted microtremor observations with large-sized arrays after the Kumamoto earthquake and constructed a three-dimensional (3D) structure model of the deep sedimentary layers (Ogawa et al. 2017;Senna et al. 2017a, b). In their model (hereinafter called SIP2017 model), the velocity structure consists of five subsurface layers (350 m/s ≤ Vs ≤ 2100 m/s) and the bedrock layer (Vs = 2700 m/s) that was not included in the J-SHIS model. The velocity profile, theoretical H/V ratio, and dispersion curve of Rayleigh wave phase velocity computed from the SIP2017 model are shown in Fig. 6b-d, respectively. The depth to the bedrock layer in the SIP2017 model is much deeper than the P-S logging Period (s) S1 S2 S6 S7 KMMH16 Fig. 8 Acceleration response spectra with 5% of damping ratio derived from recordings of event EQ18 in Mashiki. Yellow, purple, blue and red traces indicate the response spectra at sites S1, S6, S7 and KMMH16. Dashed line shows the geometric mean of the response spectra derived using observed accelerograms at S1-S3 and S5-S8 model and the consistency between the estimated and theoretical phase velocity characteristics are reasonable in comparison with the J-SHIS model, although the peak frequency of the theoretical H/V is slightly higher than that of our observations. Our velocity model calibrated by the existing Vs profile at the KiK-net KMMH16 site ) and dispersion curve of Rayleigh wave phase velocity can explain the H/V peak reasonably well. In the SIP2017 model the depth to the Vs = 2700 m/s layer (1417 m) is much deeper than our estimation but the depth of the high-Vs layer (Vs = 2100 m/s) above the bedrock layer (597 m), which was not included in our model, is very close to our estimated bedrock depth, indicating that the difference in starting structure models for inversions produced different solutions for the bedrock depth. The exact depth of the bedrock for this area is still controversial, but at least both observations suggest that the depth of Vs > 2000 m/s layer is much deeper than that of the KiK-net logging data (234 m) and is expected to be about 600 m. Our aftershock data and velocity model help improve knowledge of the deep sedimentary structure model in this area.

Conclusions
We conducted aftershock and microtremor observations at eight sites in and around the damaged areas in Mashiki Town, Kumamoto Prefecture, 2 months after the 2016 Kumamoto earthquakes. The aftershock data indicate that site amplifications at approximately 1 Hz are significant during the small events in the less damaged zone where a thick low-velocity surface layer exists, whereas the dominant frequencies are higher than 3 Hz in the damaged zones. These trends are similar to the results from microtremor surveys in the area, indicating that the characteristics of ground motions for the mainshock and small earthquakes show opposite trends in Mashiki. Our observations support the idea that the dominant frequencies of the heavily damaged area were shifted to 1-2 Hz which has a great effect on wooden buildings during the largest foreshock and mainshock, while in the less damaged zone, the 1-2 Hz peaks would be shifted to much lower frequency range, which had less effect on the buildings. The estimated phase velocities of Rayleigh wave using microtremor data as well as the aftershock data (EQ18) are smaller than those of the P-S logging model at the KiK-net station KMMH16. The derived microtremor H/V and earthquake R/V ratios from our data show common dominant peaks around 0.4 Hz, which is possibly related to the response of deeper structure beneath the observation area. We performed a GA inversion for the deep sedimentary structure model with the derived phase velocities and the dominant frequency peaks to estimate an optimal structure model near station KMMH16. The estimated velocity structure model indicates that the Vs values at shallower depths in the target area are substantially smaller compared to existing models. Also, the depth to the bedrock layer might be much deeper (about 600 m), in comparison with the existing model from logging data (234 m). The earthquake R/V spectral ratios at station KMMH16 also show the common peak at 0.4 Hz, suggesting that the similarity of deeper structures between the KiK-net station and our observation area in the severely damaged zone. Therefore, mainshock data recorded in the borehole at station KMMH16 could be directly used as input motions for seismic response analysis in the damaged zone.

Additional file
Additional file 1. Fig.S1: Observed accelerograms and Fourier amplitude spectra at KiK-net KMMH16 and Hi-net N.MSIH stations. a Comparisons of the observed three-component acceleration records, b Fourier amplitude spectra at KiK-net borehole station KMMH16 (black traces) and Hi-net station N.MSIH (red traces) for event EQ18. Note that the original velocity seismograms observed at N.MSIH have been differentiated after the removal of instrumental responses. Fig. S2: Four-day long power spectral density of the recorded seismic noise at site S6. Here the power is expressed as decibel (dB), which corresponds to 10log 10 (m 2 /s 4 /Hz) of spectral acceleration.