Site Ampliﬁcation during Strong Earthquakes Investigated by Vertical Array Records

: A number of vertical array records during eight destructive earthquakes in Japan are utilized, after discussing criteria for desirable requirements of vertical arrays, to formulate seismic ampliﬁcation between ground surface and outcrop base for seismic zonation. A correlation between peak spectrum ampliﬁcation and V s (S-wave velocity) ratio (base V s /surface V s ) was found to clearly improve by using V s in an equivalent surface layer wherein predominant frequency or ﬁrst peak is exerted, though the currently used average V s in top 30 m is also meaningful, correlating positively with the ampliﬁcation. We also found that soil nonlinearity during strong earthquakes has only a marginal effect even in soft soil sites on the ampliﬁcation between surface and outcrop base except for ultimate soil liquefaction failure, while strong nonlinearity clearly appears in the vertical array ampliﬁcation between surface and downhole base. Its theoretical basis has been explained by a simple study on a two-layered system in terms of radiation damping and strain-dependent equivalent nonlinearity.


Introduction
In determining design motions for superstructures, site amplifications have to be evaluated in soil layers above an engineering bedrock (base layer) at which the design motion is prescribed. The bedrock is defined at stiff bearing strata with S-wave velocities of around V s = 400~700 m/s or larger, normally. The incident earthquake motions are evaluated from fault mechanisms and attenuations along wave paths from faults to particular sites for individual earthquakes. Then, the evaluation of seismic amplification due to different site conditions is essential in developing hazard maps for local governments.
In evaluating the site amplification, the soil profile is simplified by a one-dimensional model unless it is very variable in geology or topography. Horizontal shear wave (SHwave) propagating in the horizontally-layered soil deposits is first considered because of its significant effect on geotechnical engineering problems, although other wave types have to be considered in more heterogeneous or topographically complex site conditions. The appreciation of the significant role of the SH-wave in site amplification stems from the pioneering work by Kanai et al. in 1956 [1] and 1966 [2], wherein simultaneous earthquake observations were conducted at the ground surface and the subsurface level below in a tunnel, demonstrating that seismic horizontal ground motions can be idealized essentially by the SH-wave vertically propagating in one-dimensional horizontal layers.
The site amplification in the 1D soil profile seems to depend on S-wave velocities, soil densities, internal damping of the individual layers. Furthermore, strain-dependent nonlinear properties may affect the amplification in soft soils in focal areas during strong Figure 1. 1D response analysis of horizontally layered soil model (a), and peak amplification between base and surface of 2As/2Ab versus Vs-ratio (b) (Shima 1978 [3]).
This finding further implies that amplification at a site relative to another may be evaluated from the ratios of the top layer Vs at the two sites if they share the common base layer of identical Vs-value. This seems to be compatible with a basic concept in currently employed formulas using near-surface S-wave velocities in the seismic zonation. However, a Vs-value not necessarily of the top surface layer but an average Vs-value of multiple layers to a certain depth is chosen in practice, because a surface layer is sometimes too thin to represent a site properly.
Thus, a couple of empirical formulas were proposed using Vs -values near the ground surface, as summarized in Figure 2. The empirical formula by Joyner and Fumal (1984) [4] evaluates maximum surface acceleration (PGA), maximum velocity (PGV) or 5% damping response spectrum using average Vs in surface layers based on strong motion data obtained in the United States. Here, the Vs-value is averaged over the soil thickness corresponding to a 1/4-wavelength from the ground surface, assuming seismic motions with the predominant period of T = 1.0 s. Another empirical formula was proposed by Midorikawa (1987) [5] based on earthquake records in Japan, wherein the ratio of maximum velocity amplitudes between ground surface and the base was correlated with the S-wave velocity averaged over surface top 30 m, Vs30 (m/s). A similar empirical formula was also proposed by Borcherdt et al. (1991) [6] in the US incorporated records during the 1989 Loma Prieta earthquake in California and associated small shocks, wherein the spectrum amplifications averaged over the period T = 0.4-2.0 s, were correlated with the top 30 m S-wave velocity Vs30 (m/s). Thus, there have been a couple of site amplification evaluation methods currently used in seismic zonation practice.  (Shima 1978 [3]).
This finding further implies that amplification at a site relative to another may be evaluated from the ratios of the top layer V s at the two sites if they share the common base layer of identical V s -value. This seems to be compatible with a basic concept in currently employed formulas using near-surface S-wave velocities in the seismic zonation. However, a V s -value not necessarily of the top surface layer but an average V s -value of multiple layers to a certain depth is chosen in practice, because a surface layer is sometimes too thin to represent a site properly.
Thus, a couple of empirical formulas were proposed using V s -values near the ground surface, as summarized in Figure 2. The empirical formula by Joyner and Fumal (1984) [4] evaluates maximum surface acceleration (PGA), maximum velocity (PGV) or 5% damping response spectrum using average V s in surface layers based on strong motion data obtained in the United States. Here, the V s -value is averaged over the soil thickness corresponding to a 1/4-wavelength from the ground surface, assuming seismic motions with the predominant period of T = 1.0 s. Another empirical formula was proposed by Midorikawa (1987) [5] based on earthquake records in Japan, wherein the ratio of maximum velocity amplitudes between ground surface and the base was correlated with the S-wave velocity averaged over surface top 30 m, V s30 (m/s). A similar empirical formula was also proposed by Borcherdt et al. (1991) [6] in the US incorporated records during the 1989 Loma Prieta earthquake in California and associated small shocks, wherein the spectrum amplifications averaged over the period T = 0.4-2.0 s, were correlated with the top 30 m S-wave velocity V s30 (m/s). Thus, there have been a couple of site amplification evaluation methods currently used in seismic zonation practice.
More recently, the world has experienced quite a few destructive earthquakes. Most have been recorded at the ground surface of different geologies, including outcropping engineering bedrock as horizontal array records (e.g., BSSC 2003 [7]). In addition, vertical array earthquake observation systems have been increasingly deployed, particularly in Japan, wherein site-specific amplifications can be investigated at the same sites along with the depth, having various geological profiles between the ground surface and downhole base. Among them, a number of vertical arrays named KiK-net have been systematically deployed at more than 700 sites all over Japan after the 1995 Kobe earthquake by NIED (National Research Institute for Earth Science and Disaster Resilience in Tsukuba, JAPAN). Each KiK-net recording site consists of a pair of three-dimensional accelerometers (EW, NS, UD), one at the ground surface and another downhole deeper than engineering bedrock, enabling to measure the site amplification between the two [8]. It is generally accepted that the engineering bedrock can serve as a common scale to define relative site amplification among sites sharing the same bedrock, though deeper geological structures may have some effect on absolute amplification (e.g., Trifunac 1990 [9]). More recently, the world has experienced quite a few destructive earthquakes. Most have been recorded at the ground surface of different geologies, including outcropping engineering bedrock as horizontal array records (e.g., BSSC 2003 [7]). In addition, vertical array earthquake observation systems have been increasingly deployed, particularly in Japan, wherein site-specific amplifications can be investigated at the same sites along with the depth, having various geological profiles between the ground surface and downhole base. Among them, a number of vertical arrays named KiK-net have been systematically deployed at more than 700 sites all over Japan after the 1995 Kobe earthquake by NIED (National Research Institute for Earth Science and Disaster Resilience in Tsukuba, JAPAN). Each KiK-net recording site consists of a pair of three-dimensional accelerometers (EW, NS, UD), one at the ground surface and another downhole deeper than engineering bedrock, enabling to measure the site amplification between the two [8]. It is generally accepted that the engineering bedrock can serve as a common scale to define relative site amplification among sites sharing the same bedrock, though deeper geological structures may have some effect on absolute amplification (e.g., Trifunac 1990 [9]).
A factor potentially influencing the site amplification during nearfield destructive earthquakes may be the nonlinearity of soil properties in soft soil sites. During the 1964 Niigata earthquake in Japan, acceleration records showed extraordinary low-amplitude long-period motions (Kanai 1966 [10]), presumably reflecting the base-isolation of intensively liquefied ground. During the 1995 Kobe earthquake, Port Island vertical array recorded evident de-amplification of horizontal acceleration due to extensive liquefaction of reclaimed soils (Kokusho & Matsumoto 1998 [11]). Thus, the soil nonlinearity effect has been recognized to be considerable during strong earthquakes. There have been several studies to investigate the soil nonlinearity effect by utilizing KiK-net records of strong and weak motions (e.g., Kokusho & Sato 2008 [12], Kokusho 2013 [13], Régnier et al., 2013 [14]), though it is not clearly taken into account in a hazard map or seismic zonation due to absence of sufficient observational evidence.
In the present paper, eight strong earthquakes in 2000-2008 recorded at many KiKnet vertical array stations throughout Japan are addressed to investigate how their spectrum peak amplifications between surface and engineering bedrock are correlated with the soil profiles and associated Vs there. In preparing for that, how to properly make use of the vertical array data in conducting the amplification study is first discussed. Then, empirical relationships between the amplification and Vs-ratios are explored by  A factor potentially influencing the site amplification during nearfield destructive earthquakes may be the nonlinearity of soil properties in soft soil sites. During the 1964 Niigata earthquake in Japan, acceleration records showed extraordinary low-amplitude long-period motions (Kanai 1966 [10]), presumably reflecting the base-isolation of intensively liquefied ground. During the 1995 Kobe earthquake, Port Island vertical array recorded evident de-amplification of horizontal acceleration due to extensive liquefaction of reclaimed soils (Kokusho & Matsumoto 1998 [11]). Thus, the soil nonlinearity effect has been recognized to be considerable during strong earthquakes. There have been several studies to investigate the soil nonlinearity effect by utilizing KiK-net records of strong and weak motions (e.g., Kokusho & Sato 2008 [12], Kokusho 2013 [13], Régnier et al., 2013 [14]), though it is not clearly taken into account in a hazard map or seismic zonation due to absence of sufficient observational evidence.
In the present paper, eight strong earthquakes in 2000-2008 recorded at many KiK-net vertical array stations throughout Japan are addressed to investigate how their spectrum peak amplifications between surface and engineering bedrock are correlated with the soil profiles and associated V s there. In preparing for that, how to properly make use of the vertical array data in conducting the amplification study is first discussed. Then, empirical relationships between the amplification and V s -ratios are explored by incorporating a number of vertical array records of mainshocks during the eight earthquakes. Finally, strain-dependent soil nonlinearity effect on the amplification is investigated by combining mainshocks and aftershocks recorded at the same sites, and an intriguing result obtained there will be examined in a simple theoretical study to find a global picture of soil nonlinearity effect on site amplification.

