Shear-wave velocity determination by combining data from passive and active source field investigations in Kumamoto city, Japan

We present the 1D subsoil structure and local site effects at KUMA strong ground motion station in Kumamoto City, Japan. We analyze data from a field campaign conducted in the framework of the Blind Prediction BP1 test of the 6th IASPEI/IAEE International Symposium: Effects of Surface Geology on Seismic Motion. In parallel with other participants of the BP1 test, we process data from passive and active source measurements aiming to determine the shear-wave velocity, Vs, structure and the site response at KUMA station. Passive measurements are associated to five microtremor arrays. In each array, seven seismometers have been deployed in a common-center triangle shape, recording microtremors simultaneously for 1 to 2 h. The vertical component of microtremors was analyzed using the spatial autocorrelation (SPAC) method. Cross-correlation coefficients were computed for all station pairs available for each array. By fitting the average SPAC’s coefficients to the first-kind zero-order Bessel function, J0, and assuming that microtremors primarily comprise fundamental mode Rayleigh waves, phase velocity dispersion curves were determined. Phase velocity values for frequencies > 15 Hz were obtained from data of a close-by active source geophone profile. We integrated the data with those of the passive measurements and obtained an experimental phase velocity dispersion curve. The resulting curve shows low velocity variation, from 150 to 200 m/s, in the surface layers, whereas significant dispersion appears in frequencies below 2.5 Hz. By inverting this curve, we achieved to determine the 1D shear-wave velocity structure at KUMA station. Site response characteristics were determined by applying the Horizontal-to-Vertical-Spectral-Ratios method. Significantly amplified peaks in the frequency range between 0.3 to 1.5 Hz dominate HVSR spectral ratios. These peaks correspond to resonant frequencies of soils and originate from different impedance contrasts within the substratum of the site.


Introduction
Local soil structure plays an important role in shaping the seismic wave-field, leading to amplification or/and de-amplification of different parts of the seismic motion spectrum and overall complicating wave propagation.So, to predict the characteristics of strong ground motion and study the seismic hazard of an area it is important to determine the structure and extract the velocity model of its local geology.This need is even more elevated when it comes to alluvial basins beneath or close to urban areas (e.g., Borcherdt 1970;Chávez-García and Bard 1994;Bonilla et al. 2002;Uetake and Kudo 2005;Lanzo et al. 2011).
One such deep alluvial basin is the Kumamoto Plain, in Kumamoto Prefecture, Japan.Within the Plain, many urban areas have been developed and seismic safety of such a densely populated environment calls for a good knowledge of the soil structure down to the bedrock.The greater area has been widely investigated, since it was the epicentral area of the strong 2016 Kumamoto earthquake sequence.It expands over 120 km from north to south and 40-60 km from east to west (GSJ-AIST 2004).Kinbo (Kinpo) Zan and Aso volcanoes surround the area to the north and east, with elevations above 1000 m (GSJ-AIST 2004).To the southeast of the Plain, two active tectonic features dominate, the Futagawa and the Hinagu fault zones, which are considered responsible for the 2016 Kumamoto earthquakes (e.g., JSEG 2016;GEER 2017;Mukonoki et al. 2016;Yamanaka et al. 2016).The Shirakawa, Midorikawa and Kiyama rivers, flowing through the Kumamoto Plain, are forming a deep alluvial terrain of different composition bed deposits (Ishizaka et al. 1995;Imanishi and Tamura 1958).Volcanic materials of the nearby volcanoes contribute to the stratigraphy of the area, making it quite complex (Mukunoki et al. 2016).
Due to the considerably active geological and seismotectonic regime of the Kumamoto region, permanent and temporary networks including strong and weak motion instruments have been implemented in

