Unified approach for evaluation of horizontal site amplification factors with special reference to history of studies on the effects of surface geology on seismic motion

Following the 1923 Kanto earthquake in Japan, Japanese researchers noticed the strong effects of surface geology on seismic motion (ESG) and began to investigate these effects to quantify the site amplification factors (SAFs) associated with soft surface sediments. On the other side of the Pacific Ocean, ESG received limited attention until the 1985 Michoacan, Mexico earthquake revealed significant long-period amplification inside Mexico City that manifested as the source of devastating damage to high-rise buildings. Since then, seismologists and earthquake engineers have performed a lot of studies on various ESG issues worldwide. We have not yet reached common conclusions on how to quantitatively predict SAFs over a broad frequency band of engineering interest, 0.1 to 20 Hz, for moderate to strong input from different types of earthquakes in different tectonic settings. However, we found here several basic guidelines useful for successfully modeling ground motions as a common approach to ESG studies. First, in this letter, we briefly review our history of understanding ESG, which is closely related to the key settings required for reliable quantifications of SAFs, and then introduce various emerging techniques for broadband quantitative evaluations of SAFs based on the vast amount of observed ground motions primarily from dense Japanese strong-motion networks. Based on the findings of our investigation and the physical relationships behind the parameters, the authors would like to recommend that researchers on ESG and related topics would refer to the five basic guidelines proposed in the conclusions for the successful implementation of techniques to delineate SAFs in a specific region of interest, such as the use of Fourier spectra instead of response spectra. We have started applying the proposed techniques to regions outside Japan. The implementation of the statistical validation exercises will follow.


Introduction
It is now apparent to every researcher in the field of seismology or earthquake engineering that the characteristics of the observed seismic motions strongly depend on the underlying geophysical structures or surrounding topographies adjacent to the target site.However, these mechanisms might not have been so obvious 100 years ago.At the time, the physical mechanisms underlying earthquake occurrence were poorly understood and even a magnitude scale had not been developed yet.Therefore, identifying a strong correlation of areas with heavily damaged buildings and houses in Tokyo and Yokohama with their surface geology represented a considerable discovery (e.g., Kitazawa 1926;Omote 1949) during the 1923 Kanto earthquake of M7.9. Figure 1 shows a map of collapsed houses (with blue and red circles) in downtown Tokyo according to the investigation by Kitazawa (1926).Collapsed houses were concentrated along the buried valleys of the Kanda River and Sumida River as white zones (classified as outcrop areas of alluvium).In contrast, very few collapsed houses were found in the yellow hilly zones (classified as outcrop areas of diluvium).Then, Sezawa and his colleagues (e.g., Sezawa 1930;Sezawa and Kanai 1932) developed theoretical tools for predicting the characteristics of seismic motions in the form of the site amplification factor (SAF) in the frequency domain, and Ishimoto (1931) reported observational evidence of the effects of surface geology on seismic motion (ESG).In this early stage of ESG studies, these were considered to be the effects of "surface geology".A recent review of the site characterization of ground motions by Trifunac (2016) also covered the early stages of ESG studies.Please note that before the intensive investigations on ESG performed after the devastating damage in Kanto, Japan, some topographic and soft-soil effects were noticed in Italy after the 1908 Messina earthquake, which were reflected in the building code in 1909 in Italy.Although other regions may have also reported similar findings in the early stage of ESG, they are beyond the scope of this article.
Later, these discoveries were effectively utilized by Takahashi and Hirano (1941), who deconvolved the theoretical one-dimensional (1D) horizontal site amplification factors (HSAF) from the observed ground motions at two sites to obtain common engineering bedrock motions for the first time.Figure 2 shows a simulation of their theoretical model, calculated HSAF, and a comparison of the deconvolved motions, along with the observed motions at Ueno on a hill (river terrace) and Shinobugaoka inside a buried alluvial valley.This first attempt presents three crucial elements of ESG studies: (1) observation of ground motions; (2) development of a theoretical method for the evaluation of ESG; and (3) quantification Fig. 1 Map of the areas with heavily damaged buildings and houses in Tokyo with their surface geology (after Kitazawa 1926).Outcrop regions of the diluvium formation are colored yellow, buried regions with thin alluvium pink and red, and thick alluvium regions white of the ground structure used in theory, as pointed out in the review of ESG studies by Kawase (1993Kawase ( , 2003)).
First, this letter briefly reviews our understanding of ESG as scientific knowledge on the critical settings for more optimized prediction of ground motions.We begin by explaining the reason why we refer to the discipline of the site effects of seismic ground motions as ESG, as referenced in the IASPEI/IAEE Joint Working Group of ESG, even though the site effects of seismic ground motions are not a direct consequence of "surface geology" only.We then introduce our recent achievements on ESG, which are primarily based on dense strong-motion network data in Japan.They will provide a reference for future research to elucidate the nature of ESG.
The IASPEI/IAEE Joint Working Group on Effects of Surface Geology on Seismic Motion (JWG-ESG) was started in 1985, led by Brian Tucker of Geo Hazard International (California Division Mines and Geology at that time) and the late Bill Iwan of California Institute of Technology.The amplifications of seismic motions at soft-sediment sites relative to stiff soil or rock sites in Japan from the beginning of the contemporary building code enforced in 1950 based on the experience of the 1923 Kanto earthquake, as mentioned above, are commonly considered.However, on the other side of the Pacific Ocean, awareness of the systematic influence of ESG was low until the 1985 Michoacan, Mexico earthquake, which revealed remarkably large long-period site amplifications inside Mexico City (e.g., Anderson et al. 1986).Therefore, 1985 was the right time to establish JWG-ESG, although this was merely a coincidence.
The first official meetings and workshops were held in 1987 during the IASPEI meeting in Vancouver, Canada.At the time, the most important and common concern among seismologists and earthquake engineers in strongmotion studies was the expression of the site term in the attenuation function of strong motions.Please note that the name as the attenuation function or attenuation formula had been used to refer to the prediction equation of ground motions until approximately 2000, but that now the ground motion prediction equation (GMPE) or ground motion modeling (GMM) is used instead.Because of that common interest in the prediction equations, they wanted to clearly define the characteristic of "surface geology" at a site that might precisely predict the site term in the attenuation function.According to the results of 35 years of investigation, the characteristics of seismic ground motions vary based on the effects of geophysical structures around the site through various types of wave propagation inside the medium and not based on the direct effects of "surface geology" (without depth information) at the site.The different results in different frequency ranges may arise from different parts of the same geophysical structure.Therefore, the word "ESG" in the context of this logical consequence refers to any site effects on seismic motions, as in the activities of JWG-ESG including the ESG symposia that started in 1992.In this article, ESG refers to the SAF in the Fourier amplitude spectra (FAS) of seismic ground motions.We present our various SAF prediction schemes in FAS as a unified approach.However, we never intend to provide an improved technique for conventional site-effect prediction methods, such as new proxies in GMPEs.This article is an extensively revised version of the keynote lecture presented during the 6th ESG symposium (ESG6) in 2021 (Kawase 2021).

