S-Wave Site Amplification Factors from Observed Ground Motions in Japan: Validation of Delineated Velocity Structures and Proposal for Empirical Correction

Abstract

We first derived site amplification factors (SAFs) from the observed strong motions by the Japanese nationwide networks, namely, K-NET and KiK-net of National Institute of Earthquake Research and Disaster Resilience and Shindokei (Instrumental Seismic Intensity) Network of Japan Meteorological Agency by using the so-called generalized spectral inversion technique. We can use these SAFs for strong motion prediction at these observation sites, however, we need at least observed weak motion or microtremor data to quantify SAF at an arbitrary site. So we tested the capability of the current velocity models in Japan whether they can reproduce or not the observed SAFs at the nearest grid of every 250 m as the one-dimensional theoretical transfer functions (TTF). We found that at about one-half of the sites the calculated 1D TTFs show more or less acceptable fit to the observed SAFs, however, the TTFs tend to underestimate the observed SAFs in general. Therefore, we propose a simple, empirical method to fill the gap between the observed SAFs and the calculated TTFs. Validation examples show that our proposed method effectively predict better SAFs than the direct substitute of TTFs at sites without observed data.

1. Introduction

The quantitative strong motion prediction with a source- and site-specific scheme is very important for mitigation of earthquake disaster and seismic design of important structures. It is especially true in Japan where large mega-thrust earthquakes are expected to occur within the coming 30 years. That is why we have a couple of nation-wide strong-motion networks in which a considerable number of strong motion records have been accumulated since 1996 [1].

There are several ways to simulate strong-motions as waveforms on the surface at a target site located at an arbitrary position. One is a theoretical Green’s function method (TGF) in which wave generation at the source, propagation from the source to the site, and amplification near the site are represented by the numerical modeling for the whole process from the source to the site. In this method we need a physical model of the medium to represent the wave propagation in the whole path. In other words, we need to calculate the theoretical Green’s function for a point source on the fault surface. The other is an empirical method in which we use observed ground motions of a small earthquake as a substitute for the Green’s function and sum up all the contributions from the elemental sources on the fault surface. It is called the empirical Green’s function method (EGF). If there are no appropriate small earthquake records to be used as the empirical Green’s function, we first generate synthetic waveforms based on many records of small earthquakes. It is called the statistical Green’s function method (SGF). Because the frequency range for the theoretical approach with coherent nature is limited to the lower end, usually below 1 Hz or lower, while the effective frequency range of EGF or SGF with inherent nature of stochasticity should be higher than that, a hybrid scheme with TGF and EGF or TGF and SGF are used naturally, as has been used in the current national project for strong motion predictions with specific sources [2].

After the deployment of the dense national strong motion observation networks, namely K-NET, KiK-net, and JMA Shindokei (Instrumental Seismic Intensity) network in Japan, a significant number of data has been accumulated. We can use these data to construct a model of SGF in a broadband frequency range. As long as we can generate the SGF for an arbitrary size of a small earthquake at an arbitrary location of a site in the frequency range of engineering interest, namely from 0.1 Hz to 20 Hz, we need not use a hybrid scheme. Thus we have been analyzing these strong-motion data in Japan by using the generalized spectral inversion technique (GIT) initially developed in 1980’s [3, 4] to delineate statistical properties of the three major terms, namely, the source term, path term, and site term [5, 6, 7, 8]. The novelty of our approach is that the hypothesized (i.e., extracted) seismological bedrock spectra at a reference site, YMGH01, are used as a reference to calculate site amplification factors at all the other observed sites. Such a separation of observed spectra into three major terms is sufficient to generate SGF at these observed sites. However, strong-motion simulations for the whole region near the seismogenic fault would be still difficult by SGF because we cannot estimate the site term at an arbitrary location other than the observed sites used in GIT. Thus, we need to develop a method to evaluate the site term at an arbitrary location as precisely as possible.