Graphical Abstract
different time periods.Most of them targeted to seismological studies, whereas others to soil structure determination and site effects estimation.In the following, we summarize the surveys that have been conducted close to our investigation site, which is the installation site of the KUMA strong motion station (Fig. 1).
NIED proposed a 3D model of the sedimentary layers in Japan, known as JSHIS model (NIED 2019).Based on this model, the investigated site is covered by thick soils with V S velocities < 1100 m/s down to 500 m depth.Layers with higher V S velocities, from 1100 to 3400 m/s, underlie the surficial sediments and continue down to 2000 m depth.Koketsu et al. (2009 and2012) proposed a 3D model for the whole island of Japan, known as JIVSM.This model characterizes subsurface layers down to 600 m depth with V S < 1300 m/s, and the deep structure by gradually increasing velocity values up to 3200 m/s.Chimoto et al. (2016) studied the shallow structure close to the target site by using small aperture microtremor measurements.Taking advantage of the temporary strong motion network deployed by Yamanaka et al. (2016) for seismological purposes, they recorded ambient noise at 26 sites within the Kumamoto Plain using triangular arrays of stations (array node locations are shown in Fig. 2).The closest arrays to KUMA site were deployed at stations KC03 and KC04.They determined shear-wave velocity over the first 30 m, V S 30, equal to 194 m/s for KC03 and 232 m/s for KC04.Yamanaka et al. (2016) analyzed weak motion records from stations KC03 and KC04 and they determined a resonance frequency of 1.66 Hz.Tsuno et al. (2017) investigated the local site effects in Kumamoto Plain along a 6-km-long north-south temporary network of seismic stations immediately deployed after the 2016 Kumamoto M j 6.4 earthquake foreshock (stations KR01 to KR06 in Fig. 1.Station KR01 was located at 240 m distance from the target site).By analyzing strong and weak motion data, they observed large amplifications in the HVSR spectral ratios of stations KR01-KR04 at frequencies of 1-3 Hz.They got similar resonant frequencies (1-2 Hz) from single-station microtremor measurements along the survey line KR01-KR06 (Fig. 1).
KUMA target site is located in the northern part of the Kumamoto Plain, very close to Shirakawa River bend (Fig. 1).The site is located 1000 m south of the Kumamoto station of JR Kyushu Railway Company (Chimoto et al. 2021).An accelerometer STR-361 (Takamisawa Cybernetics Co., Ltd) has been installed at the site since 2016.The station is located on thick Quaternary sediments of more than 500 m thickness, according to JSHIS model (NIED 2019).The accelerometer station has 3) and numerous events afterwards (Tsuno et al. 2021a).
This paper reports the results of our analysis on passive and active source measurements conducted around the KUMA station, in the framework of the Blind Prediction BP1 test organized within the frame of the 6th IASPEI/ IAEE International Symposium: Effects of Surface Geology on Seismic Motion.

Geological and seismotectonic setting of the investigated area
Kumamoto Plain is a terrain (Fig. 1) consisting of a diluvial plateau on the western slope of Mt.Aso (GSJ-AIST 2004) and an alluvial plain formed by river-bed deposits of the Shirakawa, Midorikawa and Kiyama rivers (Ishizaka et al. 1995).The plain is covered by thick alluvium layers, mainly consisting of sand, gravels, and pyroclastic materials, coming from the nearby Kinbo (Kinpo) Zan and Aso volcanoes (Imanishi and Tamura 1958;Mukunoki et al. 2016).Volcanic flows and ash of different composition and age were deposited interchangeably within soil layers (Mukunoki et al. 2016) making the stratigraphy of the area quite complex.Figure 1 shows part of the seamless geologic map of the Geological Survey of Japan (GSJ 2016) covering the broader area of the investigated site.Kumamoto city and most of the other urban developed areas are founded on the elevated Pleistocene fluvial terraces (Units 170, 171) and on pyroclastic flows (Units 95, 83).Towards northwest, the city lies on Pleistocene volcanic flows (Units 100, 101).Basement rocks are exposed mainly in the mountainous area to the north and east part.KUMA station is located on Pleistocene and Holocene marine and non-marine sediments (Unit 1).
The prevailing tectonic features in the Kumamoto Plain are the Futagawa and Hinagu faults (Fig. 1).Futagawa is an ENE-WSW trending right-lateral strike-slip fault zone consisting of three major segments (NIAIST 2012;HERPJ 2013).Its southwestern segment is located at about 12 km from our investigated site (Fig. 1).The NE-SW trending Hinagu fault zone, of 81 km length, is also divided into three main segments based on its geomorphological features and paleoseismic behavior.Both faults have been mapped and investigated in detail (e.g., Ikeda et al. 2001;Nakata et al. 2001;Nakata and Imaizumi 2002), especially after having been related to the 2016 Kumamoto earthquake sequence (among others JSEG 2016;GEER 2017;Mukonoki et al. 2016;Yamanaka et al. 2016).The 2016 sequence included a series of shallow earthquakes that started on April 15 with a magnitude MJ 6.5 event (MJ is the Japan Meteorological Agency magnitude) and crept up to the MJ 7.3 mainshock of similar faulting mechanism on April 16 (e.g., Kato et al. 2016;Asano and Iwata 2016).The earthquakes caused severe damage to houses and infrastructure in the epicentral area, especially in Mashiki, Nishihara, and Minamiaso communities, where instrumental ground motion reached the maximum level of 7 of the Japan Meteorological Agency seismic intensity scale (HERP-ERC 2016).