Vertical Array Compared to Horizontal Array
In general, there can be two different array systems to measure the earthquake amplification at a site between the ground surface and base layer as illustrated in Figure 3.

Vertical Array Compared to Horizontal Array
In general, there can be two different array systems to measure the earthquake amplification at a site between the ground surface and base layer as illustrated in Figure  3.
in the soft soil and stiff base layer respectively, as a function of angular frequency ω but abbreviated here, travelling toward positive z-direction with S-wave velocity Vs1 and Vs2, as indicated in Figure 3. The site amplification between the two different geologies can be directly evaluated 2As/2Ab if upward wave amplitude Ab in base is assumed to be basically unchanged in the area. This amplification 2As/2Ab is what we need to develop a seismic hazard map for seismic zonation. (b) Vertical array consisting of a surface and down-hole seismometers at the same place with exactly the same bedrock upward wave Ab. This can evaluate site amplification at exactly the same location but does not directly provide the amplification 2As/2Ab but another value 2As/(Ab + Bb) because the downhole record is more or less contaminated by the downward wave Bb from the overlying layers (e.g., Schnabel et al., 1972 [15], Kokusho 2017 [16]). Here, Bb are amplitudes of downward waves in the soft soil layer. Hence, modification of 2As/(Ab + Bb) in the vertical array is necessary to evaluate the amplification 2As/2Ab for seismic zonation.
In conducting the modification, it is required to pay attention to another important feature of the vertical array different from the horizontal array as described below.    Table 1. The installation levels of surface and downhole seismometers are indicated with the pair of arrows in site profiles. The V s -profiles (included in the figure) obtained by in situ S-wave logging are available in the KiK-net vertical array database. The two thin curves also overlaid in the charts represent observed Fourier spectrum ratios (the ratio of Fourier spectrums of observed motions between surface and downhole) in EW and NS directions for earthquake records obtained there. At all the sites, peak frequencies of 2A s /(A b + B b ) are mostly compatible with those of observed spectra in lower frequency peaks at least, indicating practical applicability of one-dimensional soil models in these sites to a certain extent. between surface and downhole) in EW and NS directions for earthquake records obtained there. At all the sites, peak frequencies of 2As/(Ab + Bb) are mostly compatible with those of observed spectra in lower frequency peaks at least, indicating practical applicability of one-dimensional soil models in these sites to a certain extent.
the correspondence in peak frequency is perfect in (a) and good in (b) despite the peak heights are distinctively higher in the former than the latter because of different contribution of radiation damping as mentioned later. The correspondence seems to be still tolerable in (c). However, it gets poorer in (d) and very poor in (e) and (f) because the spectrum 2A The reason may be given by examining the soil profiles. In (a) and (b), the V s -value at the base is much larger than the upper layers, and the depth of the downhole seismometer is not so deep from the boundary of clear V s -contrast indicated by a triangle. In (c) and (d), the V s -value at the base layer is not so much varied from that in the upper layer, and the seismometer depth is getting deeper from the boundary of V s -contrast again. In (e) and (f), though V s at the base layer is much greater than the upper layers, the depth of the seismometer is too deep from the boundary of clear V s -contrast to properly detect the response of the upper layers. This observation tells us the importance of appropriate in-advance planning in deploying a vertical array system considering site-specific soil profiles.
In order to clarify a basic principle of how peak frequencies in 2A Normalized frequency f /f 1 Surface layer Seismometers (a) 3-layers model , the depth of seismometer installation works as a virtual rigid boundary with no wave energy radiation into the underlying layer even if it is in the midst of a uniform layer according to a simple multireflection theory of SH-wave (e.g., Schnabel et al., 1972 [15], Kokusho 2017 [16]).
In (c), where the impedance ratios α 12 , α 23 are parametrically changing (keeping α 12 × α 23 = 0.16 constant) with a constant thickness ratio and 2A s /2A b are quite different in peak frequencies for α 23 ≈ 0.64 or larger, though they coincide to each other if α 23 ≈ 0.4 or smaller. This indicates that a sharp impedance contrast is preferable at the base boundary near the down-hole seismometer for better matching between the two transfer functions.
Consequently, in deploying vertical array systems, the downhole seismometers should be installed with enough care in a stiff base layer not too deep from the layer boundary. The clearer the impedance contrast at the base boundary to the overlying layer, the better matching in peak frequencies is possible. If the depth of the down-hole seismometer is twice or further deeper than the surface layer thickness, H 2 /H 1 > 2.0, from the boundary of major impedance contrast, virtual peaks with quite different frequencies from those in for vertical arrays unless the impedance contrast is clear enough. Thus, there should be an appropriate installation depth for the down-hole seismometer closely dependent on individual site conditions. It is never the deeper, the better. The installation depth of down-hole seismometers may sometimes be predetermined by the economy or other specifications. In the case of KiK-net, too, due care was necessary to tackle this problem, and only appropriate data have been chosen in this research for the amplification analysis as explained below.