When we look at the site term as a function of frequency evaluated at K-NET, KiK-net, and JMA Shindokei network, we found that they show strong spectral fluctuations from 1 to 10 as a normal range of fluctuations and from 1 to 50 at tens of extraordinary sites with various peak frequencies. Several attempts have been made to correlate the primary characteristics of the observed site amplification factor (SAF) with a site proxy or proxies such as the S-wave velocity (Vs) averaged over top xx m, Vs_xx (e.g., Vs30) or the depth to the layer with the S-wave velocity higher than y.y km/s, Z_y.y (e.g., Z1.0) [9, 10, 11, 12, 13], trying to reproduce primary characteristics of SAF such as the fundamental peak frequency f₀ and its peak amplitude A₀. Unfortunately, these extracted characteristics are not sufficient to reproduce synthetic seismograms needed in the SGF summation. We should find a different strategy.

In what follows, we first introduce the fundamental characteristics of the observed SAF in the horizontal component (hSAF) derived from GIT [7, 8]. Then, we show comparisons of these hSAFs with the 1D theoretical ones calculated from the recently established unified velocity model (UVM) of the National Research Institute for Earth Science and Disaster Resilience (NIED) in the Kanto and Tokai regions. Next, we obtain the modification ratios to reduce the gap between them at the observation points and propose a scheme to evaluate hSAF at an arbitrary point by using the theoretical hSAF and the interpolated modification ratios, named FMR as the frequency modification ratio and AMR as the amplitude modification ratio. Finally, we propose an interpolation scheme to get hSAF in every 250 m grid point and validate the scheme at selected sites.

2. Observed SAF from GIT

In this section, we briefly introduce the observed horizontal SAF (hSAF) and vertical SAF (vSAF) derived from GIT [7, 8]. Here we introduce only their basic aspects because we are using their results as a starting point.

They restricted events and sites with the Japan Meteorological Agency (JMA)‘s magnitude M_JMA ≥ 4.5; source depth ≤ 60 km; hypocentral distance ≤200 km; peak ground acceleration ≤ 2 m/s²; and a number of observation sites triggered simultaneously for one event ≥3. These selection criteria resulted in 150,468 event-station pairs at 2,593 sites for 1,734 events. Only a relatively short duration of acceleration record from the onset of the S-wave was analyzed (5 s if 4.5 < M_JMA ≤ 6; 10 s if 6 < M_JMA ≤ 7; 15 s if 7 < M_JMA ≤ 8). A Parzen window of 0.1 Hz was used for a minimum level of smoothing. Note that the mainshock of the 2011 Off the Pacific Coast of Tohoku earthquake was excluded because the durations of those records are extraordinarily long.

As mentioned above, the most important feature of their GIT is that they determined the S-wave velocity structure at the reference site (YMGH01) using the transfer function (the spectral ratio and the phase difference) between the surface and the borehole 200 m below and that the observed Fourier spectra on the surface were deconvolved (divided by the amplification factor) to obtain the hypothesized outcrop spectra on the seismological bedrock with Vs of 3,450 m/s. The resultant S-wave velocity profile determined by the matching of the theoretical transfer function to the observed transfer function is quite similar to the original P-S logging data published by NIED [1], only with higher bedrock velocity of 3,450 m/s from 3,100 m/s. After the determination of the velocity profile, Nakano et al. [7] corrected (divided) all the observed spectra at YMGH01 by the calculated 1D S-wave site amplification factor on the surface with respect to the outcrop motion on the bedrock (=twice of the input) and used as the reference spectra in the subsequent GIT analyses, as if they were observed at YMGH01. Thus, their separated site terms, hSAF and vSAF, are considered to be the site amplification factors with respect to the outcrop seismological bedrock, on which there is virtually no site effect. Nakano et al. [7] successfully separated the source spectra and path terms as evidenced by their correspondence to the ω⁻² source spectra shapes and Q values similar to the previous studies in Japan.

Figure 1 shows examples of the observed hSAF at four sites in the Tokai region. We can see significant differences from site to site. Although we did not show vSAF here, the amplitude and its fluctuation of vSAF is much smaller than hSAF, especially below 3 to 4 Hz, that is, below the fundamental peak frequency of vSAF. That is why the earthquake horizontal-to-vertical spectral ratio, eHVSR, which is equal to hSAF/vSAF, tends to be similar to hSAF until the fundamental peak frequency of vSAF. However, to get hSAF from eHVSR, we need to correct vSAF, as recently proposed by Ito et al. [14]. Please note that precisely speaking, vSAF in this paper should be referred to as vSAF* as in [14] because we use the same reference condition for both hSAF and vSAF as the seismological bedrock spectra in the horizontal component so that we need to have correction due to the vertical-to-horizontal spectral ratio on the seismological bedrock on top of the vertical-to-vertical (i.e., P-wave) site amplification.