Site classification or site proxy approaches
Approaches based on site classification or indices, called site proxies, are still used in most site-effect evaluation schemes in the framework of GMPEs.Initially, as mentioned in the Introduction, surface geology or soil-type classification, such as soil and rock, was used to describe the site term in GMPEs, based on the ergodic assumption in the field of probabilistic seismic hazard assessment.This setting means that all sites with the same site category were assumed to have a common site term, which is far from reality (Trifunac 2016).Aki (1988) reviewed local site effects and reported that the broad classification (e.g., soil versus rock) of site conditions fails to capture the essential factor that controls the site effect as frequency-dependent characteristics.Other site proxies were then searched to better characterize the site term in the GMPEs.The time-averaged S-wave velocity (Vs) of the top 30 m (Vs30) has been the most widely used single-site proxy since its introduction by Borcherdt (1994) and Boore et al. (1997).According to Strasser et al. (2009), the overall variability in GMPEs with respect to the observed values has not been significantly reduced for more than three decades, as reproduced in Fig. 3.The effectiveness of Vs30 as a site-specific proxy for ESG was discussed during the 4th Symposium on the Effects of Surface Geology on Seismic Motion held at the University of California at Santa Barbara in 2011.At the symposium, Abrahamson (2011) suggested that Vs30 was a suitable site parameter in his GMPE, while Zhao (2011) investigated the modeling of site effects using the fundamental period of shallow layers (T S ) in addition to Vs30 (see some follow-up work in Zhao and Xu 2013).
Despite the increasing numbers of observed data and explanatory variables, one reason for the persistent variability of approximately 0.3 on the log10 scale (~ 0.7 in the natural log scale) in Fig. 3 is coming from the fact that almost all the GMPEs use a strength index such as peak ground acceleration (PGA), peak ground velocity (PGV), or response spectra as a target value of ground motions.These strength indices are complex functions of the SAF in FAS, which is the only physical value of the site effect that is directly influenced by the geophysical structure around the site.Another major reason for the failure to reduce the variability, even as a single-station sigma (e.g., Atkinson 2006), is the insufficient use of information for the geophysical structure, which is discussed later in detail.
Another reason for the persistent variability seen in Fig. 3 is the so-called ergodic assumption.Under the ergodic assumption of site amplification in GMPEs, a failure of prediction with significant bias at one site will be treated as the variability as a whole.In other words, under the ergodic assumption, a real site factor at one specific location is considered to have inevitable variability within the same site classification or the same Vs30 group.Therefore, an emerging field of research in the GMPE community is the use of site-specific GMPEs in response spectra or Fourier spectra (e.g., Abrahamson et al. 2019;Pinilla-Ramos et al. 2022).In this regard, our unified approach proposed here uses SAF of FAS derived from the generalized inversion technique (GIT) or observed spectral ratios as the site-specific term.Therefore, the overall variability of predictions should be smaller than that of the conventional GMPE predictions with the ergodic assumption.
As mentioned in the previous section, the use of the predominant period for site classification in Japan to allow for strong-motion evaluation was based on the developments by Sezawa, Kanai, and their colleagues (e.g., Sezawa and Kanai 1932;Kanai 1952Kanai , 1961;;Kanai et al. 1954) following the implementation of the soil-type spectral correction factor in the building code enforced in 1950.Initially, the soil-type classification was based on the surface geology of sediments (e.g., hard, medium, and soft sites).However, correspondence of the predominant period from the microtremors, Tp (e.g., 0.2 s > Tp, 0.2 < Tp < 0.6 s, and 0.6 < Tp for hard, medium, and soft sites, respectively), as determined by the zero-crossing method, was subsequently proposed (e.g., Kanai and Tanaka 1961) and has since been implemented in daily building design activities in Japan.After the proposal of Nakamura (1989), together with digital signal processing instead of the zero-crossing method, soil classification for building design in Japan has been primarily based on the predominant period in the microtremor horizontal-tovertical spectral ratio (MHVR) of the Fourier spectra.We should understand that at that time, using the predominant period as a representative value for soil classification in Japan was a direct consequence of the pioneering studies by Sezawa, Kanai, and their colleagues mentioned above.
Therefore, before the implementation and extensive application of Vs30 in the United States, it was natural in Japan to use a predominant period for site classifications, as evidenced by Fukushima and Tanaka (1990), who determined the relative amplification ratios for four categories of soil types in their GMPE of peak ground accelerations (PGAs).Later, Zhao et al. (2006) integrated their soil-type classification into their GMPEs based on the Tp (SC I to SC IV with Tp < 0.2 s, 0.2 < Tp < 0.4 s, 0.4 < Tp < 0.6 s, and 0.6 < Tp), derived from the earthquake horizontal-to-vertical spectral ratio (EHVR) of response spectra.Although the proposed correspondence in these studies is based on the conventional ergodic assumption, the inclusion of the Tp for the site classification of GMPEs has been proven effective from the inception of development.
Recently, Cadet et al. (2012b) introduced the fundamental predominant frequency (f0) to the site proxy, which was defined by Cadet et al. (2012a).In both studies, the authors used KiK-net operated by National Research Institute for Earth Science and Disaster Resilience (NIED 2019) as their data source.The integration of f0 with other proxies was also proposed by Castellaro and Mulargia (2014) primarily based on the theoretical consideration, by Hassani and Atkinson (2016) for strong motions in the eastern United States, and by Derras et al. (2017) for KiK-net.In these examples, explanatory variables such as the depth to the Vs = 800 m/s or larger layer (Z800) and topographic slope were introduced in addition to Vs30 and f0.Derras et al. (2017) identified Vs30 as the optimal single proxy for short periods (T < 0.6 s), whereas the performance of f0 and H800 was more optimized at relatively long periods.
More recently, Bergamo et al. (2021) compared the performance of various combinations of site proxies using neural network analyses to predict site effects from strong motion data from Swiss and Japanese (KiK-net) stations.Their targeted site effects were bedrock-impedance and borehole-effect corrected surface-to-borehole spectral ratios (SBRs) from KiK-net in Japan and the strong motion network in Switzerland.They used the quarter-wavelength method of Poggi et al. (2011) for the former (originally proposed by Boore and Joyner 1997) and the technique of Cadet et al. (2012a) for the latter as a semi-empirical method.They found that Vs30 was acceptable as a single-site proxy within the frequency range of 1.67 to 6.67 Hz, which is consistent with the report by Kawase and Matsuo (2004) based on the GIT for K-NET and KiK-net data.Kawase and Matsuo (2004) found that Vs10 and Vs20 outperformed Vs30 in that frequency range, while Bergamo et al. (2021) found that Vs20 and Vs10 were more optimized only at 5 Hz and 6.67 Hz, respectively, for the KiK-net data.Please note that there may be no clear physical reasons that a certain Vsxx (xx = 10, 20, 30) value would have a specific frequency range as a better proxy, as reported in these papers.These findings that a shallower Vsxx is better for the high-frequency proxy represent phenomenological findings; therefore, additional assumptions are required on the average velocity profile to interpret their physical mechanisms.According to Bergamo et al. (2021), the quarter-wave-length averaged velocity (VsQWL) as a function of frequency and its impedance contrast (IcQWL) was the optimal combination for a multiple site proxy.This is reasonable because the predominant peak frequency in 1D sedimentary layers is physically determined by the quarter-wavelength theory and the peak amplitude by the impedance contrast at the boundary of the corresponding depth.
Immediately after the proposal by Bergamo et al. (2021), Zhu et al. (2023) proposed a new effective use of various combinations of site proxies, including f0, to quantitatively reproduce the SAF at each site as a sitespecific value based on the SAFs obtained by Nakano et al. (2015).
Including the observed fundamental frequency f0 as a site proxy for GMPE is not considered appropriate as a general-purpose SAF prediction scheme because including the extracted value from the observed earthquake spectra at a targeted site leads to a circular discussion.In this circular case, f0 is assumed to be obtained from earthquake data such as SBR or EHVR.As long as f0 is estimated by independent information, such as a theoretical response of the S-wave velocity profile or an observed MHVR, then the inclusion of f0 is not exactly circular.Actually, Bergamo et al. (2021) showed the effectiveness of f0 as a single-site proxy in the middle-frequency range (1.67 to 5 Hz).In other words, the observed f0 from seismic motions cannot be an independent variable from the target site amplification for prediction.The same discussion is presented in a recent paper for nonergodic evaluation of the site term in a GMPE based on MHVR (Pinilla-Ramos et al. 2022).The use of observed ground motions to extract the SAF at that site is an optimized method of comprehensively utilizing the information of the observed data, as later shown in the next chapter.However, such a method cannot be used for sites without observation.

Empirical approach to determine the effects of surface geology on seismic motion
Two major approaches are used as empirical methods for ESG.One is based on the spectral ratio between the borehole and surface sensors (e.g., Kitagawa et al. 1988;Archuleta et al. 1992;Kobayashi et al. 1992;Satoh et al. 1995a;Steidl et al. 1996;Bonilla et al. 2002) as the SBR.The hypocentral distance in this configuration is virtually the same for both the target (surface) and reference (borehole) sites.Another empirical method is based on the spectral ratio of a soft sedimentary site to a nearby site on the rock, sometimes called the standard spectral ratio or SSR (e.g., Borcherdt 1970;Lermo and Chávez-García 1993;Field and Jacob 1995;Steidl et al. 1996;Bonilla et al. 1997).The effect of the distance between the two sites can be corrected using the attenuation function.The empirical method by SBR has difficulty eliminating the effects of reflected waves contained in the observed borehole seismograms (e.g., Steidl et al. 1996;Tao and Rathje 2020).Therefore, a good velocity structure with Q values is required to delineate the effects of reflected waves, as reported by Satoh et al. (1995a).Cadet et al. (2012a) proposed a simple semi-empirical method for correcting the differences in outcrop (surface-to-input) and within-layer (surface-to-borehole) site effects.However, they neglected site amplification from layers below the borehole station.Zhu et al. (2020) used SBRs for selected KiK-net sites with a shear-wave velocity above 800 m/s at the borehole sensor layer.Conversely, the empirical method that uses the data observed at a rock site as input into the denominator often has difficulty identifying a nearby site where the intact rock is exposed.Once such a rock site is identified, whether the site can be viewed as a site without soil amplification must be verified (e.g., Steidl et al. 1996;Cadet et al. 2012a).This issue is discussed in the Discussion section.
The GIT is a natural extension of these spectral ratio techniques with a considerably larger observed dataset (e.g., Andrews 1986;Iwata and Irikura 1988;Castro et al. 1988;Kato et al. 1992;Field and Jacob 1995;Harmsen 1997).The technique was introduced to determine, on a stable basis, the SAF at each site relative to an appropriately selected reference site or sites, which is contingent on numerous earthquakes observed in multiple sites, together with the source and path spectral terms.Owing to the dense strong-motion observation networks of K-NET, KiK-net, and the JMA seismic intensity (Shindokei) network in Japan, Kawase and Matsuo (2004) separated the basic properties of the source, propagation path, and site amplification from the observed strong-motion data.A subsequent study by Kawase (2006) focused on the nonlinear site amplification effect in records with a high PGA.Using the same strategy, Nakano et al. (2015) conducted spectral inversion analyses for Fourier spectra of more than 77,000 source-station pairs to determine the stress drop dependence on magnitude and depth.
A significant feature of the spectral inversions common to these studies is that they determined the SAF with respect to the outcrop equivalent of the seismological bedrock at a reference site, with an estimated S-wave velocity of 3,450 m/s.They obtained spectra by eliminating the effect of weathered rock layers from the surface records observed at the reference site before GIT.Such a clear definition of the reference site with an S-wave velocity higher than 3.0 km/s is indispensable to obtain stable SAF without trade-offs with the source spectra.Figure 4 shows examples of the separated horizontal site amplification factor (HSAF) and vertical site amplification factor (VSAF) at two sites in solid lines from a more recent study by Nakano et al. (2019) and Nakano (2021).The VSAF shown here is the averaged amplification factor of the vertical spectrum on the surface relative to the vertical bedrock spectrum, not the VSAF directly derived from the GIT.The latter is actually VSAF*V B H B R, as the relative amplification of the vertical spectrum on the surface to the horizontal bedrock spectrum, H B .V B H B R is the ratio between the vertical and horizontal components of the seismological bedrock spectra.VSAF*V B H B R is later represented as VSAF* to distinguish it from the VSAF relative to the vertical bedrock spectrum, V B .Significantly different characteristics were observed between the two sites.The dotted lines in the figure represent the theoretical HSAF and VSAF values at these sites, which are mentioned in the subsequent section.HSAF and VSAF are obtained only from the main S-wave part with a relatively short duration of motion (5 s for 4.5 < M JMA ≤ 6; 10 s for 6 < M JMA ≤ 7; 15 s for 7 < M JMA ≤ 8) and later dubbed sHSAF and sVSAF for S-wave portion.VSAF* is a significantly smoother function of frequency with a smaller amplitude than that of HSAF.Ito et al. (2020) proposed their averaged characteristics as a vertical amplitude correction function (VACF) to obtain the HSAF directly from the observed eHVSR.
Figure 5 shows distribution maps of HSAF at two frequencies (0.5 Hz and 1 Hz) throughout Japan interpolated by GMT (Smith and Wessel 1990).Amplification is prominently observed inside the major basins and major volcanic calderas in Kyushu (the southwestern Fig. 4 HSAF and VSAF at two sites obtained for the S-wave part by Nakano et al. (2019) and theoretical 1D site amplifications based on the unified velocity structures by Wakai et al. (2019) island) and Northern Honshu (the main island).Notably, the northeastern side from the Itoigawa-Shizuoka tectonic line has systematically higher amplification than the southwestern side.The HSAF pattern at 0.5 Hz is remarkably similar to that at 1 Hz, which is not a coincidence but a natural consequence of HSAF as a result of the 1D resonance from the seismological bedrock to the surface.
The separated HSAFs for borehole sensors are also quite useful for delineating the reflected-wave effects for borehole records and the amplification effects with different S-wave velocities of the layers where the borehole sensors are installed.Figure 6 compares (a) HSAF of the surface sensors relative to the seismological bedrock (HSAF_surface, which is virtually equal to HSAF); (b) HSAF of the borehole sensors relative to the seismological bedrock (HSAF_borehole); and (c) the ratio between them (SBR).Significant fluctuations can be observed in HSAF_borehole, although the maximum amplitude of the averaged HSAF_borehole is only 2.0 at approximately 1 Hz, and the amplitude in the higher frequency range is close to 1. Large fluctuations can be regarded as a function of frequency in all of these site factors, particularly in HSAF_surface and SBR.Therefore, deriving a simple method for empirically representing these site factors using a small number of site proxies, as discussed in the previous section, may be quite difficult, if not impossible.Despite such strong site-to-site variabilities, the averaged HSAF_surface and SBR were remarkably smooth functions of frequency, with peak amplifications of approximately 6.6 and 5.5, respectively, at a frequency of approximately 10 Hz.Notably, a common peak at approximately 0.3 Hz in HSAF_surface and HSAF_borehole is caused by the observed amplification at the reference site, YMGH01, below which we found a boundary at a depth of > 5 km with a small impedance contrast.We will consider the effect of this reference site amplification in our GIT in the near future.