Earthquake Records and Soil Profiles
KiK-net records corresponding to 8 earthquakes (EQ1 to EQ8) utilized in this paper are listed in Table 1. The moment magnitudes M w listed together with M J (Japan Meteorological Agency magnitude) are varied from 6.6 to 7.9. A total of 118 mainshocks (MS) records during the eight earthquakes observed in vertical array sites near focal areas were firstly addressed. As shown in Figure 6a, the PGA-values in the vertical axis are spanning from 0.4 to 2.5 g, while maximum accelerations at the base are 0.01-1.0 g, leading to approximately 2-16 times amplification in maximum horizontal acceleration in terms of 2A s /(A b + B b ) between surface and base. Among them, those with filled symbols have been selected for the amplification analysis later on.
The depths of the down-hole seismometers in the vertical arrays addressed here vary mostly from 100 m to 330 m. S-wave velocities in the top and bottom of all the vertical arrays are plotted in the vertical and horizontal axes, respectively, in Figure 6b. The bottom velocities V sb in many earthquake recording sites diverge from 400 m/s to 3000 m/s, though all of them are stiff enough to be considered as the engineering bedrock while some are approaching seismological bedrock corresponding to the earth crust of hard rock with a high density ρ ≥ 2.7 t/m 3 and high S-wave velocity V s ≥ 3000 m/s. In a good contrast, the majority of surface velocities V s in most of the sites are within a limited extent; about 100 to 400 m/s, around 200 m/s on average. Among the plots, those of filled symbols have been selected for the later analysis.
Meteorological Agency magnitude) are varied from 6.6 to 7.9. A total of 118 mainshocks (MS) records during the eight earthquakes observed in vertical array sites near focal areas were firstly addressed. As shown in Figure 6a, the PGA-values in the vertical axis are spanning from 0.4 to 2.5 g, while maximum accelerations at the base are 0.01-1.0 g, leading to approximately 2-16 times amplification in maximum horizontal acceleration in terms of 2As/(Ab + Bb) between surface and base. Among them, those with filled symbols have been selected for the amplification analysis later on. Fourier spectrum ratio 2As/(Ab + Bb) between surface and base is calculated using the surface and downhole records in each KiK-net vertical array site. A typical result in EW and NS direction for the mainshock and several aftershocks at a site during EQ3 is shown together with associated soil profile data in Figure 7a,b. A similar result at another site during EQ3 is shown in Figure 8a,b. The aftershocks (AS) were chosen among the KiK-net dataset to have distinctively low PGAs and hypocenters not far compared to the mainshock (MS) as indicated in Table 1. It is noted in Figures 7 and 8 that the spectrum peak frequencies are almost matching among multiple aftershocks in EW/NS directions in the lower frequency peaks, at least in these two sites. Also noticed is a clear difference in the amplitude and peak frequency of the spectrum ratios between the mainshock and aftershocks. It presumably reflects strain-dependent nonlinear soil properties during strong shaking, which appears to be more pronounced in higher-order peaks due to the larger involvement of shallower layers.   Fourier spectrum ratio 2A s /(A b + B b ) between surface and base is calculated using the surface and downhole records in each KiK-net vertical array site. A typical result in EW and NS direction for the mainshock and several aftershocks at a site during EQ3 is shown together with associated soil profile data in Figure 7a,b. A similar result at another site during EQ3 is shown in Figure 8a,b. The aftershocks (AS) were chosen among the KiK-net dataset to have distinctively low PGAs and hypocenters not far compared to the mainshock (MS) as indicated in Table 1. It is noted in Figures 7 and 8 that the spectrum peak frequencies are almost matching among multiple aftershocks in EW/NS directions in the lower frequency peaks, at least in these two sites. Also noticed is a clear difference in the amplitude and peak frequency of the spectrum ratios between the mainshock and aftershocks. It presumably reflects strain-dependent nonlinear soil properties during strong shaking, which appears to be more pronounced in higher-order peaks due to the larger involvement of shallower layers.