Field measurements
A field campaign with passive and active source measurements was conducted by the Local Organizing Committee of the ESG6 symposium on April 2019 in the vicinity of KUMA station (Fig. 2).The Organizing Committee distributed data of this campaign among the participants of the BP1 test.Passive measurements included microtremor records from five arrays.Each array comprises seven triaxial seismometers, which record simultaneously microtremors for a short period of time (1 to 2 h).An active source survey was conducted along a linear geophone profile in the nearby area (Fig. 2).24 vertical geophones were recorded the wave-field generated by the repeated impacts of a sledge-hammer to soil (source excitation), at both ends of the profile.A detailed description of the field campaign is given by Chimoto et al. (2021).

Array microtremor measurements
The geometry used for the passive arrays was that of common-centered equilateral triangles.Five arrays KUM-SS1, KUM-S, KUM-SM, KUM-M and KUM-LL with side lengths ranging from 1 to 960 m, were deployed around KUMA station (Fig. 2).Seven triaxial seismometers of 10 s natural period (SE-321, Tokyo Sokushin), connected to 24-bit digitizers (LS8800, Hakusan Corporation) were used in each array and microtremors were simultaneously recorded for 1 to 2 h with sampling rate of 200 Hz.An external GPS antenna synchronized stations every 1 h.
We analyzed the vertical component of microtremors, using the SPAC method (Aki 1957).From our analysis, we determined the phase velocity dispersion curve of the fundamental mode of Rayleigh waves.Basic signal processing was applied to raw data, including baseline correction and band-pass filtering between 0.1 and 25 Hz.Analysis was performed on 1-h time windows, even for the longer duration arrays (KUM-SM, KUM-M, KUM-LL), to facilitate direct comparison of results.Each 1-h record was divided into 34 smaller time windows of 200 s duration with 50% overlap of consecutive windows.The only exception to this division scheme applied to records of KUM-SS1 array, which were divided into 31, equally spaced, 100-s windows with 50-s overlap.Each of these time windows was passed through different narrow frequency band filters (Butterworth, 8 pole, zero-phase).The results were checked against their consistency and we verified that the solution does not depend on the applied versions of the filter.
For each station pair, corresponding to a specific interstation distance, SPAC's coefficients were computed using the filtered traces.The computations were straightforward.We averaged the cross-correlation coefficients as a function of frequency for all time windows, common to a given station pair.From the average values of each station pair we discarded those that do not follow the shape of the first-kind zero-order Bessel function, as well as those presenting low values (lower than 0.75) because they contradict the fundamental hypothesis that the recorded ambient noise wave-field is common to correlated stations.In addition, we discarded coefficients that did not tend to unity at low frequencies, attributing such deviation from unity to different sources of ambient vibration affecting different stations, something that was verified by visual checking data at corresponding time records.Correlation coefficients that did not go through zero, were considered as indicating that we should have larger inter-station distances to obtain useful information at examined frequencies.Such behavior was observed in correlation coefficients of most station pairs of the KUM-SS1 array.The reliability of the average SPAC's coefficients of each station pair was further checked with those of other station pair having the same inter-station distance belonging to same or other array.
Figure 3 shows example results in the form of crosscorrelation coefficients as a function of frequency for two station pairs of different inter-station distances as indicated in each plot.Each solid line corresponds to SPAC's coefficients determined for the 34 individual time windows.The coefficients correspond to station pairs S1-S2 of KUM-S array and LL1-LL2 of KUM-LL array.The majority of individual cross-correlation coefficients of these pairs resemble the shape of the J0 first-kind zeroorder Bessel function.Open circles are the average of the 34 values corresponding to a specific distance, each with its associated standard deviation.The first-zero crossing of SPAC's coefficients shifts from higher frequency values to lower frequency values with increasing distance between the stations, which is theoretically expected.SPAC's coefficients between two stations are expected to present large values only for frequencies for which the corresponding wavelength is large compared to their inter-station distance.Consequently, small inter-station distances better resolve high frequencies, whereas large inter-station distances constrain the lower frequency part of the distribution.
The average SPAC's coefficients as a function of frequency and distance, ρ(r, ω 0 ), were inverted to obtain the phase velocity dispersion curve, c(ω 0 ).The procedure includes the minimization of the difference between average SPAC's coefficients and J 0 (ω 0 /c(ω 0 )r) from the following equation reported in Aki (1957): By choosing a priori model parameters, an initial phase velocity dispersion curve was generated.The procedure was repeated many times using different values for the a priori model.The inversion was ended when phase velocity estimates remained constant and independent of the model parameters used.We made a joint inversion of the average cross-correlation coefficients simultaneously for all distances available in each array.This procedure allowed us to determine their phase velocity dispersion curve of the fundamental mode of Rayleigh waves.From these curves, we accepted phase velocities for wavelengths between the two limits imposed by Henstridge (1979).
Figure 4 shows the five phase velocity dispersion curves, which we determined from the five corresponding passive array measurements.In each dispersion curve plot, we have superposed the two lines corresponding to Henstridge's (1979) criterion.The first line corresponds to 2 times the minimum inter-station distance dmin of the particular array and the second line corresponds to (1) ρ(r, ω 0 ) = J 0 ω 0 c(ω 0 )r .15.7 times its maximum inter-station distance dmax.The figure shows that small distance arrays (KUM-SS1, KUM-S) provide information primarily in the high frequency range, whereas larger distance arrays (KUM-SM, KUM-M, KUM-LL) provide information at low frequencies.In overlapping frequency ranges, individual curves coincide and this supports the validity of our results.