One-dimensional amplifications
As mentioned in the Introduction, one-dimensional (1D) site amplification has been considered to be the most fundamental way of representing site amplification characteristics from as early as the 1940s.All textbooks on theoretical seismology, including Aki and Richards (1980) and Udias (1999), have sections that teach us how to represent the wave propagation phenomenon within a 1D layered structure.Observational pieces of evidence Fig. 5 Observed HSAF maps at frequencies of 0.5 Hz and 1.0 Hz relative to the seismological bedrock based on the GIT by Nakano (2021).Red dotted curves are approximate locations of the Itoigawa-Shizuoka tectonic line.Interpolation is done by GMT (Smith and Wessel 1990) have been provided in several past reviews (e.g., Aki 1988Aki , 1993;;Kawase 1993Kawase , 2003)).In addition, the author needs to mention that the fundamental validity of the 1D approach for HSAF evaluation was proven by the blind prediction exercise during the first ESG symposium for Ashigara Valley (Kudo 1992) as well as the blind prediction during the sixth ESG symposium for Kumamoto (Matsushima et al. 2021(Matsushima et al. , 2022)).Figure 7 compares the predicted results in the pseudo-velocity response spectra (Midorikawa 1992;Kudo and Sawada 1992) as reported during the first ESG symposium.At the time, the evaluation of the submitted predictions sounded somewhat negative, emphasizing the residuals from the observations and the variations among the submitted predictions.However, considering the inherent variations for multiple observations at the same site, which were lacking at the time, but are now known to represent less than one-half or twice of the average (as shown later in Fig. 14), their overall matching with the observations was statistically satisfactory.Also, the variability among predictions in Fig. 7 was reasonable because it is less than the inherent variations for multiple observations (please note that the variations in Fig. 7 among predictions are quartiles, not ± one standard deviation).
More recent simultaneous comparisons from different groups of researchers for linear and nonlinear 1D evaluation of surface-to-borehole spectral ratios (SBR) were provided by Régnier et al. (2016aRégnier et al. ( , 2018)).Several studies have been conducted on the appropriateness of the 1D theoretical prediction of site amplification (mostly in SBR) at KiK-net sites these days (e.g., Thompson et al. 2012;Pilz and Cotton 2019;Tao and Rathje 2020).These studies attribute the discrepancy between observation and theory to 2D or 3D site effects.However, we need to adjust the velocity structure at each site to make it possible to reproduce HSAF before the comparison.We have shown at more than 100 K-NET and KiK-net sites the level of velocity profile adjustment to reproduce the observed EHVRs (Kawase et al. 2018).Expecting velocity exploration by in situ PS logging to constantly provide an appropriate velocity structure for simulating the observed SAF at the site is imprudent.In addition, the site amplification for the main S-wave portion between the 1D and 2D/3D responses, except for the vicinity of the basin edge in the theoretical calculation, showed remarkable correspondence (Kawase and Aki 1989).
Figure 4 shows a good match between the 1D theoretical responses from the unified velocity structures above the seismological bedrock (Wakai et al. 2019) and the empirical HSAFs and VSAFs by GIT (Nakano et al. 2019).When we judge the matching level of prediction with observation, we need to take into account the inherent variability of multiple observations, as shown in Fig. 14 later.More optimized matching can be obtained by tuning the velocity structure based on the EHVR mentioned above or by applying empirical frequency and amplitude modification functions as proposed by Ito et al. (2021a).However, it is important to point out that basically the empirical SAF for the S-wave part from the seismological bedrock to the surface could be modeled by the entire 1D structure from the seismological bedrock to the surface.As reviewed by Ito et al. (2020), the large variations of SAF amplitude among the observations or the 1D theoretical predictions in certain studies are primarily associated with an inappropriate choice of reference.The choice of an engineering bedrock or rock formation with an S-wave velocity exceeding 760 m/s frequently used in previous studies is a good example.As shown in Fig. 6b, the borehole sensors embedded in the rock formation still exhibited significant site amplifications with strong frequency fluctuations.

Two-and three-dimensional amplifications
Past 2D or 3D site amplification studies on strong motions appear to have started after the 1968 Tokachi-Oki Japan earthquake (M7.9).As always, large, devastating earthquakes present us new information.During the 1968 Tokachi-Oki earthquake, a strong-motion seismometer of the Strong-Motion Accelerometer Committee successfully recorded accelerograms at a site in Hachinohe Harbor with a PGA exceeding 200 cm/s 2 , which was the highest value at that time in Japan.In the three-component records, they found very long-lasting (180 s) and long-period (2.5 s) ground motion.Seismologists and earthquake engineers have tried to understand the cause of these long-period ground motions, which turned out to be the basin-induced surface waves generated at the edge of the Hachinohe Plain with a thickness of approximately 300 m and an average S-wave velocity of approximately 500 m/s (Kagami et al. 1976;Midorikawa and Miura 2010).
In the Kanto Basin in Japan, long-period (6-8 s) surface waves were observed on the western side of the basin when the waves originated from shallow crustal earthquakes in the Izu Peninsula, which are considered to be basin-transduced Love waves (e.g., Kudo 1980;Yamanaka et al. 1989).Interestingly, Imamura (1931) reported Fig. 7 Cross-section of the geological structure in the east-west direction of the Ashigara Blind Prediction Test Site during the 1992 1st International Symposium of ESG (on the left) and comparisons of the predicted results in the pseudo-velocity response spectra in cm/s (average and average ± one quartile) with the concealed observed spectra reported during the first ESG symposium in red (after Midorikawa 1992 andKudo andSawada 1992) the observed seismograms during the foreshock of the 1930 Kita-Izu earthquake at the University of Tokyo, as shown in Fig. 8; however, he did not mention this prominent later phase with a period of 6 to 8 s only in the N-S component (i.e., Love-wave type basin-transduced surface waves from the Tanzawa Mountains in the west of Tokyo).
As previously mentioned, the impact of the devastating damage to high-rise buildings in Mexico City caused by long-lasting, long-period ground motions inside the soft sedimentary basin during the 1985 Michoacan, Mexico earthquake was so strong that seismologists and earthquake engineers unanimously realized the importance of quantifying ESG on strong motions for hazard mitigation and seismic risk assessment.The analysis of the doublelayered trapezoidal basin response by Kawase and Aki (1989) was one of the attempts to reproduce long-lasting, long-period ground motions as a consequence of basininduced surface waves, as shown in Fig. 9.Although time-domain responses were studied by Bard and Bouchon (1980) before the work of Kawase and Aki (1989), the simulation results in the latter study appear to be sufficient to show that the effects of multiple layers may account for the long-lasting ground motions in Mexico City.
Exactly 10 years later, the other surprisingly severe impact was caused by the 1995 Hyogo-ken Nanbu (Kobe), Japan earthquake, which created 100,000 collapsed or heavily damaged buildings and houses in a concentrated area called "the damage belt" along the downtown area of Kobe City.The cause of this damage belt became the central issue of investigations immediately after the earthquake, and Kawase (1996) provided a simple yet concrete answer to the damage belt based on the 2D velocity structure of the Osaka Basin underneath Kobe and the deconvolved input motion on the seismological bedrock from the observed records.This special amplification near the sharp edge of a basin is named "the edge effect", which represents the constructive interference between the direct S-wave from the bottom of the basin and the edge-diffracted S-waves and edge-induced surface waves.The major characteristics of the observed ground motions inside the damage belt were a relatively long period of motion (~ 1 s) and a high amplitude in peak ground velocity (> 100 cm/s).Later, we proved that such a strong velocity pulse is sufficient to cause ductile buildings to collapse or be heavily damaged, regardless of their natural frequencies (Nagato and Kawase 2004).Because of the number of near-field data during this earthquake and the importance of testing our simulation capability for reproducing the edge effect, the waveform simulations in Kobe and Osaka were the subjects of the open prediction experiment during the second ESG symposium held in Yokohama in 1998.Iwata et al. (1999) summarized the experimental settings and data provided, while Kawase and Iwata (1999) summarized the submitted results.On average, the simulated results were quite good at reproducing the observed ground motions.The most important finding of this open prediction experiment was that the variation among simulations by different parties was considerably smaller when the observed data were open, and considerably larger when the answer was not provided (i.e., at sites without data).Figure 10 compares the submitted simulations by participants at Mainshock records were all off-scale.We can see a very prominent phase of surface waves in the later part (from about 4 times the S-P time), especially in the NS component one site (RKI inside the Rokko Island) with the observations (Kawase and Iwata 1999).
As for the effects of the 2D/3D amplification, we noticed a persistent misunderstanding of the 2D/3D basin effects, to which we ought to attribute the observed discrepancy between the HSAF derived from observed ground motions and 1D theoretical HSAF based on the velocity profile at the site.Except for a site in very close proximity to the basin edge or a very narrow valley with a very high impedance contrast, the observed HSAF is the result of 1D amplification with 2D/3D additional amplification as a contribution of basin-induced surface waves (e.g., Bard and Bouchon 1980;Kawase and Aki 1989).Therefore, the fundamental peak frequency of the observed HSAF should be explained by the 1D response, and the amplitude of that peak should be close to that of the 1D response, which is only slightly smaller than the observed value.If the frequency fails to match, we should  Kawase and Iwata 1999).The top trace is the observed velocity waveform in the NS component.The vertical axis is in cm/s with a 100 cm/s shifting for each simulation consider modifying the velocity profile first.The effects of 2D/3D additional amplification are not a magic tool that bridges the gap between observation and theory.