S-Wave Velocity Ratio versus Site Amplification
A flow chart is illustrated in Figure 9a to develop a reliable correlation between site amplification and S-wave velocity based on the vertical array database in this research. Totally 39 mainshock records have been selected out of 118 during eight earthquakes initially addressed as listed in Table 1. The following two criteria was employed in the selection to reduce various potential errors involved during the data reduction process.

S-Wave Velocity Ratio versus Site Amplification
A flow chart is illustrated in Figure 9a to develop a reliable correlation between site amplification and S-wave velocity based on the vertical array database in this research. Totally 39 mainshock records have been selected out of 118 during eight earthquakes initially addressed as listed in Table 1. The following two criteria was employed in the selection to reduce various potential errors involved during the data reduction process. (1) Similarity of transfer functions between 2As/(Ab + Bb) and 2As/2Ab calculated by the soil model in terms of peak frequencies as already discussed in Figure 4. Among the six charts in Figure 4; for example, (a) to (c) are acceptable because the first peak frequency (the lowest peak frequency) at least is coincidental in contrast to very displaced peaks in (d) to (f).
(2) The similarity of the first peak frequency in 2As/(Ab + Bb) between the observed spectrum ratio and the computed transfer function so that the soil model is reliable enough.
Here, the acceptable relative difference in peak frequencies was regulated within 1/1.5 to 1.5 between the two functions.
Then, fundamental mode frequencies f of the layered soil system were calculated by a 1/4-wavelength formula, Equation (1), in order to specify equivalent soil layers generating peak frequencies in the spectrum ratio, based on the soil profile together with  (1) Similarity of transfer functions between 2A s /(A b + B b ) and 2A s /2A b calculated by the soil model in terms of peak frequencies as already discussed in Figure 4. Among the six charts in Figure 4; for example, (a) to (c) are acceptable because the first peak frequency (the lowest peak frequency) at least is coincidental in contrast to very displaced peaks in (d) to (f). (2) The similarity of the first peak frequency in 2A s /(A b + B b ) between the observed spectrum ratio and the computed transfer function so that the soil model is reliable enough.
Here, the acceptable relative difference in peak frequencies was regulated within 1/1.5 to 1.5 between the two functions.
Then, fundamental mode frequencies f of the layered soil system were calculated by a 1/4-wavelength formula, Equation (1), in order to specify equivalent soil layers generating peak frequencies in the spectrum ratio, based on the soil profile together with the S-wave velocities of individual layers provided by the KiK-net database.
Here, as depicted in Figure 9b, H i and V si are the thickness and S-wave velocity of the i-th layer sequentially numbered from the top, and H i /V si is summed up layer by layer down to the base. This is an extension of the well-known 1/4 wavelength formula in a two-layer system for roughly determining predominant frequency in a multi-layer system based on the same principle of first resonance. The calculated frequency is compared one by one with the peak frequencies in the observed spectrum ratio such as in Figure 7a or Figure 8a to identify equivalent surface layers of thickness H = ∑ H i consisting of one or more layers generating the fundamental mode frequencies calculated by Equation (1) as tabulated in Figure 7b or Figure 8b.
In Figure 10a, the first peak frequencies f calculated by Equation (1)  f in Equation (1) is also compared with first peak frequencies of the transfer functions 2As/(Ab + Bb) and 2As/2Ab, respectively, computed by the 1D soil model based on the KiKnet soil profile data. Though some plots are dispersed beyond the lines f* = 0.8-1.2 f, the difference is well within 1/1.5-1.5 times as initially regulated. Then, the average S-wave velocity s V for the equivalent surface layer is determined by Equation (2)    In Figure 10b  In Figure 10b, V s -values (m/s) thus determined from the first peak frequencies are compared with V s30 (m/s), S-wave velocities averaged over the top 30 m from ground surface often used in current practice, which is calculated by the next equation Equation (3) where T 30 (s) is the time of S-wave propagation for the top 30 m based on V s -logging data.