Figure 1.

Observed horizontal site amplification factor, hSAF, extracted from strong motions at K-NET, KiK-net, and JMA Shindokei network by GIT after [7].

3. Unified velocity model of NIED

To theoretically reproduce the observed hSAF and vSAF from GIT we need a velocity structure at each site from the seismological bedrock to the surface because they are the spectral ratios with respect to the outcrop spectra on the seismological bedrock. Note again that the seismological bedrock here is the surface of the crust on which we can assume no site amplification, whose S-wave velocity should be equal to or higher than 3 km/s. We have been delineating velocity structures in the deeper- and shallower-parts separately, primarily because we need to use different methods to explore the velocity structures in different depths. The boundary between them is the so-called engineering bedrock, whose Vs would be in between 350 m/s to 450 m/s. This is so beause we are using plenty of borehole information to constrain velocity structures in the shallower-parts, for which engineers need to gather information for their construction works. They usually need information only down to the layer with Vs in between 350 m/s to 450 m/s. However, it is apparent that higher-mode contributions of reverberated S- and P-waves within the whole basin above the seismological bedrock should show up in the frequency range higher than the fundamental peak frequency [15]. Therefore, we need a unified velocity model that integrates both shallower- and deeper-parts.

To that end, the National Research Institute for Earth Science and Disaster Resilience (NIED) has developed a unified velocity model (UVM) by integrating shallower- and deeper-parts of the structures above the seismological bedrock in the Kanto and Tokai regions [16, 17, 18]. The procedure to develop the model was based on the “concept of creating a subsurface structure model” released by the government agency, the Headquarters for Earthquake Research Promotion (HERP) [19]. Details of the procedure can be found in the papers referenced above, but the following is a brief description of the procedure.

An initial model of the shallow structure from the ground surface to the engineering bedrock was created based on existing studies and continuously collected SPT values in the boring data. Meanwhile, an initial model of the deep structure from the engineering bedrock down to the seismological bedrock was created based on the velocity models developed in existing studies by HERP. Then, an initial UVM was created by connecting them at the engineering bedrock. Next, the initial UVM was adjusted by using S-wave velocity structures at the strong-motion sites in the regions and those of spatially uniform and dense array microtremor explorations conducted as a part of Japan’s national Strategic Innovation Promotion project. Finally, the adjusted UVM was verified by using earthquake data at the strong-motion sites of NIED.

Examples of the important features of the resultant UVM are shown in Figures 2 and 3. Figure 2 shows the depth to the seismological bedrock with the Vs of 3.1 km/s (Z3.1) in the Tokai region, while Figure 3 shows the depth to the engineering bedrock with the Vs of 0.35 km/s (Z0.35), which is the interface between the shallower- and deeper-parts of the UVM in the Tokai region [18]. A similar map can be seen for the Kanto region in Senna et al. [16, 17].

Figure 2.

Depth contour in meters to the seismological bedrock with Vs of 3.1 km/s after [18].

Figure 3.

Depth contour of the engineering bedrock with Vs of 350 m/s after [18].

Please note that in the following sections when we calculate theoretical one-dimensional (1D) S-wave SAF, we use the following Q values for intrinsic and scattering attenuation:

$$\begin{array}{l}Q=Q_0\ast 0.5\mathrm{if}f\le 0.5\mathrm{Hz}\\ Q=Q_0\ast f\mathrm{if}0.5\mathrm{Hz}<f<5\mathrm{Hz}\\ Q=Q_0\ast 5\mathrm{if}5\mathrm{Hz}\le f\end{array}$$E1

The Q₀ values for the deeper part are listed in Table 1, while those for the shallower part we assume Q₀ = Vs/10.