Direct estimate of horizontal site amplification factors from earthquake horizontal-to-vertical spectral ratio
In both the GIT and GMPE analyses, the entire datasets must be analyzed simultaneously to determine the horizontal site amplification factors (HSAFs) and vertical site amplification factors (VSAFs*) with the common reference.It is also difficult to apply the GIT or GMPE to events so small that their magnitude or source locations cannot be determined accurately or their seismograms cannot be recorded simultaneously at multiple sites.To fully exploit the ease of a single-station observation, Ito et al. (2020) recently proposed a vertical amplification correction function (VACF) to directly estimate HSAF from the horizontal-to-vertical spectral ratio of earthquakes (EHVR).This correction is possible because VSAF* is a relatively gentle function of frequency common to almost all sites.As mentioned previously, VSAF* includes amplitude correction at the seismological bedrock between the vertical and horizontal components, as Ito et al. (2020) used.Therefore, we should use VSAF* for the VACF calculation because the common site amplification factor relative to the same horizontal seismological bedrock motion (H B ) is needed to obtain HSAF from the EHVR.However, the VSAF in Fig. 4 is corrected with the bedrock EVHR (which means that V B H B R was extracted) to compare the 1D P-wave theoretical amplification directly.
In the GIT, the S-wave part Fourier spectrum of horizontal motion, F S_ij , of earthquake i observed at site j is decomposed into the logarithmic sum of the source term S S_i , the path term P S_ij , and the site amplification factor H S_j (which corresponds to sHSAF as the S-wave part HSAF), as shown in the following equation: Likewise, the S-wave part Fourier spectra of the vertical motion, G S_ij , is broken down as follows: This equation assumes that ground motion is propagated as S-waves until it reaches the seismological bedrock immediately below the observation site.Then, a part of these S-waves is converted to P-waves by the scattering within the propagating medium.These P-waves propagate to the ground surface and emerge as vertical (1) log F S_ij = log S S_i + log P S_ij + log H S_ j . (2) motions on the surface.V S_j is the vertical site amplification factor for the P-wave.The third term in the righthand side of Eq. ( 2), V B H B R, is a coefficient for converting S-wave amplitude into P-wave amplitude on seismological bedrock.It is the inverse of EHVR on the seismological bedrock in the diffuse field concept (DFC, Kawase et al. 2011;Ito et al. 2020).
This coefficient arises because we want to use the same reference condition for both horizontal and vertical components.Our primary focus is the amplitudes of S-wave that are propagated and scattered through the medium from the hypocenter to a point on the seismological bedrock immediately below the site.Using such a value as a common reference, we can obtain the S-wave amplification factor on the surface relative to the amplitude of seismological bedrock S-waves.
On the other hand, P-wave amplitudes on the same seismological bedrock immediately below the site should be defined as the values relative to the same S-wave reference on the seismological bedrock.Once we convert the S-wave amplitude to the P-wave amplitude on the seismological bedrock, we can obtain the P-wave amplification factor on the surface relative to the amplitude of seismological bedrock P-waves, V S_j .Notably, H S_j and V S_j *V B H B R in these equations are the specific terms derived from the GIT, whereas sHSAF and sVSAF* with a prefix of s defined here are the general terms referring to site amplifications in the horizontal and vertical directions for the S-wave part.The prefix "s" from now on is introduced because we will discuss later the site amplification factors for the whole duration of motion.The separation of V S_j and V B H B R is only conceptual, and in the actual GIT analysis, sVSAF* is directly obtained as (V Then the relationship between the horizontal-tovertical spectral ratios of the S-wave part, sEHVR, and sHSAF are discussed.From Eqs. (1) and ( 2), the sEHVR is expressed as follows: sHSAF is expressed as follows: where < > denotes the log-averaging operation.Equation (4) indicates that sHSAF can be obtained by multiplying sEHVR by sVSAF* as the vertical amplification relative to the horizontal bedrock motion.Ito et al. (2020) determined an empirical correction function averaged over multiple sites, which is called the vertical amplification correction function (sVACF) for the S-wave part, as Kawase et al. (2019) for the MHVR.After calculating averaged sVACF as seen in Eq. ( 4), sVACF yields sHSAF at an arbitrary site using the following simple equation: where sVACF is the log-averaged spectral ratio of the vertical amplitude on the ground surface with respect to the horizontal amplitude on the outcrop of the seismological bedrock for the S-wave part of the earthquakes.The sVACF can be directly determined using the GIT for vertical motion with respect to horizontal motion on the reference bedrock.In Fig. 11, we show the sVACF for all sites without a category (Ito et al. 2020) on the left and examples of the simulated sHSAF (broken red lines) compared with the observed sHSAF (black lines) on the right.The maximum amount of correction is 2.0 at frequencies higher than 1 Hz.In Ito et al. (2020), sVACF may only be valid from 0.4 to 15 Hz.Quite recently, we proposed an improved version of sVACF as a function of Vs30 (Wang et al. 2023).

Inversion method based on the earthquake horizontal-to-vertical spectral ratio
For the 1D structure down to the seismological bedrock, we can use the EHVR inversion method (5) sHSAF = sEHVR • sVACF, implemented by Nagashima et al. (2014) based on the diffuse field concept (DFC) for the sum of the different events proposed by Kawase et al. (2011).If a borehole sensor is available, then a joint inversion of the EHVR and SBR can be performed.As shown in Nagashima and Kawase (2019), we can also determine the damping (Q values) of the layers above the borehole sensor depth in the joint inversion.The EHVR inversion method is significantly robust and more than 100 sites in Japan have verified its feasibility and efficiency (e.g., Ducellier et al. 2013;Fukihara et al. 2015;Kawase et al. 2018).We have extended the applied regions to five sites in France (Ito et al. 2021b) and 23 sites in Switzerland (Chieppa et al. 2018(Chieppa et al. , 2023)).Also, we have a proposal for a new technique from Iran based on the same DFC (Ashayeri et al. 2023).
The success of the EHVR inversion is contingent on two essential factors.One is the need to model the velocity structure down to the seismological bedrock.According to the theory of the sEHVR derived by Kawase et al. (2011), sEHVR on the surface is the ratio of the transfer function of the vertically incident S-wave for all layers above the seismological bedrock to the transfer function of the vertically incident P-wave for the same layers, with the amplitude correction factor at the seismological bedrock as follows: where α H and β H are the P-and S-wave velocities of the seismological bedrock, respectively.The left-hand side simply corresponds to sEHVR in this letter.G denotes Green's function of plain body waves on the surface, and TF is the transfer function from the seismological bedrock to the surface.Zeros within the parentheses indicate that the observed point is on the surface.The impact of modeling from the seismological bedrock should not be underestimated because we used both frequency fluctuations (peaks and troughs) and the absolute sEHVR amplitudes.We simulated the effects of the assumed layers, as shown by Kawase et al. (2018) in Fig. 12.In the sEHVR, all the resonances only within a shallow layer and those of the whole basin structure are included to represent the effects of the entire structure accordingly.If we used a bedrock layer with an S-wave velocity lower than the seismological bedrock with an S-wave velocity higher than 3 km/s in the sEHVR inversion, the resultant structure would not be an actual structure.
Another important point for successful inversion is P-wave velocity modeling.The P-wave transfer function exhibited smooth characteristics, as shown in Fig. 11a, due to the relatively small P-wave velocity contrast.However, the P-wave transfer function can somewhat affect the precision of the delineated S-wave velocity structure.Therefore, we require a suitable mapping function ( 6) from the S-wave velocity to the P-wave velocity to perform a stable sEHVR inversion.To that end, Nagashima and Kawase (2021) obtained new empirical relationships between the S-and P-wave velocities from the published data on K-NET and KiK-net, which include the water saturation in the target layers.

One-dimensional theoretical evaluation for empirical horizontal site amplification factors by generalized inversion technique
As shown in Fig. 4 et al. 2015, 2019).However, a comparison of the gross image of the matching for the entire region in terms of the residual between the 1D theoretical sHSAF by the velocity profiles and the empirical sHSAF shows the need for modifications in terms of both frequency and amplitude.This need for modification is attributed to the inappropriateness of the prescribed velocity profile at the target site, and the two-/three-dimensional (2D/3D) effects that emerged even within a relatively short (5-15 s) S-wave window, including but not limited to the near-surface topography.Therefore, Ito et al. (2021a) proposed a theoretical amplitude correction method in terms of FMR and AMR, that is, the frequency and amplitude modification ratios onto the theoretical sHSAF.Figure 13 shows spatially interpolated maps for the FMR and AMR in the Kanto and Tokai regions.The need for a relatively complicated  section all over Japan is an urgent task that should be performed in the near future.