It indicates that V s -value introduced in this research is positively and loosely correlated
The next step after determining the average S-wave velocity of the equivalent surface layer generating the peak frequency is how to evaluate 2A s /2A b from 2A s /(A b + B b ) as illustrated in the flowchart (Figure 9a) and described below.
(1) A theoretical transfer function, 2A s /(A b + B b ) is calculated based on a soil model in each site using the given soil/V s -profile with frequency-independent damping ratio, D, tentatively assumed as 2.5% throughout the layer.
is compared with the spectrum ratio of observed records, as exemplified in Figure 11. If a peak frequency in the transfer function can be found within the allowable difference in the spectrum ratio of observed motions, it is identified as the corresponding peak, and the damping ratio assumed as D = 2.5% previously is modified uniformly by D = Q 1 /Q 2 × 2.5% to have the identical peak value, where Q 1 is the peak amplitude of the transfer function, and Q 2 is that of spectrum ratio of the observed records. (3) Another theoretical transfer function 2A s /2A b is computed using the modified damping ratio D (averaged between EW and NS directions) based on the same soil model, and the peak value is read off as the peak amplification for the corresponding horizontal array as also indicated in Figure 11.      In previous research conducted by the present author (Kokusho and Sato 2008 [13]), wherein out of the same KiK-net vertical array data, those belonging to EQ3, EQ4, EQ 5 were analyzed, the correlation 2A a /2A b = 0.685 × V sb /V s + 0.175 (R 2 = 0.828) was obtained. Another previous research (Kokusho 2013 [14]), where the same KiK-net data of 8 earthquakes was utilized, 2A a /2A b = 0.626 × V sb /V s + 0.369 (R 2 = 0.792) was yielded. In these previous papers, unlike this paper the strict data selection with regard to peak frequency correspondence between 2A s /2A b and 2A s /(A b + B b ) using the criteria already mentioned, was not implemented. Furthermore, not only the first peak but also higher-order peaks were also analyzed in a similar manner. The correlation obtained here seems quite reasonably better than the previous ones with higher R 2 -values, though the previous correlations superposed in Figure 12a with thin lines are not so much different. Figure 11. Example how to determine peak value of 2As/2Ab from 2As/(Ab + Bb) compared with observed spectrum ratio.  ( ) In previous research conducted by the present author (Kokusho and Sato 2008 [13]), wherein out of the same KiK-net vertical array data, those belonging to EQ3, EQ4, EQ 5 were analyzed, the correlation ( ) In these previous papers, unlike this paper the strict data selection with regard to peak Linear approx. V s -ratio: In Figure 12b, the same peak values in 2A s /2A b are plotted versus velocity ratio, V sb /V s30 . It is approximated by the next equation with R 2 = 0.722, wherein the plots are evidently more scattered than in Figure 12a. This indicates that the S-wave velocity of the top 30-m V s30 may somehow evaluate the site amplification to a certain extent despite that the effect of site-specific soil-profiling is ignored.
Thus, the site amplification for zonation study relative to an outcropping base layer has been formulated for first spectrum peaks corresponding to the predominant frequency in situ, of many vertical array sites during eight strong earthquakes as either Equations (4) or (5). The former may be recommended because of the larger R 2 -value for higher accuracy in predicting the site amplification, though both formulas may be useful practically. It should be noted that they may be applicable to arbitrary engineering bedrock with various V sb and depth because they have been developed from soil profiles including a wide variety of base layers with V sb = 400 m/s to 3000 m/s as indicated in Figure 6b.
With respect to Equation (4), site amplification may be readily assessed by the aid of microtremor measurements, even if no information on soil profiles is available, as (1) Decide predominant frequency f of a site using H/V spectrum ratios calculated from Horizontal and Vertical microtremor records (Nakamura 1989 [17]). (2) Estimate the thickness H of soft soil or equivalent surface layer where the predominant frequency is exerted, which is sometimes possible based on geological maps or highdensity soil logging data available in city areas. (3) Then, the average S-wave velocity in the equivalent surface layer where the predominant frequency is exerted can be calculated by V s = 4H f . (4) Calculate the S-wave velocity ratio V sb /V s from V sb of a base layer, which can be normally given by geological database of the area. (5) Equation (4) using the V sb /V s -value will give the site amplification with respect to the common base layer. The relative amplification at two sites sharing the same base layer can be obtained by directly comparing them. If different base layers are chosen in two sites, the relative amplification can still be readily adjustable in Equation (4) by considering the difference in V sb . The same procedures can also be followed when Equation (5) for V sb /V s30 is used in place of Equation (4) for V sb /V s .