No.	Vp (m/s)	Vs (m/s)	Density (kg/m^3)	Q0
1	1600	350	1850	70
2	1600	400	1850	80
3	1700	450	1900	90
4	1800	500	1900	100
5	1800	550	1900	110
6	2000	600	1900	120
7	2000	650	1950	130
8	2100	700	2000	140
9	2100	750	2000	150
10	2200	800	2000	160
11	2300	850	2050	170
12	2400	900	2050	180
13	2400	950	2100	190
14	2500	1000	2100	200
15	2500	1100	2150	220
16	2600	1200	2150	240
17	2700	1300	2200	260
18	3000	1400	2250	280
19	3200	1500	2250	300
20	3400	1600	2300	320
21	3500	1700	2300	340
22	3600	1800	2350	360
23	3700	1900	2350	380
24	3800	2000	2400	400
25	4000	2100	2400	420
26	5000	2700	2500	540
27	4600	2900	2550	580
28	5500	3100	2600	620
29	5500	3200	2650	640

Table 1.

Assumed layer profiles for the deeper part of UVM after [18].

To connect the bottommost layer of the shallower part Lsb with the topmost layer of the deeper part Ldt, we prioritize the shallower part if the depth Lsb is deeper than Ldt. If there is a gap between the two depths and Vs of Lsb is equal to or larger than Vs of Ldt, we extend Lsb down to Ldt. If Vs of Lsb is much smaller than Vs of Ldt, then we insert three layers with a gentle gradient of increasing Vs. Table 1 shows the parameters of layers assumed commonly in the deeper part. The bedrock S-wave velocity of UVM, 3,200 m/s, is close enough to that of the hypothesized bedrock of 3,450 m/s in GIT so we can compare both SAFs directly.

4. Observed and theoretical SAF

As mentioned in the previous section, the UVM of NIED is considered to be the most reliable velocity model for the strong motion simulation because the UVM combines all the available geophysical information to date related to the velocity structures from the ground surface to the seismological bedrock as densely sampled as possible. However, the actual S-wave SAF at a specific site, as shown in Figure 1, can be significantly different from the theoretical one calculated by a simple 1D S-wave multi-reflection theory for a stack of layers [20, 21, 22]. To see the difference, we plot in Figure 4 comparisons of the 1D theoretical hSAF with the observed hSAF at the same four sites in Figure 1. We calculate the theoretical hSAF as the 1D soil response on the surface of a combined velocity structure of the shallower- and deeper-parts with respect to the outcrop of the seismological bedrock motion (i.e., twice of the input at the bottom of the deeper part). Except for the site SZOH31, the theory tends to underestimate the observation.

Figure 4.

Observed hSAF extracted by GIT and 1D theoretical hSAF from UVM after [18]. Fundamental characteristics are well reproduced, however, the theoretical hSAF tends to underestimate the observed hSAF.

The major reasons for discrepancy are twofold; one is due to an inevitable inaccuracy of the derived velocity structure, and the other is due to a too simplistic assumption of the 1D horizontally-flat layered model. In the former there are two possible causes; one is the inaccuracy of the referenced values to delineate the structure including the assumed Q-values, and the other is the rapid spatial variations within the 250 m grid. In the latter there are two possible causes; one is the additional amplification due to the basin-induced surface waves generated at the edge of two or three-dimensional (2D or 3D) basins (see for example [23, 24, 25]) and the other is the topographic effects near the surface of irregular shapes such as a hill, a valley, or a cliff (see Kawase [26]).

To account for the effects of the basin-induced surface waves inside sedimentary basins, Nakano [8] and Nakano et al. [27] proposed to use an empirical ratio called the whole-wave-to-S-wave ratio (WSR), where the spectral ratios of the whole duration with respect to the S-wave portion with relatively short duration (5 to 15 s depending on the source magnitude as used in GIT) are averaged over all the observed events at a site. They found that the WSR tends to be close to unity irrespective of frequency for a site on hard rock, whereas it can easily exceed 10 in the lower frequency range for a site inside a soft sedimentary basin. Even for such a site, WSR will converge to unity in a higher frequency range above a few Hz. Because the spatial variation of WSR at one specific frequency highly correlates with that of the basin depths, as seen in Nakano et al. [27], Nakano [8] proposed a scheme to interpolate WSRs to make it possible to calculate a scenario-type hazard map with much higher spatial density than those of the strong motion observation sites. This WSR correction is a simple, empirical way to account for the additional amplification due to basin-induced surface waves.