Basin-induced surface waves in the whole-duration horizontal site amplification factors
Our GIT studies use only observed spectra for a relatively short duration (5-15 s) to extract sHSAF and sVSAF* for the S-wave part.After separating the source and path spectra of the S-wave, we obtained HSAF and VSAF* for the whole observed duration as the average residuals with respect to the S-wave bedrock motion spectra for each source and path (Nakano et al. 2019;Nakano andKawase 2019, 2021a;Nakano 2021).We call them wholeduration HSAF and VSAF* (wHSAF and wVSAF*), which were found to have significantly larger amplitudes in the longer period range than the corresponding sHSAF and sVSAF*.This phenomenon can be attributed to the significant contribution of basin-induced surface waves in wHSAF and wVSAF*.Considering the original Eqs.( 1) and ( 2) but for the whole duration spectra, the whole duration spectra F W_ij of earthquake i observed at site j consists of the same source term S S_i and path term P S_ij as the S-wave part, but has a different site amplification factor H W_j associated with the whole duration, that is wHSAF, as shown in the following equation: Similarly, the whole duration spectra for vertical motion, G W_ij , are broken down as follows: where < V w_j × V B H B R > corresponds to wVSAF*.These formulas show that we can obtain wHSAF and wVSAF* as residuals from the input S-wave spectra at the seismological bedrock.
Equations ( 7) and ( 8) yield wHSAF and wVSAF* for the whole duration of motion, respectively, whereas Eqs. ( 1) and ( 2) do sHSAF and sVSAF* for the S-wave part at hand, respectively.From these two pairs of SAFs, the whole duration to the S-wave part ratio, the WSR, is determined (Nakano 2021;Nakano and Kawase 2023).Because we separated the SAFs for the horizontal and vertical components, we have WSRh for the horizontal component and WSRv for the vertical component.At a site without observed data, sHSAF can be calculated using a 1D S-wave theoretical ground response analysis based on the velocity structure proposed, such as a unified model by Senna et al. (2014Senna et al. ( , 2018Senna et al. ( , 2019)), and Wakai et al. (2019), together with the FMR and AMR corrections (Ito et al. 2021a).Then wHSAF could be calculated as: where sHSAF_theory is the theoretical sHSAF.Given the sEHVR for the S-wave part, sHSAF can be calculated also by using the VACF method reported by Ito et al. (2020) for the S-wave part (Kawase et al. 2021).
Considering sHSAF and wHSAF from the GIT in Eqs. ( 1) and ( 7), Fig. 14 shows examples of sHSAF and wHSAF for three K-NET sites in Akita Prefecture.Evidently, the difference in the frequency range exceeding 2 Hz is insignificant; however, the additional amplification at frequencies lower than 1 Hz is substantial.Furthermore, additional amplification was strongly site-dependent.This is natural because long-lasting basin-induced surface waves occur inside soft and deep sedimentary basins.In contrast, there are no such surface wave contributions at sites on stiff soil or rock formations.
Figure 15 shows the entire Japanese map of WSRh based on the GIT results of Nakano (2021) and Nakano and Kawase (2023) at two frequencies, 0.5 Hz and 1.0 Hz, the values of which were interpolated by the surface function of GMT (Smith and Wessel 1990).Evidently, the WSRh is larger inside the sedimentary basins and volcanic calderas and close to 1 in mountainous areas.Similar to the case for sHSAF in Fig. 5, the southwestern side of Japan shows a systematically smaller WSRh than the northeastern side of Japan.Moreover, WSRv can be calculated from sVSAF* in Eq. ( 2) and wVSAF* in Eq. ( 8), the entire map of which is shown in Fig. 15 in the same way as WSRh.WSRv has a higher amplitude than WSRh inside the sedimentary basins and volcanic calderas.Figure 16 shows examples of the spectral shapes of WSRh and WSRv at three representative sites.Again, larger amplitudes of the WSRv are observed at sites inside the sedimentary basins.Figure 17 shows the overall image of the difference in WSRh and WSRv as the logarithmic average of WSRh and WSRv for all 2593 sites, along with the one-to-one correspondence of WSRh and WSRv at four representative sites in Akita Prefecture.Remarkably, WSRv was approximately 25% larger than WSRh, on average.As observed in the standard deviation, the variability was almost the same between them.The difference between WSRh and WSRv at one site is significantly smaller than their variability from site to site; therefore, they are virtually identical on average.The extension of the VACF method proposed by Ito et al. (2020) for the wHSAF estimation was presented by Kawase et al. (2021) at the 6th ESG symposium, in which the following two equations are proposed: (9) wHSAF = sHSAF theory • WSRh, (10) where wEHVR in Eq. ( 11) is the horizontal-to-vertical spectral ratio for the whole duration.As observational evidence, EHVR of the coda part is very similar to sEHVR (Kawase et al. 2018), and Eqs. ( 10) and ( 11 11) is easier to use for a new site because the extraction of the S-wave part from the onset of the S-wave is not required as in the case of Eq. ( 10).
The WSR approach is advantageous because all the complex site amplification factors are associated  Nakano (2021).Red dotted curves are the approximate location of the Itoigawa-Shizuoka tectonic line.Interpolation is done by GMT (Smith and Wessel 1990) Fig. 16 Examples of spectral shapes of WSRh and WSRv at one site outside of a basin and two sites inside of a large basin (OSK006 is close to the shoreline of Osaka Bay and TKYH02 is on the west of Tokyo with a borehole sensor at the depth of 2753 m)."Original" means those from the observed wHSAF and sHSAF, whereas "Estimation" does those interpolated spatially by the two-step interpolation scheme described in Nakano and Kawase (2021b).As seen in Fig. 15, the amplitudes of the low-frequency WSRv tend to be larger than those of WSRh at the sites inside a large basin Fig. 17 Overall logarithmic averages of WSRh and WSRv and its deviations for all the sites used (left) and one-to-one relationships of WSRh and WSRv for the same frequency at four selected sites in Akita Prefecture (right) with sHSAF for the S-wave part, and the 2D/3D basin response contribution in wHSAF is considered to be a smooth function of multiplication.We called the ratio between wHSAF for the whole duration and sHSAF for the S-wave part as WSR.WSR would be a smooth function because the peak frequencies of the S-wave 1D resonance would correspond to the peak frequencies of the surface waves as long as the impedance contrast between the sediments and bedrock is sufficiently high.The phenomenon can be seen in Fig. 9 for the theoretical response of a simple trapezoidal basin.The success of the sHSAF, wHSAF, and WSR is largely attributed to the representation of these values as amplification with respect to the seismological bedrock with an S-wave velocity exceeding 3 km/s.
Over the long history of ESG studies, researchers have been reluctant to use the HSAF from the seismological bedrock or site proxies related to the seismological bedrock, such as Z3.2 (the depth to the layer with an S-wave velocity of 3.2 km/s) or Tz3.2 (the travel time to Z3.2) to represent the HSAF from the seismological bedrock, because the velocity structure down to the seismological bedrock could not be easily delineated because of the prohibitive cost and effort.Currently, the sEHVR inversion method described above (e.g., Nagashima et al. 2014) can be used, and it facilitates the delineation of the S-wave (and P-wave) velocity structures down to the seismological bedrock.The resultant velocity structure by the sEHVR inversion yields the observed HSAF through a simple 1D S-wave ground response analysis because sEHVR contains all the essential information required to reproduce sHSAF, as shown by the VACF correction method.

Use of microtremor horizontal-to-vertical spectral ratio in prediction of the horizontal site amplification factors for S-wave
The use of microtremors as a substitute for information not based on the target prediction information (e.g., site-specific spectra or amplification of earthquakes) should be addressed.The prediction schemes described in Sects."Direct estimate of horizontal site amplification factors from earthquake horizontal-to-vertical spectral ratio" and "Inversion method based on the earthquake horizontal-to-vertical spectral ratio" are directly based on the sEHVR at the target site; therefore, the use of these schemes is considered to be a subject of circular argument, as mentioned in Sect."Site classification or site proxy approaches".Conversely, the use of MHVR is different from the direct use of the sEHVR for sHSAF prediction of earthquakes because the spectral characteristics of the EHVR and MHVR should not be the same (e.g., Satoh et al. 2001;Horike et al. 2001) and MHVR is significantly easier to obtain at the target site than earthquakes.
There have been numerous studies on the use of the MHVR in the past after Nogoshi and Igarashi (1971) and Nakamura (1989), as reviewed by Molnar et al. (2022) as a series of studies on non-invasive methods for velocity exploration (Yong et al. 2022).After a long discussion on the optimal use of the MHVR and its theoretical representation, our community finally obtained a complete theoretical solution for MHVR (Sánchez-Sesma et al. 2011) and a subsequent inversion code to obtain the velocity structure only from a single-station measurement of microtremors (García-Jerez et al. 2016).However, our proposed approach differs from these physical solutions of the MHVR.It is based on the empirical ratio between the EHVR and MHVR at the same site, namely, EMR (Kawase et al. 2018).These authors proposed five different EMRs (Fig. 18) for sites with different f0 values from MHVRs, namely for sites with 0.2 ≤ f0 < 1 Hz, 1 ≤ f0 < 2 Hz, 2 ≤ f0 < 5 Hz, 5 ≤ f0 < 10 Hz, and 10 ≤ f0 < 20 Hz based on the observed earthquake and microtremor data at 100 K-NET and KiK-net sites in Japan.They found that 86 out of 100 sites yielded an EMR of 1.0 at the fundamental MHVR frequency f0.This observation indicates the same peak frequency and amplitude at and around the fundamental peak in both the EHVR and MHVR.In addition, 14 sites were not used for EMR determination because several sites had f0 values outside of the frequency range of interest (0.2 to 20 Hz), whereas the other sites had MHVRs that did not match the EHVRs around f0.Although a clear answer is not available about why there are sites without good matching at present, several examples have been provided at Mashiki town, Kumamoto, Japan, where we observed microtremors two weeks after the 2016 Kumamoto earthquake sequence and found that there were Fig. 18 Japanese EMRs for all five categories based on the fundamental peak frequency (f0) of the MHVR (after Kawase et al. 2018).The horizontal axis is the normalized frequency with respect to f0 so that we will expand or shrink the horizontal axis according to f0 for the target site and interpolate EMR before applying EMR to MHVR several sites without similar MHVRs to those of nearby sites (Sun et al. 2020).We observed microtremors again two years after the event and found a drastic improvement in the MHVR similarity to those at adjacent sites.This strongly suggests that when both peak frequency and amplitude around the fundamental frequency fail to show suitable correspondence between the MHVR and EHVR or MHVRs around the target site, MHVR should be measured again at a different date and time, because such peculiar characteristics of the MHVR may not be associated with homogeneous local noise sources around the target site.
The use of the EMR is simple but effective.We first obtain f0 from the MHVR and then multiply MHVR with the corresponding EMR, which is normalized to f0, thus necessitating an interpolation operation after the multiplication of f0 with the normalized frequency of the EMR.The resultant HVR is called pseudo-EHVR or pEHVR.The tables of EMRs with five categories, the detailed procedure for applying the EMR method, and the validation exercises are reported by Kawase et al. (2018).The effectiveness of using the pEHVR in the inversion method of the EHVR proposed by Nagashima et al. (2014) was presented (Nagashima et al. 2023) in the blind prediction experiment Step 1 during the ESG6 (Matsushima et al. 2021(Matsushima et al. , 2022;;Chimoto et al. 2021Chimoto et al. , 2022)).Comparison of the predicted (pEHVR) and observed EHVR from the aftershock records opened after the blind prediction and the resultant velocity structures based on the inversions from pEHVR and EHVR are shown in Fig. 19.These results are consistent with each other.
Field and Jacob (1995) compared four methods to delineate sHSAF at four sites in Oakland, California: SSR, GIT, EHVR, and MHVR.They reported (Fig. 4 of their study) that the frequency amplification factors from SSR, GIT, and EHVR were quite similar to each other from 0.4 to 10 Hz.When we look at the amplification factor from the MHVR, the fundamental peak frequency (~ 0.7 Hz) and amplitude (~ 4) are quite similar to those of the three methods based on earthquake ground motions.In contrast, high-frequency amplitudes and fluctuations above 1 Hz are largely deficient in the MHVR.The gap between the EHVR and MHVR reported by Field and Jacob (1995) is consistent with the EMR characteristics (Category 1) shown in Fig. 18.
The optimal use of the EMR method is the direct estimation of sHSAF from the converted pEHVR and the VACF method that can transform EHVR into HSAF, as shown in Sect."Direct estimate of horizontal site amplification factors from earthquake horizontal-to-vertical spectral ratio" (Kawase et al. 2019).As pEHVR can be a viable substitute for the observed EHVR, the double correction of both EMR and VACF can effectively convert MHVR into sHSAF.Nakamura (2019) insisted that the so-called "Nakamura method" is not a method for calculating the MHVR, but should include the interpretation of the equivalence of the MHVR and sHSAF.In other words, he insisted that MHVR could be a direct substitute for sHSAF, a concept referred to as "the Nakamura assumption" (Ito et al. 2020).Our series of analyses showed that this was not the case for most sites.As mentioned by Ito et al. (2020), the observational and theoretical studies that reported positive validation evidence for the Nakamura assumption in the past, including those by his own group, would have been positive because of their primary use of the engineering bedrock as a reference, rather than the seismological bedrock, which we used for the double correction method in Kawase et al. (2019).The insufficient amplification from the engineering bedrock to the surface might compensate for the insufficient amplitude of the MHVR compared to the real sHSAF from the seismological bedrock.
Finally, the absolute value of EMR should depend on the average velocity structure of the target basin because the EMR is derived from the difference in the EHVR as the ratio of horizontal and vertical transfer functions (Kawase et al. 2011) and the MHVR as the ratio of the imaginary parts of the horizontal and vertical Green's functions (Sánchez-Sesma et al. 2011), making it sitedependent.In the case of a basin with a small number of observation sites with EHVR and MHVR, Ito et al. (2021b) proposed a method for modifying the Japanese EMR to fit the EMR specific to the target basin in Fig. 19 Based on the distributed data of microtremors for the blind prediction experiment during the ESG6 symposium in 2021, we compare the observed MHVR (green), the converted pEHVR from MHVR (red), and the observed EHVR (blue) from aftershock records disclosed after the blind prediction on the left.On the right, we compare the resultant S-wave velocity profiles inverted for pEHVR (red) and EHVR (blue), shown in Nagashima et al. (2023).The profile with an orange line is the preferred model distributed by the committee of the blind prediction (Matsushima et al. 2021(Matsushima et al. , 2022) ) constructed by the boring data down to 38 m and the deep structure proposed by Senna et al. (2018).The inclined dashed and dotted lines on the velocity models show resonant frequencies of the corresponding layer interfaces based on the one-quarter wavelength law Grenoble, France.They showed that such a basin-specific EMR could work for other sites inside the target basin with only MHVRs.Before applying the EMRs shown in Fig. 18 to a targeted basin that is significantly different from typical basins in Japan, we need to have a minimum number of sites (~ 5) with earthquake and microtremor records to confirm their applicability.

Use of response spectra versus Fourier spectra
As shown in the GMPE list by Douglas (2022), most of the previously proposed GMPEs used response spectra, either pseudo-velocity response spectra or pseudo-acceleration spectra with 5% damping.The response spectra are used mainly because they are convenient for engineering purposes, such as probabilistic seismic hazard analysis, if one can predict the response spectra directly.From a theoretical viewpoint, however, response spectra are not ideal as a target for the quantitative evaluation of strong motion characteristics, either for the source term, path term, or site term, because the value at one frequency is a complex function of the values at other adjacent frequencies.Bora et al. (2016) attribute this to the transfer function with a peak at the resonant frequency of a single-degree-of-freedom system with 5% damping; however, this system has a large side lobe on both sides of the resonant frequency.As a result of the interrelationship between adjacent frequencies, the prediction variability of the response spectra could be difficult to reduce, regardless of the use of complex formulae in the GMPE (e.g., Strasser et al. 2009; Fig. 3).Please note that the response spectra are the maximum response of a singledegree-of-freedom system with a given damping rested on the ground, not the response of the ground itself.
However, the problem associated with the strong smoothing effect in the response spectra due to the influence of the adjacent frequency energy would not be a significant problem, as long as our target is the sHSAF with a short duration.Evaluating the site amplification of the whole duration of motion, wHSAF, in response spectra would be really challenging.As shown in Figs. 14,15,16,17, we observed a strong effect of the observed long duration of motion in the lower frequency range at sites inside sedimentary basins such as the Kanto and Osaka Basins.The primary cause of such a long duration of motion is basin-induced surface waves converted from incident S-waves at the edge of the basin (Kawase and Aki 1989;Kawase and Sato 1992;Kawase 1993).We plot the sHSAF and wHSAF at selected sites in Fig. 20 for both the Fourier spectra and pseudo-velocity response spectra.A significant reduction of the wHSAF amplitudes is observed in the response spectra compared to those of Fourier spectra, because of the fundamental nature of the response spectra, which is the maximum response of a single-degree-of-freedom system with 5% damping and therefore, not very sensitive to the whole duration of motion; as long as it includes 3 to 4 packets of wavelets at the resonance frequency we can obtain a saturation state (e.g., Kramer 1996).In the regression analysis to construct GMPEs, one would usually use response spectra for the whole duration.However, the resultant SAFs do not comprehensively capture the effects of site response coming from the whole duration, as is clearly shown in Fig. 20.This deficiency is the most important reason for using Fourier spectra for ESG studies, instead of response spectra.In addition to the deficiency in the wHSAF, the response spectra are incapable of capturing the decreasing amplitudes of Fourier Spectra in both lower and higher frequency ranges in sHSAF.This is because the amplitude of sHSAF at a specific frequency in response spectra is not determined solely by the contribution from that frequency, but is determined as a result of convolution with contributions from other adjacent frequencies.

Reference issue: rock outcrop or seismological bedrock
As shown in the GMPE list reported by Douglas (2022), most of the previously proposed GMPEs used rock outcrops or stiff soil sites as a reference for their GMPEs.Numerous studies use 760 m/s or higher, but usually less than 1.0 km/s in terms of Vs30, to select the reference sites for their GMPEs.The rock outcrop sites with such a low Vs as a reference are used primarily because it is convenient as a real observation site or site.It is almost impossible to find a rock-outcrop site with a surface Vs higher than 3.0 km/s, which is the condition that we can call the seismological bedrock, where we can assume no site amplifications.To overcome such difficulty, Kawase and Matsuo (2004) and Nakano et al. (2015) used the deconvolution technique to calculate the outcrop motions at one rock site of KiK-net, YMGH01, with the estimated Vs of 3.45 km/s as a reference of their GIT analyses.Owing to the use of the corrected outcrop spectra on the seismological bedrock, the resultant HSAFs at the other sites can be considered as absolute HSAFs relative to the input S-wave ground motions observed on the surface of a half-space with seismic sources.
To delineate the level of the side effects on the evaluation of sHSAF using sites with Vs30 higher than 760 m/s, we plotted the average sHSAF for 76 sites in K-NET and KiK-net with 760 < Vs30 < 1000 m/s in Fig. 21.The Vs30s of the K-NET sites are extrapolated values using an empirical formula (Boore et al. 2011).The geometrical mean of sHSAF at these 76 sites shows non-negligible amplitudes higher than 2 in the frequency range from 1 to 10 Hz.Several sites showed significantly higher amplifications; therefore, careful selection is necessary when we use a group of sites as a reference in the GMPE construction.However, this level of sHSAF at the reference site (or sites) is carried over to the estimated source spectra if we use a site or sites with certain amplifications, thus providing higher corner frequencies and estimated stress drops.This is because the constraint condition of the fundamental equations in Eqs. ( 1) and ( 2) or ( 7) and ( 8) should determine the relative distribution ratio between the source and site terms, even if the path term is fixed.

Direct use of spectra, not extracted values like f0
Historically, most of the GMPE developments have used a site condition proxy or proxies to consider the effects of ESG on the observed strong motions such as Vs30 and the depth to the prescribed Vs and Z Vs .Derras et al. (2017) and Zhu et al. (2020) reviewed the effectiveness of various recent proposals for additional site condition proxies rather than a single proxy of Vs30.Numerous reports recommended the fundamental period T0 or frequency f0 (reciprocal of T0) for the optimal proxy in addition to Vs30.This makes sense, as reviewed in Sect."Site classification or site proxy approaches" because T0 or f0 Fig. 20 Comparison of sHSAF and wHSAF for Fourier spectra derived by the GIT analysis (top three panels) and sHSAF and wHSAF for pseudo response spectra with 5% damping derived in the same way as Fourier spectra.We should note that the amplitudes of sHSAF for response spectra are close to that for Fourier spectra, with stronger smoothing effects, however, the amplitudes of wHSAF for response spectra are strongly underestimated, especially in the lower frequency range should reflect the most important frequency characteristics to predict the sHSAF, especially for amplitudes in the frequency range below 1 or 2 Hz.However, when the MHVR, EHVR, or SBR are used to detect site-specific f0 to better predict the sHSAF, the comprehensive use of the frequency values of one of these observed spectra becomes very effective for the quantitative prediction of the entire frequency-dependent sHSAF, as evidenced in the recent studies by Zhu et al. (2022) and Zhu et al. (2023).Similarly, Z1.0 or Z800 is a good additional proxy for evaluating the sHSAF in the lower frequency range; however, if we have a whole velocity profile down to the layer with Vs = 800 m/s or 1 km/s, using only the depth to the layer with a prescribed Vs is not a logical choice for obtaining the entire frequency-dependent sHSAF.The use of a vector value of the entire velocity profile or the entire theoretical sHSAF estimate based on the 1D ground response analysis from the velocity profile as site condition proxies would be more optimized than the sole integration of Z800 or Z1.0 as an additional proxy.Therefore, the inverted velocity structure from the seismological bedrock based on the EMR method proposed by Kawase et al. (2018), the direct use of the empirical VACF method proposed by Kawase et al. (2019) for the MHVR, or Ito et al. (2020) for the EHVR should be more optimized choices for predicting the sHSAF for the entire frequency range of interest (e.g., 0.12 to 20 Hz).

Nonlinearity
Nonlinearity in sHSAF is an important factor for scenario-based source-and site-specific prediction of ground motions in the near-field region of a sufficiently large magnitude earthquake.Traditionally, the nonlinearity of the soil response in a shallow part of the ground has been modeled by the nonlinear (including equivalent linear) ground response analysis, where the level of the shear strain of each layer is used as the soil nonlinearity parameter in the evaluation of the shear rigidity degradation due to soil nonlinearity (e.g., Hardin and Drnevich 1972;Schnabel et al. 1972;Joyner and Chen 1975;Mohammadioun and Pecker 1984).In the early stages of soil nonlinearity evaluation studies after the 1987 Loma Prieta and 1994 Northridge, California earthquakes, as well as the 1995 Kobe, Japan earthquake, weak and strong sHSAF or SBR comparisons were performed to delineate the levels of nonlinearity at observation sites with high PGA or PGV (e.g., Chin and Aki 1991;Satoh et al 1995b;Aguirre and Irikura 1997;Field et al. 1997;Trifunac and Todorovska 1998;Beresnev and Wen 1996;Hartzell 1998).We refer to the reviews of Régnier et al. (2016aRégnier et al. ( , 2018) ) for more recent studies on empirical and theoretical approaches to soil nonlinearity.
There are also a number of studies on GMPEs that consider nonlinear soil behavior in the site amplification coefficient (e.g., Walling et al. 2008;Sandıkkaya et al. 2013;Seyhan and Stewart 2014;Kamai et al. 2014;Loviknes et al. 2021).However, approaches that directly include the nonlinear site effects in GMPEs have "the delusion of nonlinearity" problem as averaged characteristics because the emergence of nonlinearity is quite limited in terms of the number of observed data, while its frequency characteristics are strongly dependent on the site-specific velocity profile, in addition to the PGA, PGV, or effective shear strain levels.This necessitates the extraction of the emergence of common nonlinear characteristics using the linear site-specific sHSAF as a reference.From this point of view, the approach proposed by Régnier et al. (2013) and Régnier et al. (2016b), and later extended by Derras et al. (2020), where the relative Fourier amplitude ratios between the nonlinear HSAF with respect to the linear HSAF are used, makes practical sense for the representation of soil nonlinearity.
Quantitative extraction of soil nonlinearity from a high-amplitude record of a single event at a specific site is significantly more difficult than that from many smaller events because of the time-varying characteristics of the observed records from a specific event at a specific site.To overcome this challenge, even in the era without theoretical formula for the EHVR (Kawase et al. 2011), researchers have used the EHVR to delineate the soil Fig. 21 sHSAFs of all the 76 sites with the Vs30 range of 760 < Vs30 < 1000 m/s (gray lines) and their geometrical mean and mean ± one standard deviation (red lines).The average relative amplifications of these rock sites at 1 Hz and 10 Hz are 2.2 and 3.2, respectively nonlinearity between weak and strong observed motions (Dimitriu et al. 1999;Dimitriu 2002;Wen et al. 2006;Noguchi and Sasatani 2011;Régnier et al. 2013;Ren et al. 2017).As discussed in Sect."Direct estimate of horizontal site amplification factors from earthquake horizontal-tovertical spectral ratio", we can use the EHVR for sourceand site-specific strong motions with large amplitudes based on the following two assumptions: (1) no strong nonlinearity is observed in the vertical component; and (2) the same DFC assumption in Eq. ( 6) for EHVR can be used for a single-event record.The nonlinear sHSAF can then be delineated from the sEHVR for strong motions based on either the empirical VSAF* from weak motions or theoretical VSAF* from the Vp profile inverted by the EHVR.After obtaining the sHSAF, we calculate (deconvolve) input motions on the seismological bedrock (e.g., Sun et al. 2021;Nagashima and Kawase 2022;Nagashima et al. 2022;Fukutake et al. 2023) in the same manner as proposed for linear sHSAFs.Such backwardly calculated input motions can be used to reproduce surface ground motions at arbitrary sites based on inverted velocity structures from sEHVR.
The approach adopted by Wang et al. (2021Wang et al. ( , 2023) is different from these site-specific methods but is applicable to any synthetic linear ground motion, although the same two assumptions are used.First, Wang et al. (2021Wang et al. ( , 2023) ) determined the frequency modulation factor α and amplitude modulation factor β by comparing weak and strong motion sEHVRs.The combination of α and β was named as a nonlinear correction factor.Figure 22 shows an example of the determination process of α and β by a simple grid search to match the modified sEHVR from weak motions with the observed nonlinear sEHVR throughout the frequency range (0.5 to 20 Hz) on the left and the resultant nonlinear sEHVR using the best nonlinear correction factors for the specific spectra on the right.After determining the nonlinear correction factors from the averaged sEHVR of weak motions at all the sites with a strong motion (where PGAs are higher than 100 cm/s 2 ), they plotted these factors as a function of the observed PGAs of strong motions and obtained regression curves for α and β, as shown in Fig. 23.They observed a very systematic increase in α (positive α means decrease in peak frequency) as PGAs increased and a rather scattered but positive tendency in β (positive β means increase in peak amplitude).This approach is advantageous because it can predict nonlinear sEHVR (as a function of PGA) from the weak motion sEHVR upon the determination of the empirical relationships of nonlinear correction factors in Fig. 23.Then sEHVR can be converted to sHSAF using the VACF method mentioned in Sect."Direct estimate of horizontal site amplification factors from earthquake horizontal-to-vertical spectral ratio".Because α and β are functions of the surface PGA including nonlinearity, iterations of deconvolution and convolution need to be performed if we first obtain a simulated linear ground motion by the statistical Green's function method from the source, path, and site terms separated by GIT (e.g., Nakano and Kawase 2021b).

Borehole horizontal site amplification factors for the S-wave
As mentioned in Sect."Site classification or site proxy approaches", researchers have been using the SBR (surface-to-borehole ratio) to delineate the velocity profile and damping estimate whenever there is a surfaceborehole combination (vertical array).As long as the main portion of the incoming wave is an S-wave, we can effectively invert these parameters between the two sensors.However, numerous studies have been published based on the assumption that borehole sensors can be good candidates for reference if they are embedded inside hard rocks with sufficiently large Vs.As has been reported in the textbook (e.g., Kramer 1996) or in recent publications (e.g., Thompson et al. 2012;Tao and Rathje 2020), the sHSAF inside the borehole should not be treated as a good candidate for the reference, even if the embedded layer of a borehole sensor has sufficiently large Vs, such as Vs > 3.0 km/s.
There are two problems with borehole sHSAF: the cancellation effect of the reflected waves at the ground surface and arriving at the borehole sensor.At the resonance frequency, the upgoing wave is reflected at the surface and returns as the downgoing wave with the opposite phase, so that periodic interferences occur at the borehole depth.In the absence of intrinsic or scattering damping, the amplitudes at these cancellation (pseudo-resonance) frequencies are zero.As shown later, the actual cancellation effect was not significant owing to the scattering effects.This cancellation effect should exist regardless of the S-wave velocity, except for the case with a very deep borehole sufficient to absorb the downgoing wave energy or a very strong input sufficient to liquefy shallow soil layers.The second is the emergence of a common resonance effect between the interface below the borehole level and the surface.The resonance frequency owing to such deep interfaces should be lower than the fundamental peak frequency associated with the cancellation effect.
In the past, we did not attempt to separate sHSAFs from borehole sensors of KiK-net based on GIT (e.g., Nakano et al. 2015Nakano et al. , 2019;;Nakano and Kawase 2021a) because this would double the weight of KiK-net sites on the source and path terms.The sHSAF of the borehole sensor can be obtained from the sHSAF of the surface sensor and SBR.However, we decided to perform the GIT analysis for both surface and borehole sensors of KiK-net (Kawase et al. 2022), because there are an increasing number of researchers who underestimate these two effects mentioned above and attempt to use borehole spectra as reference.
Figure 24 shows an example of the surface and borehole sHSAFs at the AKTH19 KiK-net station, together with the SBR.As can be seen, sHSAF on the surface (sHSAF_ surface, the same as sHSAF) shares common characteristics with SBR above the fundamental peak frequency of the SBR.However, the amplitude of sHSAF_surface was slightly larger than that of the SBR.The sHSAF of the borehole (sHSAF_borehole) shows troughs at the peak frequencies of the sHSAF_surface.Although the amplitudes of these troughs are not close to zero (~ 1.0), the ratio of the amplitude of the surrounding frequencies on both sides of the trough is approximately 1/2.This means that the peak amplitude of the SBR is formed as a result of both the real soil amplification and the cancellation effect at the borehole.The slightly larger sHSAF_surface amplitude than the SBR is attributed to the gentle amplification of sHSAF_borehole, except for the troughs.Such amplification is significant in a frequency range below the fundamental peak frequency in the SBR.This is what we mentioned above as the second problem associated with sHSAF_borehole.In Fig. 24, we plot the theoretical sHSAF_surface and sHSAF_borehole spectra based on the 1D ground response analysis of the inverted Vs profile using the EHVR method (Nagashima et al. 2014), along with the assumed damping of Q = Vs/10.We can see that the observed sHSAF_surface is well reproduced by the 1D ground response analysis, whereas the theoretical sHSAF_borehole has extremely deep troughs at the cancellation frequencies.This means that the actual amplitude of the downgoing waves is significantly smaller than that of the upgoing waves, or that the phases of the reflected waves are not very coherent.In any case, borehole records should not be regarded as outcrop records of the corresponding layer with the same Vs at the borehole sensor levels.

Conclusions
In this letter, a brief history of ESG studies is first provided, and then emerging techniques, primarily based on empirical observations in Japan, are developed and described.Considering the history of ESG studies, we may conclude that the studies have been performed in a straight direction of expansion from the surface geology, to the shallow layers down to the engineering bedrock (or only 30 m from the surface), and finally throughout the basin structure down to the seismological bedrock.As has been emphasized repeatedly in this review, the use of the information for the entire structure from the surface down to the common seismological bedrock with an S-wave velocity of 3 km/s or higher is a logical consequence for any ESG study because the seismological bedrock is the only place where we should not see any site effects by definition.The use of a rock formation with an S-wave velocity higher than 760 m/s in one region, for example, may work equally well as the seismological bedrock in that region.However, the scheme using such a rock formation as a reference is not guaranteed to work in other regions.In addition, the various levels of site amplification below such soft layers will significantly contribute to the uncertainty of the prediction, and their average will be mapped onto the source spectra, particularly in the high-frequency range.
The emerging techniques introduced here all targeted the ESG from the seismological bedrock.The VACF method provides an easy way of utilizing EHVR to obtain HSAF directly.For quantitative theoretical modeling of HSAF, we use the EHVR inversion method to delineate an effective S-wave velocity structure for HSAF of the S-wave part.The WSR method can account for the contributions of basin-induced surface waves in a broadband frequency range without performing rigorous calculations of waveforms through numerical modeling of complicated 2D/3D shallow crustal structures.The EMR method can effectively transform the MHVR into pseudo-EHVR, which can be used for velocity inversion as well as the direct estimation of HSAF using the VACF method.We can also use the EHVR for weak and strong motions to grasp the degree of nonlinearity in HSAF and empirically model the average nonlinear features that appeared only in the strong motions.Based on these outcomes and subsequent discussions, we emphasize the following five important points necessary for a stable and quantitative evaluation of ESG: 1) The use of Fourier spectra, rather than the response spectra.
2) The use of the seismological bedrock outcrop as a reference, rather than just a rock outcrop.
3) The direct use of the observed spectra in the target frequency range, rather than their scalar proxies as a subset of the entire velocity profile.4) The use of relative ratios between strong motion and weak motion, rather than strong motion spectra directly in nonlinearity evaluations.5) The use of extracted outcrop motions or spectra, rather than the borehole motions or spectra as a reference.
Here are brief descriptions of the rationales for these basic rules: 1) As clearly shown in Fig. 20, the response spectra cannot capture the decreasing amplitudes of Fourier spectra in both lower and higher frequency ranges and the effects of long durations for the whole duration of motion.2) As clearly shown in Fig. 21, a potential error will be surfaced with the use of soft rock sites, instead of the seismological bedrock outcrop as a reference.3) Multiple proxy approaches are becoming increasingly complicated and approaching the use of the entire velocity profile at a site, which apparently supports the optimal result with the use of either the theoretical or empirical HSAF that reflects the entire velocity structure underneath.4) Deviations (shifts) from the linear regime due to the stronger inputs should be explored because they are the only physical entity representing the nonlinearity in the HSAF.5) As clearly shown in Fig. 24, the spectral characteristics of borehole records are contaminated by both reflected body waves from the surface and amplifications of layers below borehole sensors.
As for rule No. 3, we would like to refer to the recent papers of Zhu et al. (2022) and Zhu et al. (2023), although their comparisons were not based on the direct regression as in the ordinary GMPEs but based on the AI technique within a framework of the same HSAFs obtained by GIT (Nakano et al. 2015).Therefore, their results might be biased toward better evaluation for site proxies than the direct regression on the GMPE for site proxies.