Soil Nonlinearity Effect
In order to investigate the effect of soil nonlinearity on the site amplification using actual earthquake data, aftershock records were compared with corresponding mainshock records in the same sites. For the comparison, 28 recording sites were selected, as summarized in Table 1 out of 39 where mainshocks were analyzed. Four aftershocks were chosen in individual sites (112 aftershock records in total) with their PGAs distinctively different, as small as 1% to 10~65% maximum to the mainshocks as summarized in Table 1  Bb) (vertical array) with open circles. It is noted that the absolute amplification values of 2As/2Ab are far smaller than 2As/(Ab + Bb) mainly because the former is strongly influenced by radiation damping by underlying layers, while the radiation damping is ineffective in the latter. It is also remarkable is that 2As/2Ab shows almost identical amplification between MS and AS, while 2As/(Ab + Bb) seems to be quite different, wherein the mainshock amplifications can be as low as half of the aftershocks. In Figure 13b, the same peak spectrum amplifications, 2As/(Ab + Bb) and 2As/2Ab, are plotted versus s V -ratio ( sb s V V ) for mainshocks and aftershocks of the 8 earthquakes with filled and open symbols. For 2As/2Ab, a clear correlation can be recognized between the peak amplification and s V -ratio, which can be approximated by Equation (4), again. It should also be noted that the difference in peak amplifications between the mainshock and corresponding aftershocks are really marginal, indicating a negligible effect of soilnonlinearity. In contrast, the plots for 2As/(Ab + Bb) are not so well correlated to s V -ratio.
Furthermore, the peak amplifications of 2As/(Ab + Bb) for the mainshocks (filled triangles) are evidently smaller than the corresponding aftershocks (open triangles) in most sites. This indicates that the soil nonlinearity effect is evidently seen in vertical array earthquake records but far less pronounced in surface array records. The soil nonlinearity effect was also clearly recognized independently by Régnier et al., (2013), where the borehole site response corresponding to 2As/(Ab + Bb) was investigated using KiK-net vertical array records, though outcrop site response 2As/2Ab was not addressed in the same study. This indicates that the soil nonlinearity effect is evidently seen in vertical array earthquake records but far less pronounced in surface array records. The soil nonlinearity effect was also clearly recognized independently by Régnier et al. (2013), where the borehole site response corresponding to 2A s /(A b + B b ) was investigated using KiK-net vertical array records, though outcrop site response 2A s /2A b was not addressed in the same study.
where ( )  In taking into account the effect of strain-dependent soil properties on the amplification, the shear modulus degradation G/G0 from small-strain shear modulus G0, and the damping ratio D in the surface layer are parametrically changed; G/G0 = 1.0, 0.65, 0.25 and D1 = 2.5, 5 and 15%, for strain level of 5 × 10 −6 , 1 × 10 −4 , 1 × 10 −3 , respectively, assuming a typical degradation curve for sand (Seed & Idriss 1970 [18]) as indicated in Figure 14b, while in the base layer D2 = 0. The calculated results of 2As/(Ab + Bb) and 2As/2Ab are shown in Figure 15a,b, respectively. Here, the amplification between surface and base is shown versus the normalized frequency, f/f1, where f1 = first mode frequency of the surface layer for small strain properties (G/G0 = 1.0). In taking into account the effect of strain-dependent soil properties on the amplification, the shear modulus degradation G/G 0 from small-strain shear modulus G 0 , and the damping ratio D in the surface layer are parametrically changed; G/G 0 = 1.0, 0.65, 0.25 and D 1 = 2.5, 5 and 15%, for strain level of 5 × 10 −6 , 1 × 10 −4 , 1 × 10 −3 , respectively, assuming a typical degradation curve for sand (Seed & Idriss 1970 [18]) as indicated in Figure 14b Figure 15a,b, respectively. Here, the amplification between surface and base is shown versus the normalized frequency, f /f 1 , where f 1 = first mode frequency of the surface layer for small strain properties (G/G 0 = 1.0).
Obviously, nonlinear soil properties have a great effect on peak frequency and amplification. In particular, the peak frequency becomes evidently lower with modulus degradation both in 2A s /(A b + B b ) and 2A s /2A b . However, with regard to peak amplification, it is noted that the soil nonlinearity has a much smaller effect in 2A s /2A b than in 2A s /(A b + B b ) for the first peak in particular. This is because radiation damping associated with impedance ratio α = ρ 1 V s1 /ρ 2 V s2 tends to dominate 2A s /2A b overwhelming the soil nonlinearity effect. In Equation (7), the amplification is expressed as |2A s /2A b | ≈ |1/α * | when k 1 * ≈ π/2H at the first peak, while no effect of impedance ratio α * is involved in 2A s /(A b + B b ) as indicated in Equation (6). Under the paramount effect of radiation damping, the difference in the amplification due to strain-dependent properties becomes obscure in 2A s /2A b . Furthermore, the impedance ratio α = ρ 1 V s1 /ρ 2 V s2 , which becomes smaller with degraded shear modulus in the surface layer, tends to give larger amplification, compensating for the increasing effect of internal soil damping during strong earthquakes. Thus, the difference in soil properties between mainshock and aftershocks tends to have less influence on the amplification in 2A s /2A b than in 2A s /(A b + B b ), as demonstrated in Figure 13 by actual earthquake records. Figure 14. Two-layers model calculating soil nonlinearity effect on site amplification (a), and Strain dependent variations of normalized shear modulus G/G0 and damping ratio D at three strain levels (b).
In taking into account the effect of strain-dependent soil properties on the amplification, the shear modulus degradation G/G0 from small-strain shear modulus G0, and the damping ratio D in the surface layer are parametrically changed; G/G0 = 1.0, 0.65, 0.25 and D1 = 2.5, 5 and 15%, for strain level of 5 × 10 −6 , 1 × 10 −4 , 1 × 10 −3 , respectively, assuming a typical degradation curve for sand (Seed & Idriss 1970 [18]) as indicated in Figure 14b, while in the base layer D2 = 0. The calculated results of 2As/(Ab + Bb) and 2As/2Ab are shown in Figure 15a,b, respectively. Here, the amplification between surface and base is shown versus the normalized frequency, f/f1, where f1 = first mode frequency of the surface layer for small strain properties (G/G0 = 1.0). Obviously, nonlinear soil properties have a great effect on peak frequency and amplification. In particular, the peak frequency becomes evidently lower with modulus degradation both in 2As/(Ab + Bb) and 2As/2Ab. However, with regard to peak amplification, Summarizing the above discussions, there are two definitions of site amplification; the horizontal array and vertical array, as shown here again in Figure 16a, where the horizontal array is directly applicable to the seismic zonation. Figure 16b conceptually shows the accelerations at the stiff base, the point B (outcropping) or B (downhole) taken in the horizontal axis versus those at the soil surface, the point A, in the vertical axis. The acceleration on the soft soil surface shows a strong nonlinearity with respect to increasing downhole-based acceleration at B due to nonlinear soil properties as schematically illustrated with the dashed shaded belt. In the horizontal array, however, the surface acceleration is remarkably linear with the outcropping base acceleration at B, as illustrated with the solid shaded belt in the same diagram. It has been actually evidenced in Figure 13a that there is no significant amplification difference between mainshock and corresponding aftershocks at many recording sites. This is because the base acceleration in the horizontal axis defined at the outcropping base B tends to be larger than that at the downhole base B where the effect of downward waves from the overburdened soils cannot be ignored. it is noted that the soil nonlinearity has a much smaller effect in 2As/2Ab than in 2As/(Ab + Bb) for the first peak in particular. This is because radiation damping associated with impedance ratio tends to dominate 2As/2Ab overwhelming the soil nonlinearity effect. In Equation (7), the amplification is expressed as at the first peak, while no effect of impedance ratio * α is involved in 2As/(Ab + Bb) as indicated in Equation (6). Under the paramount effect of radiation damping, the difference in the amplification due to strain-dependent properties becomes obscure in 2As/2Ab. Furthermore, the impedance ratio , which becomes smaller with degraded shear modulus in the surface layer, tends to give larger amplification, compensating for the increasing effect of internal soil damping during strong earthquakes. Thus, the difference in soil properties between mainshock and aftershocks tends to have less influence on the amplification in 2As/2Ab than in 2As/(Ab + Bb), as demonstrated in Figure 13 by actual earthquake records.
Summarizing the above discussions, there are two definitions of site amplification; the horizontal array and vertical array, as shown here again in Figure 16a, where the horizontal array is directly applicable to the seismic zonation. Figure 16b conceptually shows the accelerations at the stiff base, the point B (outcropping) or B′ (downhole) taken in the horizontal axis versus those at the soil surface, the point A, in the vertical axis. The acceleration on the soft soil surface shows a strong nonlinearity with respect to increasing downhole-based acceleration at B′ due to nonlinear soil properties as schematically illustrated with the dashed shaded belt. In the horizontal array, however, the surface acceleration is remarkably linear with the outcropping base acceleration at B, as illustrated with the solid shaded belt in the same diagram. It has been actually evidenced in Figure  13a that there is no significant amplification difference between mainshock and corresponding aftershocks at many recording sites. This is because the base acceleration in the horizontal axis defined at the outcropping base B tends to be larger than that at the downhole base B′ where the effect of downward waves from the overburdened soils cannot be ignored. A thick solid curve overlaid in Figure 16b represents a similar relationship of PGAvalues at soft soil sites and stiff rock sites proposed in the United States (Idriss 1990 [19]). It was estimated in the framework of the horizontal array from actual smaller earthquake records of similar epicenter distances during the 1985 Mexican earthquake and 1989 Loma A thick solid curve overlaid in Figure 16b represents a similar relationship of PGAvalues at soft soil sites and stiff rock sites proposed in the United States (Idriss 1990 [19]). It was estimated in the framework of the horizontal array from actual smaller earthquake records of similar epicenter distances during the 1985 Mexican earthquake and 1989 Loma Prieta earthquake and from numerical analyses for strong motions as well. The curve is obviously nonlinear, reflecting the soil nonlinearity, and the amplification of soil sites to rock sits to be lower than unity for PGA larger than around 0.4 g, though the PGA and the spectrum peak amplification may possibly be somewhat different in reflecting the soil nonlinearity effect.
The above findings based on the actual records as well as the simple model analysis may indicate that soil nonlinearity is not so pronounced in the horizontal array, unlike the vertical array in terms of spectrum peak amplification, including strong earthquakes. The mainshock records utilized here were with PGA around 0.6 g and even as high as 2.5 g. As for the extreme case of PGA = 2.5 g, no clear difference in the peak amplification of 2A s /2A b was observed by one of the plots in Figure 13. The near-surface soil profile there was not so stiff, comprising 2 m thick layer of V s = 150 m/s underlain by 18 m thick layer of V s = 430 m/s, followed by stiff strata of V s = 1000 m/s or greater down to 100 m deep base layer of V s = 1500 m/s. Nevertheless, it should be remarked here that no sites addressed here experienced soil liquefaction during the earthquakes. It is obvious that when soil shear moduli are lost substantially due to liquefaction, the amplification will be way down from unity (deamplification) because most wave energy cannot arrive at the soil surface. Such cases were demonstrated by earthquake records during the 1964 Niigata earthquake (Kanai 1966 [10]), 1987 Imperial Valley earthquake (Adalier et al., 1997 [20]), 1995 Kobe earthquake (Kokusho & Matsumoto 1998 [11]), and many others, and may be interpreted as a kind of base-isolation wherein wave energy arriving at soil surface tends to decrease dramatically (Kokusho 2014 [21]). This indicates that, at some points when shear stiffness is nearly lost due to pore-pressure buildup making S-wave untransmissible, the downward waves in the overburden soils diminish accordingly, and the downhole accelerations A b + B b at B' ultimately become identical with the outcropping acceleration 2A b at B, merging the two curves in the horizontal and vertical arrays as illustrated in the diagram. Up to that point, the site amplification defined by the horizontal array 2A s /2A b may possibly occur more linearly than normally anticipated.