However, other than the contributions of basin-induced surface waves through the modeling of WSR on top of hSAF as proposed by Nakano [8], it is difficult to account for the physical cause of the discrepancy between observation and theory at every sites as seen in Figure 4. We have been attempting to fill the gap by inverting the velocity structures from observed horizontal-to-vertical spectral ratios (HVSRs) of earthquakes under the diffuse field assumption [28, 29, 30], which is quite successful to reproduce observed HVSRs (and consequently also hSAF). This approach can provide us an equivalent 1D structure that will reproduce the observed hSAF at the target site quite precisely, however, the method is valid only for a site with either earthquake or microtremor records. We need a different strategy to evaluate hSAF as precisely as possible at an arbitrary site without any records. Because the velocity structure in the UVM of NIED is obtained with the spatial density of the 250 m grid, we want to reflect the fundamental characteristics of the theoretical transfer function of that structure, yet the resultant hSAF for strong motion simulation should be close enough to the observed hSAF.

5. Frequency and amplitude modification ratio

To overcome the difficulty to obtain better velocity models through physical parametrization of the underground structure at the sites without observed records, we would like to propose a different but simple approach here.

Suppose that a theoretical 1D hSAF from the UVM at a certain site deviates from the true hSAF by one of the aforementioned causes or their combination, the difference of the theoretical hSAF would be manifested in both the frequency axis and the amplitude axis. If we have stiffer or thinner layers in reality than the assumed profile in UVM, then the frequency characteristics would be all shifted towards the higher side. Or if we have a stronger velocity gradient within layers in reality, then the peak amplitudes would be all shifted towards the higher side. Thus, we have two different ways to adjust the theoretical hSAF to make it closer to the observed hSAF, one in the frequency axis and the other in the amplitude axis. Here is a simple way of correction for the theoretical hSAF, HSAF_the:

$$HSA{F}_{\mathit{mod}}\left(f\right)=AMR\ast HSA{F}_{\mathit{the}}\left(f/FMR\right)\text{,}$$E2

where FMR is the frequency modification ratio and AMR is the amplitude modification ratio. HSAF_mod is the resultant hSAF after both modifications as a function of frequency f. We need to determine FMR and AMR to make RES, the residual between HSAF_mod and the observed hSAF (HSAF_obs), minimum:

$$RES\left(FMR,AMR\right)={\displaystyle \sum _{f=\mathit{fmin}}^{\mathit{fmax}}\frac{\left|log10\left(HSA{F}_{\mathit{mod}}\right)-log10\left(HSA{F}_{obs}\right)\right|}{f}}$$E3

where fmin and fmax are the minimum and maximum frequencies of interest and we set them 0.12 Hz and 20 Hz, respectively. Here we use frequency f as a weight because the higher the frequency the denser the data in the linear space.

5.1 Grid search scheme

Because the calculation of Eq. (2) is quite easy, we use the grid search to obtain the best FMR and AMR combination. However, after several experiments, we found that the evaluation function in Eq. (2) seems too weak to determine FMR and AMR in a reasonable range because there is a trade-off between them. Therefore, we introduce the correlation function between HSAF_mod and HSAF_obs as an additional constraint. Then the target function to be maximized, GOF, becomes.

$$GOF=RE{S}_{\mathit{min}}/RES\left(FMR,AMR\right)\ast COR\left(FMR\right)/CO{R}_{\mathit{max}}$$E4

where RES_min is the minimum residual in the searching range, RES(FMR, AMR) is the residual shown in Eq. (2) as a function of FMR and AMR, COR(FMR) is the correlation coefficient between HSAF_mod and HSAF_obs, which is a function of only FMR, not a function of AMR, and COR_max is the maximum correlation coefficient in the searching range. Thus 1.0 is the maximum of GOF.

We set the searching range for FMR depending on the original correlation without frequency modulation, which is COR(1.0), as follows:

$$\begin{array}{c}\hfill \mathrm{If}0.6\leqq COR\left(1.0\right),0.80\le FMR\le 1.25\hfill \\ \hfill \mathrm{If}0.4\leqq COR\left(1.0\right)<0.6,0.67\le FMR\le 1.50\hfill \\ \hfill \mathrm{If}COR\left(1.0\right)<0.4,0.50\le FMR\le 2.00\hfill \end{array}$$E5

We use these searching ranges because if the correlation of the original model is sufficiently high, we should not modify its frequency characteristics so much. For AMR we set the searching range to be 0.333 ≤ AMR ≤ 3.00 irrespective of COR(1.0) because there are a few tens of sites with those amplitude differences as high as 3 times or as low as 1/3 and AMR does not alter COR(1.0). To efficiently search the best FMR and AMR with the precision of two digits, we employ the two-step grid search; first with every 0.1 increments, then with every 0.01 increments around the best FMR and AMR in the first step.

5.2 Evaluated FMR and AMR

Figure 5 shows examples of the resultant HSAF_mod in comparison with HSAF_obs and HSAF_the at the same K-NET, KiK-net, or JMA Shindokei network sites shown in Figures 1 and 4. As we can see, HSAF_mod determined by the grid search are quite close to HSAF_obs in both frequency fluctuations and amplitudes.

Figure 5.

Modified hSAF (HSAF_mod) by using FMR and AMR in comparison with observed hSAF (HSAF_obs) extracted by GIT and 1D theoretical hSAF (HSAF_the) from the velocity structure taken from UVM after [18] at the same sites in Figures 1 and 4. The optimal values of FMR and AMR for each site are shown in the upper-right corner.

Figure 6 shows the resultant optimal values of FMR and AMR at all the 546 sites used. We can see a very weak correlation between them. As for the searching range of AMR, namely 1/3 to 3, looks sufficient. On the other hand, we see a significant concentration of sites near the searching range boundary, 1/2 or 2 for FMR. This means that we could obtain better residuals and correlations if we search the optimal FMR in the frequency range wider than those limits. However, when we increase the searching range for FMR too much, we will see some cases where the reverberated fluctuations by the 1D resonance within the sediments seem to be shifted to the next overtone.

Figure 6.

Resultant values of FMR and AMR after the grid search for 546 sites in the Kanto and Tokai regions. Red triangles are those for the sites in Figure 5.

5.3 Correlation and residual improvement

If no improvement to the matching with the observed hSAF is achieved, our correction method does not have any merit. Therefore, we need to check if we can see significant improvement or not.

Figure 7 shows the distributions of the obtained correlation and the averaged residuals between HSAF_mod and HSAF_obs. We can see most of the site shows residuals less than 1.5 and correlations higher than 0.5. Table 2 shows the percentage of the sites in different categories in terms of the goodness of fit to HSAF_obs. When we compare the matching seen in Figure 5 and the distribution of these data in Figure 7, we can see that the average residual is more important than the correlation in terms of the quality of the modified hSAF because the correlation can be deteriorated easily by small fluctuations at different frequencies.

Figure 7.

Distributions of the correlation of HSAF_mod versus the averaged residual (the averaged spectral ratio between HSAF_mod and HSAF_obs) for the same 546 sites shown in Figure 6. Triangles are values at the sites shown in Figure 5.

Residual	Before	After	Correlation	Before	After
<1.25	5	43	<0.0	84	7
1.25–1.50	175	389	0.0–0.4	227	86
1.50–2.00	286	108	0.4–0.6	116	186
2.00–3.00	77	6	0.6–0.8	85	193
3.00<	3	0	0.8–1.0	34	74
Total	546	546	Total	546	546

Table 2.

Residual and correlation improvement in terms of the number of sites in each range before and after the correction.

Figure 8 shows significant improvements in the correlations from the original ones to the modified ones. There is no data with decreased correlations, however, there still remain 7 sites with the correlations less than zero, which was decreased from 84. Figure 9 also shows correlation improvements but as a function of FMR. We can see a clear concentration of FMR near the boundaries at 0.8 and 1.25, the boundaries of the searching range if the original correlation is larger than 0.6 as shown in Eq. (4). We can see a smaller improvement in this FMR range for higher correlation sites, in comparison to those with lower correlations.