Fig. 2
Fig. 2 Theoretical 1D two-layered model on the top left, derived theoretical site amplification with respect to the normalized frequency on the bottom left, and deconvolution analysis at two sites with different damping factors for the alluvial site, Shinobugaoka on the top right, based on the observed weak motions on the bottom right (after Takahashi and Hirano 1941)

Fig. 3
Fig. 3 Total variability of PGA from GMPEs developed in 40 years (after Strasser et al. 2009)

Fig. 6
Fig. 6 Comparisons of the empirical site amplification factors: a HSAF of the surface sensors relative to the seismological bedrock (HSAF_surface), b HSAF of the borehole sensors to the seismological bedrock (HSAF_borehole), and c the ratio between them (SBR), at 16 KiK-net sites in Aomori Prefecture and the overall geometrical means of all the 696 sites used (red broken lines)

Fig. 8
Fig. 8 Observed seismograms during the foreshock of the 1930 Kita-Izu earthquake of M7.3 at the University of Tokyo (after Imamura 1931).Mainshock records were all off-scale.We can see a very prominent phase of surface waves in the later part (from about 4 times the S-P time), especially in the NS component

Fig. 9
Fig.9Time-domain responses on the surface of the double trapezoidal basin for a vertically incident SH wave with a period of 4 s, which corresponds to the fundamental resonant period of the basin (afterKawase and Aki 1989)