Active source measurements
Seismic records obtained from the linear geophone profile (Fig. 2) were analyzed in the frequency domain.For this active source survey, 24 vertical geophones with natural frequency of 4.5 Hz (GS-11D, Geospace) connected to a 24-ch recording unit (Geode, OYO corp.) were deployed along a 34.5-m-long profile.Geophones were placed every 1.5 m distance and were set to record at a sampling rate of 1000 Hz.Seismic waves were artificially generated by vertical impacts of a sledge-hammer at two shot-points, located at 0.3 m from the two edges of the survey line (shot_0m and shot_34.5 m in Fig. 5).Each shot-point was triggered 10 times to produce repeated datasets.In our analysis we used seismograms of all datasets produced by the two shots (in total, 20 datasets) and of all geophones except those presenting amplitude saturation.Such behavior was shown by geophones g1-g4 adjusted to shot_0m and geophones g20-g24 adjusted to shot_34.5 m, and were thus excluded from the analysis.
Prominent dispersed surface waves recorded in seismograms were analyzed following McMechan and Yedlin's (1981) approach.Seismic records were baseline corrected and processed using the code of Herrmann (2004).Phase velocities were estimated using a p-ω stack, where p is the phase slowness and ω is the angular frequency.From the analysis of all datasets, we Fig. 3 SPAC's coefficients as a function of frequency for two different station pairs.Black lines correspond to individual cross-correlation coefficients computed for the 34 time windows.The distance between stations is reported within the parenthesis for each subplot.Open circles and associated error bars (red solid lines) indicate the mean SPAC's coefficients ± 1 standard deviation computed at each frequency, assuming normal distribution for each station pair obtained 20 phase velocity curves.Individual curves within the same shot-point dataset and between the two different shot-points, present good similarity.This suggests that subsoil structure along the geophone profile, is characterized by horizontal layers.We used an average phase velocity dispersion curve as representative of all datasets produced by the two shot-points.
Figure 5a shows seismograms recorded by the active source survey and the result of their analysis (Fig. 5b).Seismograms produced by the 9 th hammer impact when the source excitation was located at 0 m (shot_0m).We discarded seismograms of the first four geophones (g1 to g4) due to their amplitude saturation (Fig. 5a, inset), so we analyzed the remaining 20 seismograms (from g5 to ).Seismic traces are dominated by a large amplitude wavetrain of Rayleigh surface waves.Velocities of this dispersed wavetrain were calculated, in the time-distance domain, in the range of 155-170 m/s using the Goldstein and Minner (1996) software.Figure 5b presents results from the processing of the g5 to g24 seismograms in p-ω stack.Red colors denote the largest phase velocity stack values.The black curve connects the symbols inside the red-colored area and corresponds to the dispersion curve of the fundamental mode of Rayleigh waves.This curve suggests phase velocities from 150 to 170 m/s in the period range from 0.032 to 0.065 s, similar to those estimated in time-distance domain.Higher modes of Rayleigh waves were not detected.