Figure 8.

Improvement from the original correlations of HSAF_the to the modified ones of HSAF_mod. Triangles are values at the sites shown in Figure 5.

Figure 9.

Improvement from the original correlations of HSAF_the (before) to those of HSAF_mod (after) as a function of FMR. Data after the correction are concentrated near the boundaries of the searching range (0.8 and 1.25) when the original correlation is more than 0.6 because of the setting used. Note that FMR in the vertical axis of this figure is common for both before and after the correction.

5.4 Spatial interpolation

Now we have more than 500 sites in the Kanto and Tokai regions where we have determined modification ratios for frequency and amplitude, FMR and AMR. There are various ways to utilize these ratios for the prediction of site amplifications with much denser spatial resolutions based on the UVM in the 250 m grid. One way is to establish relationships of these modification ratios with respect to a site proxy or proxies such as Vs30 or Z1.0 as mentioned in the introduction. As seen in a lot of previous studies for site effects based on such relationships, however, we need to accept large deviations from the average relationships from site to site because it is the nature of the site amplification. We also face the possibility of the inaccurate choice of a site proxy in UVM for an arbitrary site used for modification ratios.

Therefore, we decide to use a direct spatial interpolation scheme as Nakano [8] proposed for WSR. In this scheme, we employ first GMT’s “surface” function [31, 32] in which the curvature minimization scheme is used together with the smoothing constraint of an elastic shell with the tension factor of 0.25. In this Step-1 interpolation, we use the 3 km equal-spaced grid. Then in Step-2, we use the 250 m grid to interpolate further by using the modified Shepard’s method [33].

Figure 10 shows the comparison of interpolated values in Step-1 with those used as targets for (a) FMR and (b) AMR. Interpolation in AMR is better because its spatial variation is smoother than FMR as shown later. Although we see tens of sites in Figure 10a away from the 1:1 line, the average deviation from unity for FMR is about 10% (1.1 or 0.9 times) and 91% of the interpolated values are within the range between 1/1.25 and 1.25 times of the original FMR. The number of outliers is misleading because we use the “blockmedian” command of GMT to refer to only the median value if we have plural sites in the same 3 km grid.

Figure 10.

Comparison of the original FMR and AMR used as targets of interpolations and those of interpolated values in Step-1 at strong motion observation sites. Red broken lines are linear regression lines whose inclinations and coefficients of determination R² are listed inside. (a) FMR and (b) AMR.

Similarly, Figure 11 shows the comparison of interpolated values in Step-2 with those used as targets from Step-1 for (a) FMR and (b) AMR. In Step-2 the interpolation is performed from the 3 km grid in Step-1 to the 250 m grid. The linearity of interpolation in Step-2 is much higher than that in Step-1 in the case of FMR, whereas that of Step-2 is as high as that of Step-1 in the case of AMR.

Figure 11.

Comparison of Step-1 FMR and AMR used as targets of interpolations in Step-2 interpolation and those of interpolated values in Step-2 at strong motion observation sites. Red broken lines are linear regression lines whose inclinations and coefficients of determination R² are listed inside. (a) FMR and (b) AMR.

We can see the spatial stability of the interpolation scheme as a gross picture shown in Figure 12 for FMR and AMR. These are the results of Step-2 with a spatial resolution of 250 m. Apparently, AMR is much smoother than FMR in terms of spatial variability so that the interpolation for AMR should be much easier and precise than FMR. On average, the Kanto region needs smaller correction values in FMR than those in the Tokai region, although it needs relatively larger correction values in AMR inside the whole soft-sedimentary areas in the north of Tokyo Bay.

Figure 12.

Interpolated contour maps of FMR (left) and AMR (right) after the Step-2 interpolation with the spatial resolution of 250 m.

5.5 Validation

So far we show that a simple two-step scheme of interpolation works to calculate modification ratios for both frequency and amplitude, namely, FMR and AMR from 546 strong motion stations in the Kanto and Tokai regions in a grid as small as 250 m. Because the UVM in these regions has a spatial resolution of 250 m, we can directly use these interpolated modification ratios once we calculate 1D theoretical hSAF at any of these grid points. To validate the method, we take four sites shown in Figure 5 out from the control points used for interpolation and let the program interpolate the modification ratios there and see how it works.