Fig. 11
Fig. 11Obtained correction function for S-wave, sVACF, from all the sites with more than 10 events, together with the average plus/minus one standard deviation (left, Panel a) and comparisons of the simulated sHSAF (red broken line) and observed sHSAF (black line) at four example sites (right, Panel b).sEHVR used for correction is also shown with blue lines (afterIto et al. 2020) , the matching with the unified model from the seismological bedrock to the surface reported by Senna et al. (2014), Senna et al. (2019), and Wakai et al. (2019) with a spatial resolution of 250 m in the Kanto and Tokai regions is suitable for reproducing the empirical sHSAF derived from the GIT (Nakano

Fig. 12
Fig.12Parametric study on the effect of the bottommost layer on the theoretical sEHVR for the best-fit model at MYG006 with 14 layers byNagashima et al. (2014).When we omit the two layers in one step from the original velocity model (right), not only the amplitude of the lowest frequency peak but also the amplitudes of peaks and troughs in the higher frequency range are strongly affected (afterKawase et al. 2018)

Fig. 13
Fig. 13 Contour maps of FMR (left) and AMR (right) after the interpolation with the spatial resolution of 250 m byGMT (after Ito et al. 2021a).Note that FMR tends to have a stronger spatial variability with a smaller coherence distance, while AMR tends to have smaller spatial variability with a larger coherence distance.Larger AMR areas tend to be located inside large basins, while the larger FMR areas are located both inside and outside of the basins ) support the (11) wHSAF = wEHVR • wVSAF* = wEHVR • sVACF • WSRv, similarity of WSRh and WSRv.As mentioned by Kawase et al. (2021), Eq. (

Fig. 15
Fig. 15WSRh (top)  and WSRv (bottom) maps for the whole Japan at 0.5 Hz and 1.0 Hz based on the GIT byNakano (2021).Red dotted curves are the approximate location of the Itoigawa-Shizuoka tectonic line.Interpolation is done by GMT(Smith and Wessel 1990)

Fig. 22
Fig. 22 The grid search result for the matching function as the function of α and β on the left (Panel a) and on the right (Panel b) the observed EHVR in the linear (black) and nonlinear regime (purple), together with the reproduced EHVR (green) with the highest matching parameters on the left (Panel a)

Fig. 23
Fig. 23 Binned average (〇) of α and β for all the high PGA data (・), together with the regression curves for α and β (red)