Conclusions
Site amplification formula between the soil surface and engineering bedrock for seismic zoning has been developed by utilizing a number of KiK-net vertical array records during recent eight strong earthquakes. The first peak amplifications of horizontal array transfer functions 2A s /2A b have been determined from those of vertical arrays 2A s /(A b + B b ) compatible with the observations and correlated with S-wave velocity ratios between the base layers and equivalent surface layers, yielding the following major outcomes.
(1) It was found that the transfer function and its peaks for the vertical array 2A s /(A b + B b ) may sometimes differ considerably from those of horizontal array 2A s /2A b needed in seismic zonation, as a significant potential drawback of the vertical array. It is attributed to the improper installation depth of the downhole seismometer with respect to engineering bedrock, indicating a caveat "never be the deeper the better". Criteria in deploying the downhole seismometer in proper depth with regard to soil profiles has been provided based on simple analyses. (2) Average S-wave velocity V s for an equivalent surface layer, generating a peak in spectrum amplification, has been incorporated here to evaluate site amplification. The first peak amplifications of 2A s /2A b determined from vertical array observation data are correlated well with V sb /V s (V sb = S-wave velocity at bedrock) in Equation (4) with the determination coefficient R 2 = 0.876 for the mainshocks recorded at 39 sites during eight strong earthquakes.
(3) The velocity V s of the equivalent surface layer, if V s -logging data is unavailable, may be readily determined from predominant frequency f of a site using H/V spectrum ratio in microtremor measurement together with the thickness of soft soil or Holocene layer H by V s = 4H f . (4) The same peak amplifications 2A s /2A b are also correlated positively with V sb /V s30 (V s30 = average velocity for top 30 m often used in the current seismic zonation practice) in Equation (5) with a slightly smaller coefficient of R 2 = 0.722. Both Equations (4) and (5) may be applicable in evaluating site amplifications relative to base layers with arbitrary V s values V sb . The relative surface amplification between two different sites can also be readily determined from the correlations.
Furthermore, the effect of nonlinear soil properties has been studied by comparing peak amplifications for mainshock and small shocks recorded at the same sites, revealing the following; (2) According to a basic study on a simple two-layer equivalent linear system, the soil nonlinearity effect becomes minor on the peak amplification in 2A s /2A b mainly because of the paramount effect of radiation damping in the bedrock in contrast to 2A s /(A b +B b ) where no effect of radiation damping is involved. Thus, site amplification in seismic zonation may generally be considered in the realm of linear soil properties.
(3) However, if extensive soil liquefaction occurs, the base-isolation effect due to a drastic decrease of soil modulus may completely change this trend. Up to that point, the site amplification defined by the horizontal array 2A s /2A b to be used in seismic zonation may occur more linearly than generally anticipated.
Author Contributions: Conceptualization, methodology, data selection, investigation, analyses, writing and editing, visualization, supervision were carried out by T.K. Data selection and analyses were supported by T.I. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.
Data Availability Statement: All earthquake data presented in this study are available for free for research purposes by accessing the website (https://www.kyoshin.bosai.go.jp/kyoshin/, accessed on 6 December 2021) of NIED, Tsukuba, Japan.