Figure 13 shows contour maps of FMR and AMR without four points used as examples in Figure 5. We can see smooth interpolation is achieved at these four points. Figure 14 shows the correspondence of original and interpolated values of FMR and AMR for the cases with and without four points. In case of FMR, the original and interpolated values are close to the 1:1 line and the pure interpolation values at three sites out of four are close to the original ones. We should note that the worst site of the interpolated FMR, AIC001, and the JMA site ___E34 were very close to each other as shown in Figure 13. In case of AMR, at three sites the interpolation values without four sites were not as good as those with four sites. Still the deviation from the original values are within the range of 1.5 or 1/1.5 times.

Figure 13.

Interpolated contour maps of FMR (top) and AMR (bottom) after the Step-2 interpolation without four points shown by triangles. Inside the triangles interpolated values at these sites are shown by color-coded circles.

Figure 14.

Comparison of the original FMR and AMR used as targets of interpolations and those of interpolated values in Step-2 with and without four example sites in Figure 5. Black crosses are original values at four sites and red circles are interpolated values without referring to those original values, that is to say, purely interpolated values. Except for one site for FMR (AIC001), our interpolation scheme works.

Finally, we compare the corrected hSAF at those four sites not used in the spatial interpolation as references but used as the targets of the interpolation with the observed hSAF in Figure 15. Although the corrections by the original FMR and AMR seen in Figure 5 are much better than the corrections by these interpolated FMR and AMR in this figure, the interpolated corrections still make theoretical hSAF closer to the observed hSAF.

Figure 15.

Modified hSAF (HSAF_mod) by using FMR and AMR after interpolation (without using these four site) in comparison with observed hSAF (HSAF_obs) extracted by GIT and 1D theoretical hSAF (HSAF_the) from the velocity structure taken from UVM at the same sites in Figures 1,4 and 5. The values of FMR and AMR by interpolation for each site are shown in the upper-right corner.

6. Conclusions

In order to evaluate an equivalent 1D S-wave site amplification factor at an arbitrary point, we propose an empirical method of correction on to the theoretical site amplification factor calculated from the unified velocity model of NIED for the Kanto and Tokai regions where the shallower- and deeper parts of the velocity structure are combined. First, we check how well the current unified velocity model in Japan can reproduce horizontal site amplification factors derived from the observed strong motions in the form of the equivalent 1D S-wave theoretical transfer functions at the nearest grid of every 250 m. The observed site amplification factors were obtained by GIT relative to the reference spectra extracted as the outcrop motions on the seismological bedrock. To be consistent with these observed site amplification factors, the theoretical transfer functions are calculated relative to the outcrop motions (twice of the input) on the seismological bedrock. We find that at about one-half of the sites the calculated 1D amplification factors show more or less acceptable fit to the observed ones, however, they tend to underestimate the observed amplifications in general. Therefore, we propose a simple, empirical method to fill the gap between the observed site amplification factors and the calculated ones based on the frequency and amplitude modification ratios. Once we obtain these modification ratios at 546 observed sites, we can interpolate them in space to obtain the modification ratios at an arbitrary point. Validation examples show that our proposed method effectively predict better site amplifications than the direct substitute of theoretical amplification factors at a site without observed data.

In the future investigation, we will apply the proposed correction method to the sites where we have observed records of either earthquakes or microtremors but we do not include them in the delineation of the modification ratios in order to further validate the effectiveness of the proposed method.

Acknowledgments

This study has used the strong-motion observation records from the K-NET and KiK-net of the National Research Institute for Earth Science and Disaster Resilience (doi: 10.17598/NIED.0004) as well as the seismic intensity (Shindokei) network of the Japan Meteorological Agency (https://www.jma.go.jp/jma/indexe.html, last accessed November 2020). A part of this study was supported by the Japan Society for the Promotion of Science (JSPS) Kakenhi Grant-in-Aid for Basic Research (B) Number 19H02405. Continuous support from Hanshin Consultants Co., Ltd. is highly appreciated.

Conflict of interest

The authors declare that they have no conflict of interest.