Shear-wave velocity profile
Overall, the analysis of microtremors recorded in the different aperture arrays allowed the determination of phase velocities in the frequency range of 0.5 to 14.2 Hz.This is the range resolved by wavelengths imposed by Henstridge's (1979) criterion.Phase velocity estimates in higher frequencies, from 15 to 32 Hz, were determined by the active source measurements.We integrated the total of six extracted dispersion curves into one, to construct the most possible complete dispersion curve of the fundamental mode of Rayleigh waves for the investigated site.
Figure 6a shows the experimental dispersion curve proposed by our analysis for the KUMA station.The curve has an inflection point at frequency around 2.5-3.0Hz.Small variation velocities from 150 to 200 m/s characterize the curve for frequencies > 2.5 Hz, whereas significant dispersion appears for frequencies lower than the frequency of the inflection point.The curve was inverted, using the iterative inversion code of Herrmann (2004).An initial soil model consisting of a number of depth intervals (layers) was introduced in the code, described by layer thickness, V P and V S velocities, and density.Prior information in the neighbor of an investigated site, always helps to justify the parameters of the initial model, i.e., the number of layers constituting the soil model, V P and V S velocities, density ρ and V P -to-V S ratio of each layer.Some of the model's parameters may be appropriately changed from their initial value under user's judgment.Based on the available information around the target site (e.g., JSHIS model (NIED 2019), JIVSM model (Koketsu et al. 2009(Koketsu et al. , 2012)), we set up the soil parameters of the initial model.We started with a number of ten layers to describe the local geology.Density values were taken from the aforementioned models and fixed constant within each layer.V P and V S values of the individual soil layers were considered to increase with depth and V P -to-V S ratio was chosen equal to √3.This value is considered as average of the V P -to-V S ratio calculated for the individual layers of the aforementioned soil models in combination to empirical relations as well (e.g., Brocher 2005).The default value of damping was 1.0, and reduced to about 0.1 at any iteration step as convergence was approached.Starting from this initial model, the theoretical dispersion curve was determined and compared to the experimental one.Through an Fig. 6 Left: experimental phase velocity dispersion curve (black dots) produced by merging results of the passive and active source data analysis.The theoretical dispersion curve (solid red line) shows a very good fit to the observations.Right: V S velocity distribution with depth, determined from the inversion of the experimental dispersion curve.In the right part of the figure, the distribution of resolving kernels with depth is shown together with the damping ratio achieved iterative procedure, which included consecutive adjustments of model parameters, we finally accept the solution shown in Fig. 6.
The result of our inversion was deemed successful when (1) the theoretical dispersion curve fitted quite well the experimental one; (2) the damping value obtained quite small values, < 1.0E-7; and (3) the resolution Kernels corresponding to each layer showed an impulse like shape at the middle of that layer.The V S model of Fig. 6 describes a deep structure below KUMA station.It consists of eight soil layers with V S velocities from 151 to 1125 m/s, increasing downwards from the surface.Soft soils cover the surficial layers up to 40 m depth.The average V S over the shallowest 30 m, V S 30, was calculated at 200 m/s.Soil layers with velocities above 500 m/s appear at larger depths and down to 606 m, which marks the limit of the investigation depth of our analysis.Shearwave velocity of the half space was estimated to be larger than 2300 m/s.
To justify the results of our analysis, we compare them with those available for the target site.Phase velocity estimates and V S velocities of all partners' contribution of the BP1 test, have been statistically processed by the Local Organization Committee of the ESG6 (Chimoto et al. 2021(Chimoto et al. , 2023)).The average dispersion curve and average Vs model resulting from the statistical analysis are shown in Fig. 7a and b, respectively.Our dispersion results show a quite good coincidence with the average values for frequencies > 2.5 Hz, while more pronounced differences occur below this value (Fig. 7a).Differences in the dispersion curves predispose possible discrepancies in Vs velocities too, something that is obvious in Fig. 7b.A preferred Vs model down to 2000 m, has been proposed by Chimoto et al. (2021) and Chimoto et al. (2023).This model incorporates information of all shallow and deep structures, such as (1) the JSHIS model (NIED 2019); (2) the JIVSM model (Koketsu et al. 2009(Koketsu et al. , 2012)) 2016) close to KUMA site (stations KC03 and KC04 in Fig. 1).Down to 600 m depth, Vs velocities of our analysis are generally lower than the average ones at corresponding depths (Fig. 7b).Similar trends are shown in the preferred model relative to average one, as reported in Chimoto et al. (2021) and Chimoto et al. (2023).In conclusion, our soil model shows gradual increase of the Vs velocities of soils in depth.They show variations compared to statistically processed average Vs values of the corresponding depths and with the preferred model, which was used to predict Fig. 7 a Comparison between phase velocity dispersion curves determined by our analysis and the average dispersion curve ± 1 standard deviation resulted from the statistical analysis made in the framework of the BP1 test (Chimoto et al. 2021).Different colors used to plot different data sources.All experimental curves correspond to fundamental mode of Rayleigh waves.b Comparison between our Vs model (green line) and the average Vs velocities ± 1 standard deviation (black dots) determined by the inversion of the average dispersion curve of a weak and strong ground motion at the target area in the framework of BP2 and BP3 tests (Tsuno et al. 2021a, b;Tsuno et al. 2023).

Site response characteristics
Microtremors from the array measurements were further processed to estimate site response characteristics using the standard Horizontal-to-Vertical-Spectral-Ratios (HVSR) method introduced by Nogoshi and Igarashi (1971).This method has been applied worldwide (among others Nakamura 1989;Lermo and Chavez-Garcia 1994;Bard 1999; Chavez-Garcia and Kang 2014) and in many cases it provides a good estimate of site characteristics, in terms of resonant frequency and site amplification.The results of the method are considered reliable when the subsoil structure is relatively simple.Credibility of HVSR peaks increases with peaks' amplitude and, thus, clear peaks of amplification larger than 2 are usually considered in the interpretation of the HVSR results.
We computed HVSR spectral ratios to microtremor records of all stations (in total 35) deployed in the array measurements.We selected parts of the noise recordings, of 400 s duration and further divided them into 39 time windows of 20 s duration each, with an overlap of 50% between consecutive windows.Each 20-s window was baseline corrected, cosine-tapered (10%), band-pass filtered between 0.1 and 25 Hz and Fourier transformed.The amplitude Fourier spectrum was smoothed and the spectral ratios of the two horizontal components relative to the vertical were computed for each time window.The good consistency of the HVSR ratios of the 39 individual time windows of each horizontal component allowed us to calculate their average values.Based on the average ratios, frequencies of the main amplified peaks were determined.
In total, 65 average HVSR ratios were determined by the stations deployed around the KUMA station, eliminating those facing malfunction problems.Figure 8a and b shows as an example HVSR ratios of stations SS5 and LL6 installed in KUM-SS1 and KUM-LL arrays, respectively.In both plots, we superimposed the ratios determined for the 39 individual time windows as well as their average.Three significant peaks with variable site amplification (marked as 1st, 2nd and 3rd in Fig. 8a and  b) predominate in average ratios, in the frequency range between 0.3 and 1.5 Hz.The 1st peak is observed at frequencies 0.3-0.35Hz, the 2nd at frequencies 0.8-1.0Hz and the 3rd peak which has the strongest amplification is observed at 1.3-1.5 Hz.Although all three peaks are described sufficiently at the spectral ratios of station SS5 (Fig. 8a), the 2nd peak is absent in the HVSR ratios of station LL6 (Fig. 8b). Figure 8c summarizes 65 average HVSR ratios.In almost all of them, the three peaks are clearly described, suggesting a rather smooth spatial distribution of the soil's dominant frequency within a radius of 500 m around the KUMA station.A trough at 3 Hz is prominent to all spectral ratios, except those of station LL6, which is observed at higher frequency value, 6 Hz (Fig. 8c).
Based on the local geology of the site, we are trying to interpret the origin of the HVSR peaks.According to Vs velocities and thickness of soil layers constituting our model (Fig. 6), we concluded that the 2nd and 3rd peaks are caused by two impedance contrasts within sediments.The velocity contrast at 41 m depth is responsible for the 3rd peak at frequencies 1.3-1.5 Hz.The influence of this impedance contrast is supported by the trough/ peak frequency ratio with value around 2, constant to spectral ratios of all stations except those of station LL6 which takes double value (around 4).The 2nd peak at frequencies of 0.8-1.0Hz is associated with the impedance contrast at 96 m depth, while the 1 st peak at 0.3-0.35Hz seems to be caused by the interface between sediments and bedrock.It corresponds to the resonant frequency of the whole package soil deposit below KUMA station.
To validate our soil model, we calculated the 1D theoretical Transfer Function (TF) at the free surface of KUMA site.The TF was calculated by applying Kennett's reflectivity coefficient method (Kennett and Kerry 1979) to our soil model (Fig. 6), by making the assumption of one-dimensional horizontally layering and vertically incidence of S-waves.The amplified frequencies of the theoretical TF were compared with those of the HVSR ratios.Three significant peaks at 0.5 Hz, 1.0 Hz and 1.5 Hz, are shown in the theoretical TF (Fig. 8c).The peaks of 1.0 Hz and 1.5 Hz are comparable to frequency values of the 2nd and 3rd peaks, respectively, of the HVSR ratios.Difference in resonant frequency between empirical results, 0.3-0.35Hz, and theoretical, 0.5 Hz, supports the idea of a more thick structure than the one we determined for the KUMA site.Sediments thickness may extend deeper than 606 m, which is the investigation depth resolved by our analysis.
Our results are further compared to previous site response studies conducted in the Kumamoto Plain.The frequencies of the 2nd (0.8-1.0 Hz) and 3rd (1.3-1.5 Hz) peaks are well correlated with the results by Yamanaka et al. (2016), who processed strong motion data and determined a resonant frequency of 1.66 Hz (between stations KC03 and KC04 in Fig. 1).Tsuno et al. (2017) resulted in resonant frequencies in the range of 1.0-3.0Hz by analyzing weak motion data from a temporary, 6 km long, accelerometer network (along stations KR01 to KR04 in Fig. 1).The same authors obtained similar resonant frequencies, 1.0-2.0Hz by analyzing microtremor data.None of the aforementioned studies have determined the low-frequency peak at 0.3-0.35Hz, the resonance of soils corresponding to the whole sediments thickness.Similar frequencies to our HVSR peaks were obtained by other teams participating to BP1 test (Chimoto et al. 2021(Chimoto et al. , 2023)).

Conclusions
We have investigated the soil structure underneath the KUMA strong motion station, in Kumamoto Plain, Japan.We analyzed data from five microtremor arrays including inter-station distances from 0.5 to 962 m and seismograms from an active source, 34.5 m long, geophone seismic profile.All field measurements were provided in the framework of the Blind Prediction BP1 test of the 6th IASPEI/IAEE International Symposium: Effects of Surface Geology on Seismic Motion.
The vertical components of synchronized microtremor records were analyzed using the SPAC method (Aki 1957).We used data from seven triaxial seismometers deployed in different dimension arrays (KUM-SS1, KUM-S, KUM-SM, KUM-M and KUM-LL in Fig. 2).Cross-correlation coefficients (SPAC's coefficients) were calculated for each array and for all available station pairs, together with their average.The agreement among the curves proves the stationarity in time of the microtremor wavefield.We made a joint inversion of the average SPAC's coefficients simultaneously for all distances involved in each array and we determined the corresponding phase velocity dispersion curves.These curves correspond to the fundamental mode of Rayleigh waves.We accepted only phase velocity values for those wavelengths lying between the two limits imposed by Henstridge's (1979) criterion.
Seismograms from the active source seismic profile (24 vertical geophones installed every 1.5 m distance) were analyzed using the McMechan and Yedlin (1981) approach.We used all datasets of the sledge-hammer impacts produced by the two shot-points at the ends of the profile.We analyzed seismograms of all geophones, except those showing amplitude saturation.The processing produced 20 dispersion curves of similar shape.By combining the results of the geophone profile with those of the array microtremor measurements, we proposed an experimental phase velocity dispersion curve, representative for the target site.This dispersion curve corresponds to fundamental mode of Rayleigh waves and is characterized by phase velocities varying from 150 to 1920 m/s in the frequency range between 0.5 to 32 Hz.
This dispersion curve was introduced to iterative inversion software of Herrmann (2004) and the 1D shearwave velocity profile was determined down to 606 m depth.The analysis showed that the structure at KUMA site comprises thick sediments with gradually increasing velocity with depth.Shallow structure (0-41 m) is characterized by soft soils with Vs velocities ranging from 150 to 200 m/s.Below this depth, velocities are increasing to 500 m/s up to 96 m depth and get values of 600-700 m/s down to 606 m depth, which is the investigation depth achieved by our analysis.
We further investigate the influence of local site effects to seismic response of the site.We applied the HVSR technique to microtremor records from the microtremor arrays and we determined site predominant frequencies of soils.We analyzed multiple time windows of microtremors and examined both individual spectral ratios and their averages of the two horizontal components.Results across different time windows were found to be consistent.Three significant peaks corresponding to fundamental and higher modes (namely 1st, 2nd and 3rd) were determined in the spectral ratios.The 1st peak occurs at frequencies 0.3-0.35Hz, the 2nd peak at values of 0.8-1.0Hz and the 3rd peak with the highest amplitude, is observed at frequencies of 1.3-1.5 Hz.These peaks are clearly described at almost all stations of the arrays, suggesting a rather smooth spatial distribution of the soil's dominant frequency within a radius of 500 m around the KUMA station.
The stratigraphy of our soil model resulted to a trough/ peak structure.A trough appears at 3 Hz in the HVSR ratios of all stations deployed around the KUMA station (Fig. 8c).It shifts to frequency of 6 Hz, in the HVSR curves of station LL6.The trough at 3 Hz is very close to inflection point of the experimental dispersion curve at 2.5-3.0Hz (Fig. 6).A trough/peak frequency ratio of approximately two is calculated for all HVSR ratios, something that according to Konno and Ohmachi (1998) ensures the existence of high Poisson's ratio and high impedance contrast of the substructure.For station LL6, the value of the trough/peak frequency ratio equals 4 Hz.The shift of the trough/peak frequency ratio to higher frequencies together with the decrease of its amplitude may be related to a different soil structure below this station, something that could be supported by the location of this station from the other side of the Shirakawa River (Fig. 2).Similar assumptions were made to other studies too (e.g., Fäh et al. 2001).
We calculated the theoretical 1D Transfer Function (TF) at the free surface of the KUMA site, using the Kennett's reflectivity coefficient method (Kennett and Kerry 1979).The method was applied to the Vs model resulted by our analysis, by assuming horizontally stratified layering and vertical incidence of S-waves.We used the information determined from the HVSR with those of the theoretical TF, to validate our soil structure.The TF is dominated by three significant peaks at 0.5 Hz, 1.0 Hz and 1.5 Hz, which are comparable to frequency values of the 1st, 2nd and 3rd peaks of the HVSR spectral ratios.The 2nd and 3rd peaks are correlated to two impedance contrasts within sediments, at 41 and 96 m depth.It seems that most of the energy is trapped within surface sediments, something that is supported by the high amplification of these two peaks.The first contrast at 41 m depth, is strongly supported by the trough/peak frequency ratio of a value of two, similarly detected to all HVSR ratios of the stations deployed in the arrays, except of station LL6.The interface between sediments and bedrock seems to be responsible for the 1st peak in the HVSR spectral ratios, at 0.3-0.35Hz.The sediments-bedrock interface may lie at depths even larger than 606 m, the maximum investigation depth of our analysis, due to its difference with the corresponding value, 0.5 Hz, of the theoretical transfer function.The low value impedance contrast of the sediments-bedrock interface, seems to strongly influence the amplitude of this peak.
Our site response results accent the necessity of knowing the deep soil structure of a site and not only the surface shear-wave velocities, quite common information incorporated in seismic hazard estimates through V S 30.Local site effects are certainly controlled by the entire thickness of the sediments from surface to top of the bedrock.The shape of the spectral ratios resulted by our analysis supports the idea of a complex subsurface geology rather than a simple 1D structure under the target site.Site response at KUMA station seems to be affected, besides the local soil conditions, by the general geology and geometry of the basin (basin effects), and the existence of the geological discontinuities in the near area.This information should be incorporated to an overall study of site response in the vicinity of KUMA strong motion station.

Fig. 1
Fig. 1 Geology of the central part of the Kumamoto Plain (GSJ-AIST 2016) (https:// gbank.gsj.jp/ geona vi/).The location of the KUMA station is shown as black circle.The Shiragawa, the Midorikawa and the Kiyama rivers are drawn together with parts of the Futagawa and Hinagu faults.Locations of the field surveys related to our work, are shown as square symbols explained in the map's legend

Fig. 2
Fig. 2 Stations deployed in five microtremor arrays around KUMA station are shown as blue squares (KUM-LL in yellow, KUM-M in red, KUM-SM in green, KUM-S and KUM-SS1 out of scale, but located next to KUMA station).The adopted, common-centered, equilateral triangle geometry of the arrays has been indicatively drawn for KUM-LL array.The seismic profile location is shown as a solid white line to the north of KUMA station

Fig. 4 Fig. 5 a
Fig. 4 Phase velocity dispersion curves obtained from data of the five microtremor arrays deployed around the KUMA station.Dashed and dotted lines mark the limits imposed by Henstridge's (1979) criterion, which are associated to the minimum, dmin, and maximum, dmax, inter-station distances of each array.Dashed lines correspond to wavelengths equal to 2 times dmin, whereas dotted lines to wavelengths 15.7 times dmax.Arrays are distinguished by the use of different colors ; (3) the model proposed by Senna et al. (2018); and (4) Vs velocities determined by PS-logging measurements conducted in a shallow borehole near the KUMA station (Matsushima et al. 2021).Vs values of our analysis, as well as average statistically processed Vs velocities, differ from the preferred model, especially for the shallow layers, something that is obvious at Vs30 values; 200 m/s for our model, 305.6 ± 102.3 m/s for the average model and 260 m/s for the preferred model (based on PS-logging).Vs30 values of our model are well correlated to those determined by Chimoto et al. (

Fig. 8
Fig.8Examples of HVSR ratios computed for two stations: a SS5 of KUM-SS1 array and b LL6 of KUM-LL array.Each plot shows results from 39 individual time windows (grey lines) together with their average values (solid lines).In all plots, black line corresponds to EW component and red line to NS component.We marked with 1st, 2nd and 3rd the three main significant peaks observed in the spectral ratios.Plot c shows the average HVSR ratios calculated for all stations deployed around the KUMA station (65 in total).The theoretical 1D transfer function applying the Kennett's reflectivity coefficient method to the V S model of Fig.6is superimposed with the blue thick line